# Numerical simulation of finite dimensional multibody nonsmooth mechanical systems

Bernard Brogliato, Tony ten Dam, Laetitia Paoli, Frank Génot, Michel Abadie

► To cite this version:

Bernard Brogliato, Tony ten Dam, Laetitia Paoli, Frank Génot, Michel Abadie. Numerical simulation of finite dimensional multibody nonsmooth mechanical systems. Applied Mechanics Reviews, 2002, 55 (2), pp.107-150. 10.1115/1.1454112. hal-01349852

HAL Id: hal-01349852

https://hal.science/hal-01349852

Submitted on 29 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical simulation of finite dimensional
multibody nonsmooth mechanical systems

## B Brogliato

INRIA Rhône-Alpes, ZIRST Monbonnot, 655 avenue de l'Europe, 38334 Saint Ismier, France; Bernard.Brogliato@inrialpes.fr

## AA ten Dam

National Aerospace Laboratory NLR, Department of Mathematical Models and Methods, Information and Communication Technology
Division, PO Box 90502, 1006 Amsterdam, The Netherlands; tendam@nlr.nl

## L Paoli

Université Jean Monnet-St Etienne, Equipe d'Analyse Numérique UMR CNRS 5585, 23 rue du Dr Michelon, 42023 Saint-Etienne
cedex 2, France; paoli@anumsun1.univ-st.etienne.fr

## F Génot

INRIA Rocquencourt, BP 105, 78153 Le Chesnay cedex, France; Frank.Genot@inria.fr

## M Abadie

Schneider Electric, Avenue des Jeux Olympiques, S2E, bureau 204 38000 Grenoble, France

This review article focuses on the problems related to numerical simulation of finite dimensional nonsmooth multibody mechanical systems. The rigid body dynamical case is examined here. This class of systems involves complementarity conditions and impact phenomena, which make its study and numerical analysis a difficult problem that cannot be solved by relying on known Ordinary Differential Equation (ODE) or Differential Algebraic Equation (DAE) integrators only. The main techniques, mathematical tools, and existing algorithms are reviewed.

## 1 INTRODUCTION

This article focuses on the problem of numerical simulation of systems composed of interacting rigid bodies. It aims at showing the state of the art in this field and, at the same time, introduces the main features, difficulties, and proposed solutions to simulation of nonsmooth multibody dynamics. This topic is an important part of the multibody systems dynamics general problem, which is itself a major field of investigation [1–3]. It concerns many domains of applications: interactive computer graphics, virtual prototyping [4,5] (electrical devices, watch mechanisms design), aerospace (control of space manipulators, liquid slosh phenomena in satellites), diesel engines, landing gear, analysis of vibrations in machines [6,7] assembly and disassembly processes [8], design of optimal protection from collisions [9], granular matter—dynamics of sandpiles, gravels, planetary rings—[10–12],

dynamics of buildings and structures [13], modeling and control of buildings, hopping [14], walking machines [15], kinematic chains with backlash [16,17], manipulators performing complete robotic tasks (deburring, grinding, polishing, chamfering) [18–20], etc. It also involves various scientific communities: Mechanical Engineering, Systems and Control, Applied Mathematics, Robotics, and Physics. In real-world applications, there exists a trend to use one and the same simulator for design, engineering, and training purposes [21]. A typical example is the simulation environment for a space manipulator [22]. The challenge for model developers is to provide users of simulators with accurate and realistic responses of simulated nonsmooth systems with sufficient level of detail [23]. Issues like repeatability and uniqueness of solutions are clearly of paramount importance. If real-time simulations are desired, the numerical algorithm must also perform real-time (ie, it has to calculate the motion faster than the real process evolution), possibly in a trade-off
with accuracy requirements. The choice of model, eg, rigid body versus compliant, becomes important. We will deal with these issues in the remainder of this paper. With the growing availability of commercial and research code for nonsmooth mechanical systems, the capabilities of these codes become of interest, especially since those packages tend to limit the time needed to code an application model.

The choice between rigid body or compliant models is a problem in itself. At this stage, the reader should understand these two definitions in a very rough fashion since a rigorous definition of rigidity and compliance requires some care [24]. This will be discussed in this article. Let us just point out that the choice of the model in general strongly depends on one's goals and field of application. It may be the result of an iterative process with successive simplifications or complications. There are applications (eg, Systems and Control, Robotics) in which one needs simple enough models (eg, that allow the designer to construct a stable feedback controller). At the same time, such models have to predict the real motion reasonably well, hence keep enough physical information inside them. In this setting, rigid body models seem quite suitable in many applications, where a compromise between accuracy and tractability has to be found.

As is well known, the principal qualities of a numerical algorithm are its accuracy and its stability (conditional stability, A-stability, or unconditional stability). Accuracy is related to the order of the scheme. Stability is related to the boundedness of the numerical solution for fixed steps of integration, as time progresses. A third property is the convergence of the numerical solution towards some function, on a fixed time interval, as the integration step goes to zero. One question that comes to mind is: are these properties transported from ODEs (free-motion systems) or DAEs (bilaterally-constrained systems) towards unilaterally constrained systems? For example, does a multistep algorithm that is known to converge for index 2 DAEs ([25], Ch. VII.3), still converge for nonsmooth systems? As we shall expose throughout this paper, the answer is not trivial, and might be negative in many cases. For instance, the Euler method is known to be not very accurate for ODEs or DAEs ([26], page 247). However, the occurrence of topology modifications and/or impacts (state re-initializations) may also render the higher-order methods useless, or at least much less powerful [27–29]. Moreover, some widely used tricks like Baumgarte's stabilization of constraints become quite inadequate in the case of inequality constraints, where the sign of the Lagrange multipliers is of primary importance. To summarize, the challenge could be formulated as follows: Consider the multibody system in [25] Ch. VII.7 (a seven-body mechanism), on which six different DAE codes are compared. Add some unilateral constraints with friction (say 15, which is a reasonable number). How much is the problem complicated? What does there remain to be done before getting a reliable and accurate simulation software, starting from the proposed algorithm (DASSL, DOPRI, RADAU, etc)? The answer to both questions is: a lot! We will lay the foundation for this answer in the remainder of this paper.

In the following, we shall generically denote the impact times as $t_k$, and the set of admissible positions as

$$K=\{q: \phi(q) \ge 0\} \quad (1)$$

with $\phi(\cdot)$ a differentiable vector valued function. The general dynamics of the systems we deal with may be written as follows:

$$ \left\{ \begin{array}{l} M(q)\ddot{q} = Q(q, \dot{q}, t) + \nabla \phi(q)\lambda + P(q, \dot{q}, \lambda) + \nabla \phi_b(q)\lambda_b \\ \phi(q) \ge 0, \quad \lambda \ge 0, \quad \phi^T \lambda = 0 \\ \phi_b(q) = 0 \\ \text{State re-initialization (Impact rule)} \\ \text{Coulomb friction.} \end{array} \right. \quad (2) $$

The matrix $M(q) \in \mathbb{R}^{n \times n}$ is the symmetric positive definite inertia matrix, $Q(q, \dot{q}, t)$ accounts for inertial (centripetal, Coriolis), external (eg, control inputs), and frictional (Rayleigh dissipation) generalized torques, $P(q, \dot{q}, \lambda)$ accounts for frictional (Coulomb) generalized forces, $\phi_b(q) \in \mathbb{R}^{m_b}$ is a set of bilateral constraints, $\lambda_b$ is the corresponding Lagrange multiplier vector, $\phi(q) \ge 0$, $\phi(q) \in \mathbb{R}^m$ is the set of unilateral constraints which indicate in which domain of the configuration space of the system $q$ is allowed to evolve, and $\lambda$ is a vector of Lagrange multipliers which represent the normal part of the interaction between the system and its environment. As a convention in this paper, $\nabla$ will always denote the gradient in the Euclidean metric (ie, $\nabla \phi^T = \partial \phi / \partial q$), though we shall make use of other metrics. The complementarity (or orthogonality) between $\phi(q)$ and $\lambda$ indicates that physical interactions producing forces without contact (magnetic effects) and gluing forces are not taken into account in the model. In the following we shall generally consider that there are no permanent bilateral constraints in order to focus our attention on inequality constraints only. Systems as in (2) are sometimes named *complementary slackness* mechanical systems by analogy with the mathematical programming and convex analysis language ([30] page 226). Assume frictionless constraints, and that $\lambda$ is an impulsive term, ie, $\lambda = p_k \delta_{t_k}$ with $\delta_{t_k}$ the Dirac measure at time $t_k$. Let $f(t_k^+)$ and $f(t_k^-)$ denote the right and left limits of $f(\cdot)$ at $t=t_k$, respectively. Then the Lagrange equations at $t_k$ become [18]:

$$M(q(t_k))[\dot{q}(t_k^+) - \dot{q}(t_k^-)] = \nabla \phi(q(t_k))p_k. \quad (3)$$

The mathematical problems associated to the dynamics in (2) are not trivial. As we shall explain later in this paper, the dynamics in Eqs. (2) and (3) represent a Measure Differential Inclusion (MDI). Its well-posedness (existence and uniqueness of solutions) is still an open problem in the general case. One of the very first contributions to this field can be found in [31]. Other contributions have been made by Monteiro-Marques [32], Stewart [33,34], Mabrouk [35,36], Carriero and Pascali [37,38], Heemels et al [39,40], van der Schaft and Schumacher [41,42], Lötstedt [43], Percivale and Buttazzo [44–48], Ballard [49,50], and Schatzman and Paoli [51–57]. In particular, systems with multiple contacts and/or friction create deep modeling and analysis difficulties. It is not the goal of this survey to provide many details on the mathematical aspects of the continuous-time dynamics. However, it is worth being aware of them since it seems difficult, even hazardous, to simulate a system reliably that is not well-posed (see Section 6.2)! In particular, one should be
aware of the fact that multiple impacts generally imply dis-
continuity of the solutions with respect to the initial data
[18,32,39,58]. This may have consequences on the numerical
as well as the control aspects. A non well-posed problem has
little chance to yield a well-posed numerical algorithm when
it is time-discretized. In [59], it is shown that a modified
Coulomb’s law yields non-uniqueness of solutions in a 1D
system. It is also shown that the discretization with an im-
plicit Euler scheme leads to difference equations that possess
several solutions.

It is noteworthy that the tools and models described in this
paper concern essentially mechanical systems subject to uni-
lateral constraints with or without dry friction. However,
they may also be suitable to the simulation of other non-
smooth physical systems, like electrical networks [27,42].
Indeed, some components like diodes can be modeled with
the complementarity conditions which we describe later, and
the state may possess discontinuities in complex electrical
networks, requiring reinitialization rules similar to the me-
chanical collisions [39]. As we shall see, all these nonsmooth
dynamical systems can be recast in the framework of Mea-
sure Differential Inclusions (MDIs). Infinite-dimensional as
well as quasi-static or static problems with unilateral con-
tacts are not examined in this paper, see, eg, [60] for a sur-
vey.

In this survey, we shall first recall in Section 2 the basics
of modeling nonsmooth systems so that the paper contains
some introductory material, before passing to more specific
developments for numerical simulations. In Section 3, we
introduce some peculiarities of nonsmooth multibody sys-
tems simulations when compared to simulation of smooth
systems. Such an introduction is thought to be necessary in
order to fix and clarify the fundamental discrepancies be-
tween systems as in (2) above, and various other types of
switching, discontinuous, or impulsive dynamical systems
that are often encountered in the literature (Applied Math-
ematics, Systems and Control, Computer Science). Section 4
is concerned with a discussion on the use of rigid and com-
pliant models for simulation of contact of nonsmooth
systems.

The main focus of the paper is on numerical analysis in
Sections 5, 6, 7, and 9. In Section 5, a brief overview is
given of numerical issues for bilaterally constrained me-
chanical systems. This analysis is used in Section 6 to dis-
cuss algorithms for unilaterally constrained mechanical sys-
tems. The numerical methods that are discussed are either
so-called event-driven or time-stepping methods. A compari-
son of the methods is made. Convex analysis and mathemati-
cal programming tools are at the core of many works in this
field and are therefore recalled. In Section 8, we briefly re-
turn to mathematical issues in modeling nonsmooth systems.
An overview of commercial packages capable of simulating
nonsmooth mechanical systems is presented in Section 9.
Finally, concluding remarks can be found in Section 10.

## 2 DYNAMICS OF BILATERALLY CONSTRAINED MECHANICAL SYSTEMS—DAES

2.1 Formulation of the unconstrained
continuous dynamics

When a system contains a number of subsystems, it may not
be practical to obtain a single system of differential equa-
tions that describes the behavior of the dynamical system.
The subsystems may be of a different dynamical character.
Take for example a satellite that contains liquid (cooling or
fuel). Both the motion of the liquid and of the satellite can be
described by differential equations. Another example is a ro-
botic arm that contains flexible as well as rigid bodies. So, at
first, a system description will contain submodels that cap-
ture specific (local) dynamics that can already contain opera-
tional constraint descriptions. By adding environmental con-
straints, one obtains a closed-chain description. In addition,
feedback control laws lead to closed-loop, closed-chain sys-
tem descriptions. There are two approaches to formulate
models of unconstrained continuous dynamical systems:
those that compute the mass matrix (Lagrange or Hamilton
equations, recursive Newton-Euler), and those that do not
(Featherstone recursive algorithm). Both approaches yield
formulations that can be applied to open (serial or tree)
chains and to closed-loop chains [26]. For closed-loop chains
the trick is to cut the loop, introduce a multiplier, and then
work with the multiplier method or reduce to independent
coordinates. Recursive methods seem to have lost interest in
recent years ([26] Remark, page 273). One of the reasons
may be the availability of increasing computer power at ever
lower prices. The advent of so-called O(N) methods, where
N stands for the number of bodies in a dynamics model,
happened at a time when computers were not as powerful as
today. In many applications, special effort was paid to saving
computational effort, and also to parallelization of algo-
rithms. For large N, algorithms of O(N) became of interest
in view of there efficiency. Some theoretical analyses sug-
gested that the speed of performance of an O(N) algorithm
would become superior over that of a O(N³) algorithm for a
sufficiently large number of bodies. Featherstone [61] is
O(N) and supersedes O(N³) for more than 10 bodies. In
[22], it is reported that a special O(N) algorithm is faster,
already, for two bodies. This may be due to tailoring the
algorithm to the application at hand: a robotic manipulator
with six rigid bodies and two flexible bodies. We refer to
[26] for a summary of the different methods. A note on Jour-
dain’s principle, and its relation to D’Alembert’s principle
and the so-called Kane’s dynamical equations can be found
in [22], Appendix A.

2.2 Reduction into independent coordinates

Since the presence of (equality) constraints leads to limiting
the degrees of freedom of a system (model), several authors
have proposed to find and select independent coordinates and
reduce the number of dynamical equations. See, for example,
McClamroch and Wang’s (see [18]), Kane’s methods (both
more or less based on implicit functions), Wehage and Haug
[62]. Haug’s method does not aim at reducing the dynamics
in closed-form, but just reorders the generalized coordinates in view of correcting the drift away from the constraints [5]. The choice of (dependent) coordinates is not unique, and the reduction to independent ones is neither. In general, the reduction is a difficult task [63]. Advantages of the reduction are that the reduced-order dynamics usually require less computation time and the formulation is free of drift from the constraint surface. Reduced-order techniques are often available in commercial code-generation packages (see Section 9). A drawback is that addition of a new constraint again leads to a differential/algebraic description, and the interesting theoretical issue arises whether one should obtain a new reduced-order formulation from the original description or from the already obtained reduced-order formulation.

A decomposition into submodels simplifies modeling and makes the overall system model more flexible and easy to adapt. This modular approach to modeling facilitates the possibility of exchanging models for subsystems, which is a necessity in system design. Therefore, in the remainder of this paper we will concentrate on the Lagrange multiplier approach to modeling.

Some softwares [64] include heuristic approaches which allow one to decrease significantly, in certain cases, the number of algebraic conditions. An example of a gearbox with 11 clutches and algebraic (bilateral) constraints including 212 unknowns is presented in [64]. The algorithm reduces the number of variables to 23 unknowns. This is an NP-complete<fnref n="1" /> problem. The implementation of such fast methods that allow one to eliminate useless constraints is also employed in collision detection, see Section 6.5. Clearly, a software must incorporate such algorithms to be efficient and fast enough.

## 2.3 Lagrange multipliers methods

Lagrange multiplier formulations of constrained mechanical systems can be obtained via so-called first-principles modeling. It is of interest to note that for constrained mechanical systems, Lagrange formulations can be derived that are particularly useful for simulation. See, for example, [25,65] (proposition 9.4.1) [66], and [67].

Let us view system (2) in the absence of friction, where we assume that $Q$ does not depend explicitly on time.

$$ \begin{cases} M(q)\ddot{q} = Q(q, \dot{q}) + \nabla \phi(q)\lambda \\ \phi(q) \ge 0, \quad \lambda \ge 0, \quad \phi^T \lambda = 0. \end{cases} \quad (4) $$

For systems as in (4), information on where and how contact with the constraints can be made and release can take place can be derived off-line using the algorithms in [68], for the linear case, and in [65] (Ch 7) for the nonlinear case. For simulation studies, it is useful to first look at (4) in case of equality constraints. The bilaterally constrained mechanical system in Eq. (5) is obtained.

$$ \begin{aligned} M(q)\ddot{q} &= Q(q, \dot{q}) + \nabla \phi(q)\lambda \\ \phi(q) &= 0. \end{aligned} \quad (5) $$

The Lagrange multiplier now ensures that the motion of the mechanical system satisfies the constraint equation. Differentiating the constraint equation gives $\nabla \phi^T(q)\dot{q} = 0$. Denote $(x_1^T, x_2^T)^T := (q^T, \dot{q}^T)^T$. Then starting from Eq. (5), the following first-order formulation can be obtained.

$$ \left\{ \begin{aligned} \dot{x}_1 &= x_2 \\ \dot{x}_2 &= M^{-1}(x_1)Q(x_1, x_2) + M^{-1}(x_1)\nabla\phi(x_1)\lambda \\ \phi(x_1) &= 0 \\ \nabla\phi^T(x_1)x_2 &= 0. \end{aligned} \right. \quad (6) $$

In [65] proposition 9.4.1 (see also [29] Ch. 2.5.3), it is shown that any trajectory of the system in (6) is also a trajectory of the system in (7) and vice versa.

$$ \left\{ \begin{aligned} \dot{x}_1 &= x_2 + \nabla \phi(x_1) \mu \\ \dot{x}_2 &= M^{-1}(x_1) Q(x_1, x_2) + M^{-1}(x_1) \nabla \phi(x_1) \lambda \\ \phi(x_1) &= 0 \\ \nabla \phi^T(x_1) x_2 &= 0. \end{aligned} \right. \quad (7) $$

The formulation in (7) is similar to the formulation in (6), apart from the presence of the Lagrange multiplier $\mu$. It allows one to reduce the index of the DAE from 3 to 2 ([25], Ch VII.1 [67]). The actual analysis that leads to the formulation in (7) is beyond the scope of the present paper, and we refer to [65], also for other equivalent (overdetermined) formulations. It must be emphasized that the Lagrange multiplier $\mu$ has no physical meaning. From the analysis in [29], [65], and [66], it actually follows that $\mu=0$, but the importance of formulation (7) lies in its use to obtain numerical solutions for bilaterally and unilaterally constrained mechanical systems. The two continuous-time representations are equivalent. However, they do not result in equivalent discrete-time formulations (see also Section 5.3).

All formulations (5), (6), and (7) can be used as starting points for analysis and simulation of constrained mechanical systems. In Section 6.3, these formulations will be used in what we will call event-driven simulation schemes, where decisions whether or not contact or release of a constraint takes place must be based on physical interpretations of $\lambda$ and information on the state of the system, and not only on the position.

# 3 DYNAMICS OF UNILATERALLY CONSTRAINED MECHANICAL SYSTEMS—MDIs

## 3.1 The major problem of multiple contacts—LCPs

In multibody systems with multiple contacts, the major difficulty is that the change in one contact generally implies changes at the other contacts, and this is true for detachment conditions as well as for impacts. How to treat this within the rigid body approach? This is the subject of the next sections. In this section, we will describe the effective formalisms that are used in view of numerical simulations. Since there generally exist several manners to formulate the same problem, it is quite important to choose the most suitable one (eg, the

<footnote n="1">¹See eg. Computer Theory FAQ at http://www.cs.unb.ca/lopez-o/comp-faq/faq.html for an introduction.</footnote>
simplest one if CPU time is an issue). In Section 6.6, we shall focus on the numerical tools used to solve these problems.

**Detachment conditions (Delassus' problem):** In the case of a single frictionless contact point, the detachment condition can be watched by looking at the sign of the Lagrange multiplier: if $\lambda=0$, then a sufficient condition for subsequent detachment is that the normal acceleration be >0. Things complicate when several contact points persist on a nonzero time interval. It is well-known since Delassus [18,69] (see also [7] for a simpler example) that a priori assuming (in an arbitrary way) certain values of the acceleration (for instance, one decides that the contacts persist) may lead to a contradiction because the calculated multipliers $\lambda_i<0$. Therefore, one is led to test all the possible combinations—there are $2^m$ combinations for $m$ active contacts. Such enumeration is cumbersome! However, it is a simple matter to use the constrained Lagrange dynamics in order to express the acceleration $\dot{\phi}(q(t))$ in terms of $\lambda$, see Eqs. (10) and (35) below. Doing this way and noticing that, provided the $m$ contacts have been active on the interval $[t - \epsilon, t)$, which allows one to state that the complementarity conditions in (2) imply the same conditions with $\phi$ being replaced by $\dot{\phi}$ (see (10)), one is able to construct a Linear Complementarity Problem or LCP, of the form:

$$A\lambda + B \ge 0, \quad \lambda \ge 0, \quad \lambda^T(A\lambda + B) = 0 \qquad (8)$$

where the matrices $A$ and $B$ generally depend on $q$, and $q, \dot{q}$ respectively, and on the system's physical parameters. The unknown of the LCP in Eq. (8) is $\lambda$, and we shall denote it as LCP($\lambda$). When friction is present at the contacts and in dimension 3, one gets a Nonlinear Complementarity Problem (NCP, or NLCP) due to the 3D friction cone, which is a much more tricky object, see Section 3.3. Moreau [70,71] was the first to formulate the contact law this way, hence greatly improving Delassus' work. In case of friction it is known that unexpected phenomena can occur, which are due to the lack of physical modeling. However, these phenomena do represent the behavior of the system when the contact stiffnesses are high, and will occur when doing a simulation. For instance, a planar disc may remain wedged in an angle $2\alpha$ if the friction coefficient at both contacts $\mu$ satisfies $\mu > \arctan(\alpha)$. This is related to the pointedness of the total friction cone $C_1+C_2$, see [33]. Also, the LCP may not possess any solution. In conclusion, let us state that the contact status management is one of the major issues of multibody systems simulation. But contrary to collisions which involve a lot of physical modeling, it is mainly a mathematical problem. As we shall recall later, LCPs possess several equivalent formulations (quadratic programs, complementarity functions) that may be used for numerical implementation.

Notice that to form LCPs that will monitor the topology or transition modifications, index 1 formulations are more convenient (but this does not imply that higher index approaches cannot be settled). Hence we focus on two cases: Velocity-impulse schemes:

$$\nabla \phi^T \dot{q} \ge 0, \quad \lambda \ge 0, \quad \lambda^T \nabla \phi^T \dot{q} = 0. \qquad (9)$$

Acceleration-force schemes:

$$\ddot{\phi}(q) = \overleftarrow{\nabla} \phi^T \dot{q} + \nabla \phi^T \ddot{q} \ge 0, \quad \lambda \ge 0, \quad \ddot{\phi}^T \lambda = 0. \qquad (10)$$

It is important to keep in mind that the dynamics in (2) with the position complementarity conditions is not equivalent to the dynamics with the velocity or the acceleration complementarity conditions. Actually, if $\phi(q(\cdot))=0$ on $I \triangleq [t - \epsilon, t)$, then the position complementarity implies the velocity complementarity on $I$. If $\phi(q(\cdot))=0$ and $\dot{\phi}(q(\cdot))=0$ on $I$, then it implies the acceleration formulation of the complementarity conditions. And so on with the higher derivatives. This yields the notion of Dynamic Complementarity Problems [41]. As we shall see the index, reduction is crucial in discretized schemes since it permits to formulate complementarity problems at each step of integration. Index problems in DAEs are major problems [72]. The reader may expect that the difficulties will be magnified when unilaterally constrained systems are considered. See also Section 7.2.4 concerning the choice between position and velocity complementarity conditions.

**Remark 1 (Hyperstatic systems):** When the system is hyperstatic (ie, $\text{rank}(\nabla\phi) > n$, implying $m > n$) it may become necessary to solve an LCP at each instant (in the numerical integration) to manage the detachment events. In case of non-hyperstatic system, like when $m=1$, the Lagrange multipliers can be computed uniquely and it suffices to solve an LCP when they approach zero. As long as they keep their signs strictly positive and no constraint is attained, the system can safely be considered as a DAE and integrated as such. Unfortunately, hyperstatic systems are by far the more commonly met category in practice! An example is the rocking block [18] when both contact points of the base stick [7].

## 3.2 Complementarity formulation of restitution mappings

A multiple impact occurs each time the system collides the boundary $\partial K$ at a singularity. If the codimension of the struck subspace is $m$, the impact is called an $m$-shock or $m$-impact. In the rocking block example, the shock is a 2-impact. Multiple impacts pose deep problems: continuity of the solutions with respect to the initial data may be lost<fnref n="2" /> (hence a high sensitivity with respect to the choice of the initial conditions), and even their formulation (ie, how to get an impact rule that maps the pre-impact velocities to the post-impact ones) is not trivial. Various approaches have been presented in the literature, some of which will be described in the sequel. We can already state that the problem of multiple contacts with or without friction is one of the major problems of nonsmooth mechanics (at the date of writing of this paper). Ivanov [58,73] argues that as soon as the codimension of the struck subspace is larger than 3, impacts should be modeled in a statistical way. This may be related to the hybrid manner of designing the software by introducing a

<footnote n="2">²This is easily seen with a 2D particle striking in an angle. See eg [51] Section 3.a.</footnote>
tableau containing the probabilities of events to occur. See Section 10 for more details. Actually, the study of multiple impact rules via the limit of penalized problems or the limit of sequences of simple impacts seems to be a hard task in general [74,75]. Collecting statistical data from experiments may be an alternative path.

It is possible to formulate the restitution mappings (or more generally the mappings that rule the re-initialization of the state vector when the system jumps from one mode to another one [42]) through a complementarity law. The following plastic impact rule:

$$ \dot{q}(t_k^+) = \arg \min_{\nabla \phi^T z \ge 0} \sigma_{\dot{q}}^T M(q) \sigma_{\dot{q}} \quad (11) $$

with $\sigma_{\dot{q}} \triangleq z - \dot{q}(t_k^*)$, has been apparently first formulated (with words) by Carnot [76]. Its link with Gauss' principle and complementarity laws has been first discovered by Moreau [70,71], and used in [27]. Actually, solving Eq. (11) is equivalent to solving the following LCP whose unknown is $\dot{q}(t_k^+)$:

$$ \nabla \phi^T(t_k) \dot{q}(t_k^+) \ge 0, \quad p(t_k) \ge 0, \\ p^T(t_k) \nabla \phi^T(t_k) \dot{q}(t_k^+) = 0 \quad (12) $$

where $p_k$ is defined in (3). This can be shown using the algebraic shock dynamics in (3), and the Kuhn-Tucker's conditions (in which case the percussion $p_k$ is to be interpreted as a slack variable). Notice that for an m-shock with $m \ge 1$, $p_k^T(t_k) \nabla \phi^T(t_k) \dot{q}(t_k^+) = 0$ is equivalent to $p_j^T(t_k) \nabla \phi_j^T(t_k) \dot{q}(t_k^+) = 0$ for each $j=1, \dots, m$, with all m components of each vector non-negative (hence the set of complementarity conditions in Eq. (12) is equivalent to m 1D complementary slackness conditions). Then it is easy to see that in the one degree-of-freedom case, this corresponds to a plastic shock $e=0$. The physical validity of these impact rules is discussed in [77], where it is argued that some simple multiple impacts may not be described by such laws. Basically, this is due to the physical fact that impulsive forces may appear at contacting points that were previously lasting and which detach after the shock. Therefore, both the post-impact velocity and percussion are nonzero at the impact time, precluding any complementarity between them. However, this rule may apply in other cases, as proved in [75] where the convergence of a penalized problem is studied. Moreover, a proper definition of the pairs of complementarity variables allows one to treat detachment of previously contacting points (on a non-zero time interval) to describe such motions as rocking [78].

Notice that one can replace $\dot{q}(t_k^+)$ in (12) by a combination of $\dot{q}(t_k^+)$ and $\dot{q}(t_k^-)$ to allow for elastic impacts, see Sections 3.4.1 and 3.4.3 (eg, replace $\dot{q}(t_k^+)$ by $\dot{q}(t_k^+) + e\dot{q}(t_k^-)$ and check that this implies that $\nabla \phi^T(t_k) \dot{q}(t_k^+) \ge 0$ whenever $\nabla \phi^T(t_k) \dot{q}(t_k^-) \le 0, e \in [0,1]$).

## 3.3 Constraints with Coulomb friction

Let us assume that there is one active constraint where dry friction acts at a contact point $A$. At $A$ one attaches a local frame in which Coulomb's law is expressed. The problem

now is how to write the generalized torque $P(q, \dot{q}, \lambda)$ in (2) in order to recover a complementary slackness formulation that includes both normal and tangential reactions [7,41,79-83]. Several particular cases have to be considered: 2D versus 3D friction cones, and all possible transitions between stick and sliding modes. Let us denote the tangential velocity components at the contact point $A$ as $V_t$, the normal one as $V_n$, so that $V = \begin{pmatrix} V_i \\ V_n \end{pmatrix}$, and the reciprocal contact force as $F = \begin{pmatrix} F_t \\ N \end{pmatrix}$. Consider $m=1$ in (2), ie, $\phi(\cdot)$ is a scalar function. If $V_n = \nabla \phi^T \dot{q}$, then $\lambda = N$. Then one has [82,83]:

$$ \left\{ \begin{array}{l} V_t = 0 \Rightarrow |F_t| \le fN \\ V_t \ne 0 \Rightarrow |F_t| = fN, \quad F_t = -\alpha V_t, \quad \alpha \ge 0. \end{array} \right. \\ \\ \left\{ \begin{array}{l} V_t \ne 0 \Rightarrow F_t = -fNV_t/|V_t| \\ V_t = 0: \quad \begin{cases} \dot{V}_t = 0 \Rightarrow |F_t| \le fN \\ \dot{V}_t \ne 0 \Rightarrow |F_t| = fN, \quad F_t = -\alpha \dot{V}_t, \quad \alpha \ge 0. \end{cases} \end{array} \right. \quad (13) $$

The first set of conditions describes dry friction at the velocity level, and the second set at the acceleration level. The latter is used to monitor the transitions from stick to sliding (rolling constraints [82,89]) in what we shall call event-driven algorithms. The reader may notice that the way to go from velocity to acceleration formulation of Coulomb's law is exactly the same as for the normal direction when one replaces position by velocity or acceleration complementarity conditions. Here we assume that there is only one contact to avoid cumbersome notations, but the formulations can be readily extended to the case of multiple contacts (see, however, Sections 7.2.2 and 7.4 for comments about the discretization of multiple contacts with friction). The next two Sections, 3.3.1 and 3.2.2, deal with how to transform the models in (13) into complementary slackness forms, and Sections 3.3.3 and 3.3.4 describe how to replace $P(q, \dot{q}, \lambda)$ by a set of complementarity conditions (linear or nonlinear) that enables one to monitor all the possible mode transitions for frictional contacting points. As we will see later, this enables one to discretize the whole set of dynamical equations and conditions and construct an LCP (or an NCP, or more generally any set of equations and/or conditions that can be solved with available algorithms) at each step, which permits us to advance the numerical algorithm to the next step.

### 3.3.1 Two-dimensional friction

In this case $V_t \in R$. The transformation of the velocity conditions in Eq. (13) uses the fact that a relay characteristic between two variables $v$ and $z$ (that may be expressed as $v = \text{sign}(z)$, with $-1 \le v \le 1$ if $z=0$), can be formulated with complementary slackness conditions as [7,79,81]:

$$ \begin{gathered} z = y_1 - y_2, \quad v = \frac{1}{2}(u_2 - u_1), \quad u_1 + u_2 = 2 \\ u_i y_i = 0, \quad u_i \ge 0, \quad y_i \ge 0, \quad i = 1, 2. \end{gathered} \quad (14) $$

The result follows by observing that $F_t = f \cdot N \cdot \text{sign}(-V_t)$. In particular, if $V_t = 0$ one gets $-1 \le F_t/f \cdot N \le 1$. The quantities $y_1$ and $y_2$ are the positive and negative part of $z(-V_t)$, respectively. This derivation can also be understood from (23) (see below) noting that the projection $D$ of
the friction cone on the tangent plane at A reduces to an interval, so that the results of convex analysis concerning normal cones of polyhedral sets (see [30], Ch 23) and complementary slackness conditions can be used directly. The same transformation can be performed for the acceleration law in Eq. (13) [82,83,89] introducing the positive and negative parts of $\dot{V}_t$.

### 3.3.2 Three-dimensional friction

This case is more difficult since $V_t \in \mathbb{R}^2$ and the friction cone is not polyhedral. Consequently (23) is still true, but cannot be transformed into complementarity conditions by introducing additional slack variables as in (14). The basic idea [33,79,80,82–88] is to approximate the cone, ie, to perform a polygonalization. Some more comments on this will be made in Section 7.4. Notice that this is the same as transforming the disc $\mathcal{D}$ in (23) below into a (convex) polyhedral set $\hat{\mathcal{D}}$ [80,83]. The pyramid approximation obviously corresponds to $\hat{\mathcal{D}}$ being a square [89]. A finer approximation corresponds to adding edges to the polygon $\hat{\mathcal{D}}$. Then a result of convex analysis ([30], Ch 23) allows one to generalize (14) and to express this approximated Coulomb’s law with complementary slackness conditions of the form (we choose the acceleration formulation in (13)):

$$ 
\begin{aligned}
-\dot{V}_t &= \sum_{i=1}^{e} \mathbf{e}_i \xi_i, \quad \sigma_i = f \cdot \mathbf{N} - \mathbf{e}_i^T F_t \\
\sigma_i &\ge 0, \quad \xi_i \ge 0, \quad \sigma_i \xi_i = 0, \quad i=1, \dots, e
\end{aligned}
\quad (15) $$

where $\mathbf{e}_i$ is a unit vector normal to the ith edge of $\hat{\mathcal{D}}$, $\mathbf{e}$ is the number of edges, and $\xi_i$ are slack variables. A similar derivation can be found in [89] for the pyramid approximation, where the variables $\sigma_i$ and $\xi_i$ are directly the positive and negative parts of the corresponding quantities, provided the local frame and $\hat{\mathcal{D}}$ are properly oriented (hence the $\mathbf{e}_i$s simplify). Notice that almost all authors choose an outer approximation of the cone, although there does not seem to be any fundamental reason for such a choice that yields an overestimation of frictional effects. Actually the most important feature is to find out a formulation that involves as few slack variables as possible, so as to simplify the subsequent numerical procedure.

Let us assume that $\dot{V}_t = J_t \ddot{q} + h(t)$. At this stage one has $P(q, \dot{q}, \lambda) = J_t^T F_t$, and Coulomb’s law is expressed through complementary slackness conditions as in Eqs. (14) or (15). The next step is: how to get an LCP that allows one to calculate the unknowns (the normal and tangential contact forces) at each time?

### 3.3.3 Sliding contacts

The velocity formulation in Eq. (13) holds, which is linear in the contact force [89]. There is a direct relationship between $F_t$ and $\mathbf{N}$, whatever the dimension. Then $P(q, \dot{q}, \lambda) = J_t^T(q) f \cdot N V_t / |V_t|$. It is therefore a simple matter to express $\ddot{q}$ as a function of $q, \ddot{q}$, and $\mathbf{N}$, then to replace $\ddot{q}$ in $\ddot{\phi}(q)$ to get a LCP as (8) with unknown $\mathbf{N}$. In the case of $m_s \le m$ sliding contacts, $\lambda^T = (\mathbf{N}^{(1)}, \dots, \mathbf{N}^{(m_s)})$ in (2). Due to friction, the matrix $A \in \mathbb{R}^{m_s \times m_s}$ (see (8)) is not symmetric in general and it may also lose its copositivity property. Then existence and uniqueness of $\lambda$ at time $t$ may no longer be assured. This

yields Painlevé problems, see Section 8.1. Notice that there is no inversion of the matrix $\nabla \phi^T M^{-1} \nabla \phi$ in this procedure.

### 3.3.4 Sticking (or rolling) contacts

In the 3D case, one has to resort to the polygonalization of the friction cone to avoid nonlinear formulations, see Section 3.2.2. Obviously in the 2D case the cone is already polygonal ($\mathcal{D}$ is an interval) so we shall consider this case as a subproblem of the 3D one once the cone has been polygonalized. Notice that the acceleration friction model is highly nonlinear in the contact force, contrary to the velocity formulation for sliding contacts. Various techniques have been employed in the literature to formulate from (2) and (14) or (15) an LCP whose unknowns are slack variables (including $\lambda$ in (2)) that allow one to monitor the contact status (detachment/ contact, sticking/sliding). The method proposed in [79] does not accommodate hyperstatic systems. Glocker [7,83] uses the formalism in Eq. (15) and the fact that $\dot{V}_t = J_t(q) \ddot{q} + h(t)$, to get at each contact point an $(\mathbf{e}+1)$-dimensional LCP in $\lambda$ plus $e$ slack variables $\xi_i$ taken in such a way that hyperstatism problems are avoided. No inversion of the matrix $J_t^T M^{-1} J_t$ is needed. A similar method is proposed in [82,89] for the pyramidal approximation, together with some existence result for the LCP. Without polygonalization, a Nonlinear Complementarity Problem (NCP) of the form

$$ y \ge 0, \quad g(y) \ge 0, \quad y^T g(y) = 0 \quad (16) $$

has to be solved for some slack variable $y$ [83], or even some non-standard NCPs (generalized complementarity problems) [82,89–91], see Section 3.4.2. NCPs are more complex to solve than their linear counterparts, see Section 6.6, and potential users have to develop their own codes. However they may involve less variables than polyhedral approximations [90]. In [83], it is shown how to get a standard 4D NCP per contact. The all-rolling contacts problem is formulated with a Quasi-Variational Inequality in [90] and results on existence of solutions are given. The work in [92] is not based on an approximation of the friction cone and treats the acceleration as the unknown, see Section 6.7.1 for more details.

**Remark 2:** From a hybrid dynamical system point of view [42] Coulomb friction adds modes or discrete-event states to the system.

### 3.3.5 Shocks with friction

The problem of frictional impacts is complex, and we shall not insist on it in this article, see [18,93] for details. In the complementary slackness framework, it is of interest to extend the frictionless rule in Eq. (12) to the frictional case. Several authors proposed complementarity formulations of frictional impacts [7,27,86,94]. The basic idea of Lötstedt [27] is to extend directly Eqs. (11) and (12) as:

$$ 
\begin{aligned}
& \min \frac{1}{2} P_k^T M^{-1} P_k + P_k^T \dot{q}(t_k^-) \\
& p_k \ge 0, \quad -f p_k \le p_{t,k} \le f p_k
\end{aligned}
\quad (17) $$

where $p_k$ is the normal percussion in (3), $p_{t,k}$ is its tangential counterpart, and $P_k$ denotes the total percussion vector (ie, the right-hand-side of (3) for frictional impacts). The second line in Eq. (17) is the direct extension of Coulomb’s law at
the percussion level, a model that may be energetically in-
consistent, see Section 7.1.2. If $f=0$ then (17) reduces to
(11) using Kuhn-Tucker conditions. There are different ways
to rewrite Eq. (17), using the dual QP and (3), see also (54).
Some details on the works in [86] will be given in Section
7.4. For 2D problems, it is similar to the one in [7], where
the collision is decomposed into a compression and restitu-
tion phases. Concerning [94] see Section 7.2.4.

3.4 Convex analysis in nonsmooth mechanics

In this section we present compact forms of the MDIs in Eqs. (2) and (3). The formulations that follow in this section are based on mathematical tools from convex analysis [30,95]. It is not our goal in this article to provide details on such mathematical tools. However, it seems mandatory to make a short presentation of these approaches, since they are used to develop well-posedness and numerical issues of many algorithms. A tool that will be needed in the following is the definition of a subdifferential. The subdifferential of a convex function $f(\cdot)$, denoted as $\partial f$, is the set of subgradients of $f$ at $y$, ie, of vectors $\gamma$ satisfying $f(x)-f(y) \ge \gamma^T(x-y)$ for all $x$. This is the generalization of the derivative for non-differentiable convex functions. Let $K$ denote a convex set (not necessarily the one in (1), though this may be the case as we shall see further). The function $\psi_K(y)$ is called the indicator function of the set $K$, and is defined by $\psi_K(y)=0$ if $y \in K$, $\psi_K(y)=+\infty$ if $y \notin K$. The values inside and outside the set $K$ are consistent with the definition of subgradients.

The subdifferential of $\psi_K(\cdot)$ generalizes the notion of nor-
mal directions to $\partial K$ at $y$. Roughly speaking, this represents
the outward normal directions to the boundary of $K$ at some
point $y$, and if the boundary is smooth this is the usual out-
ward normal half line. If the boundary is nonsmooth, it gen-
eralizes to a cone, called the normal cone and is denoted as
$N_K(\cdot)$. When $y \in \text{Int}(K)$, then $\partial \psi_K(y) = \{0\}$, and when
$y \notin K$, then $\partial \psi_K(y) = \emptyset$. One also defines tangent cones
$V(y)$ as $V(y) = \{v: \forall \ x \in N_K(y), \ x^T v \le 0\}$ when $y \in K$.
Thus if $y \in \partial K$, $V(y)$ is another convex cone. In this case
tangent cones reduce to (inward) half spaces for differen-
tiable codimension one boundaries $\partial K$, see Fig. 1a,c. If $y$
$\in \text{Int}(K)$, then $V(y) = \mathbb{R}^n$, and if $y \notin K$, then $V(y) = \emptyset$.
When $K$ is defined as in (1), one may define the tangent cone
to $K$ at $q$ as [96]:

$$
V(q) = \{v \in \mathbb{R}^n : v^T \nabla \phi_i \ge 0, \ \forall \ i \in J(q)\} \quad (18)
$$

with $J(q) = \{i \in \{1, \dots, m\} : \phi_i(q) \le 0\}$. One notes that this
definition coincides with the original definition of the tangent
cone as long as $y \in K$. The fact that $V(y)$ is given a meaning
outside $K$ is useful for numerical applications in which the
admissible domain of configurations in (1) may be violated.
In order to motivate the reader, let us note that complemen-
tarity conditions between two scalar variables $\lambda$ and $y$:

$$
\lambda \ge 0, \quad y \ge 0, \quad \lambda \cdot y = 0 \tag{19}
$$

can be equivalently formulated as

<figure><img src="image_2.png" /><figcaption>Fig. 1 Tangent and normal cones, and the corner law</figcaption></figure>
$$ -\lambda \in \partial \psi_K(y), \quad \text{with} \quad K = \mathbb{R}^+. \qquad (20) $$

This means that the corner law whose graph is depicted in Fig. 1b, is actually the graph of $-\partial\psi_{\mathbb{R}^+}(\cdot)$. Evidently the role of $\lambda$ and $y$ can be reversed in (20). This generalizes to the case of complementarity conditions as in (2) whenever $K$ is convex. In this setting, (4) is equivalent to $-M(q)\ddot{q} + Q(q, \dot{q}) \in \partial\psi_K(q(t))$. Both Eqs. (19) and (20) are also equivalent to any one of the formulations [30] page 215; [97]:

$$ \lambda \in K, \quad \forall \ z \in K, \quad \langle \lambda, z - y \rangle \ge 0 \qquad (21) $$

$$ \lambda \in K, \quad \lambda = \text{proj}_K[\lambda - \rho y], \quad \text{for any } \rho > 0. $$

The first formulation in Eq. (21) is a variational formulation of the complementarity conditions. This illustrates that there are several manners to express the same physical law, and such a flexibility is likely to be very useful for various reasons (mathematical studies, effective dynamic formalisms in view of numerical applications). Quite interestingly, the Coulomb friction model can be expressed in a similar fashion, as we already saw. This comes from the fact that by introducing slack variables [7,41] one can recast Coulomb friction into a complementarity framework, see Eqs. (14) and (15). The same holds for many systems with piecewise linear characteristics [39,42]. Some works [64] rely on the parametrization of Eq. (19) by introducing an additional slack variable $s$. These facts clearly indicate why the dynamics in (2) can be rewritten as a special type of differential inclusion, containing measures. This will give rise to Measure Differential Inclusions (MDIs), as introduced by Moreau [98].

**Remark 3:** One very good reason for expressing the dynamics in terms of convex analysis is that it yields compact and powerful mathematical formulations and that these properties are preserved in the numerical schemes derived from them. For instance, the graphs of the complementarity (or Signorini) conditions and of Coulomb's law are monotone. It has also been used to derive well-posed impact rules in [99].

### 3.4.1 Moreau's sweeping process

The theoretical setting exposed in this section has been developed by Moreau [96,98,100] (we refer to [101] for an introduction to the sweeping process, including the first order case and to [102] for a general exposition of its mathematical properties). Moreau's sweeping process is a velocity-impulse formulation, and its main interest lies in its compactness. This property has consequences for its numerical analysis, see Section 7.2. Let us present it briefly in the case of frictionless constraints. What follows does not pretend to possess any mathematical rigor, but is to be considered as introductionary material only. The starting point is:

$$ -M(q)dv + Q(q, v, t) \in \partial\psi_{V(q(t))}(v(t^+)) \qquad (22) $$

The variable $v$ satisfies $q(t) = q(0) + \int_0^t v(\tau)d\tau$, ie, it equals $\dot{q}(t)$ almost everywhere in the Lebesgue measure sense. Notice that Eq. (22) is an *implicit* formulation in $v(t^+)$. In particular, if $v$ has a jump $\sigma_v = v(t^+) - v(t^-)$ at $t$ then $dv = \sigma_v \delta_t$. When $v(\cdot)$ is time-continuous, then one just replaces its right-limit by $v(t)$ in (22). Let us intuitively explain

what the differential inclusion in (22) means. The interested reader may have a look at [18,96,100] for further details and examples with figures. The set $V(q(t))$ is the tangent cone to $K$ (see (1)) at $q(t)$. When $V(q(t)) = \mathbb{R}^+$ then the interpretation of Eq. (22) is obvious from Eqs. (19)–(21) (this is, for instance, the case for a one-dimensional bouncing ball touching the ground at time $t$). The term $\partial\psi_{V(q(t))}(z)$ is a subdifferential and can be understood as follows when $q(t) \in \partial K$. If the vector $z^3$ points inward $V(q(t))$, then $\partial\psi_{V(q(t))}(z) = \{0\}$, so that (22) reduces to the time-continuous Lagrangian equations. If $z$ belongs to the boundary of $V(q(t))$, then $\partial\psi_{V(q(t))}(z)$ is the normal cone to $V(q(t))$ at $z$. If $z$ points outward $V(q(t))$ (thus outside $K$) then $\partial\psi_{V(q(t))}(z) = \emptyset$: this situation is forbidden, and one sees that the implicit formulation in (22) implies that $v(t_k^+) \in V(q(t_k))$. Now recalling that at $t=t_k$, the variable $z$ is the post-impact velocity, one obtains that (22) becomes $-M(q(t_k))\sigma_v \in \partial\psi_{V(q(t_k))}(v(t_k^+))$. This is equivalent to a quadratic problem as the one in (11) or (12) and this implies that $v(t_k^+) \in \partial K$: this is a generalization of plastic impacts for codimension one constraints. Another interpretation is in terms of *proximation*, that is $v(t_k^+) = \text{prox}[V(q(t_k)), v(t_k^-)]$, where the proximation is to be made in the kinetic metric. $v(t_k^+)$ is therefore the closest element to $v(t_k^-)$ in $V(q(t_k))$, in the kinetic metric. If $q(t) \in \text{Int}(K)$, then $V(q) = \mathbb{R}^n$ so that $\psi_{V(q)}(z) = 0$ for all $z \in \mathbb{R}^n$ and $\partial\psi_{V(q)}(z) = \{0\}$: (22) is a classical Lagrange equation.

**Remark 4:**

* From Eqs. (19) and (20), one suspects that (22) implies complementarity conditions between the contact reaction and $v(t^+)$. This is indeed the case, see [18], [96].

* It is possible to show that the set $\partial\psi_{V(q(t))}(v^+) \subseteq N_K(q(t))$. This can be seen by drawing a 2D angle and plotting the various cones and half-spaces involved, see Fig. 1c where $\partial\psi_{V(q(t))}(v_1) \subseteq N_K(q(t))$. This means that (22) is not strictly equivalent to (4).

* It is possible to include a restitution coefficient $e$ in this formulation by replacing $v$ in both sides of (22) by an averaged velocity $u = v(t^+) + ev(t^-)/1 + e$ [35,96], see Section 7.2. It is equivalent to inclusion of $u$ in Eq. (12), as pointed out after (12) (notice that in Eq. (22), $(1+e)$ in the denominator is needed to encompass continuous velocities as well).

* The ‘-’ sign in the left-hand-side of (22) is only a consequence of some standard notations in convex analysis, see (20).

* Mechanical systems with position constraints have index 3 [25]. Indeed one needs to differentiate the constraint $\phi(q) = 0$ three times to recover a set of ODEs with state $(q, \dot{q}, \lambda)$. Now consider the simple dynamics $M\dot{q} = \lambda$, $\phi(q) = 0$. Let us write it in a velocity-impulse form, ie, $M(\dot{q} - \dot{q}(0)) = p$, with $p(t) = \int_0^t \lambda(\tau)d\tau$. Then one needs only to differentiate $\phi(q)$ twice (with respect to time) to get an ODE with state $(q, p)$. Thus the constrained modes of (22) correspond to index 2 problems.

<footnote n="3">³Usually one draws the sets $q(t) + N_K(q(t))$ instead of $N_K(q(t))$ [102], so $z$ should be understood as the vector emanating from $q(t)$. Similarly for the other sets.</footnote>
It is important to notice that the essence of the dynamics in (22) is that the complementarity conditions in (2) and in (12) are stated through the subdifferential in the right-hand side of (22), which is a function of the velocity (for fixed position). Velocity spaces are linear spaces (tangent to the configuration space at each $q$), whereas configuration spaces are seldom linear. Consequently, the use of velocity spaces allows one to recover all the advantages of linear spaces (for instance the tangent cone $V(q(t))$ is always convex, whereas $K$ in (1) may not be). Mathematicians use such Nonsmooth Analysis tools in their existence and uniqueness of solutions studies, which are of high importance for numerical analysis. Indeed, on the one hand, well-posedness results are often based on time-discretization of the continuous dynamics (similarly to the fundamental result for ODEs that uses a Euler discretized scheme whose piecewise constant solutions are shown to converge to solutions of the continuous dynamics when the vector field is continuous [103]). This is what is usually called the *time-stepping method* and will be described later in this article. On the other hand, a necessary first step (but not sufficient to get a robust simulator!) is to establish the existence and uniqueness of solutions before designing a numerical simulator.

### 3.4.2 Coulomb friction and the principle of maximal dissipation

Moreau [104] has shown that a Coulomb friction model actually satisfies a maximal dissipation principle. When there is sliding between two bodies at the contact point, it is intuitively clear that this model obeys such a principle. Indeed it says that the tangent reaction $F_t$ has to be in the opposite direction to the tangent velocity $V_t$, see Eq. (13), and the product $-F_t^T V_t$ is maximal when one seeks $F_t$ inside its admissible domain $\mathcal{D}$. Here one has $\mathcal{D}=\{R \in \Pi: |R| \le f \cdot N\}$ where $\Pi$ is the common tangent plane at the contact point. Incidentally, it is understood here that $N$ is known. One sees that $\mathcal{D}$ is the projection of a section of the friction cone $\mathcal{C}$ on $\Pi$. In particular, $\mathcal{D}$ is convex, which allows one to use the tools presented in the introduction of Section 3.4, replacing $K$ by $\mathcal{D}$. With this in mind, it is possible to show (see, eg, Eqs. (14), (19), and (20)) that Coulomb friction can be formulated as

$$-V_t \in \partial \psi_{\mathcal{D}}(F_t) \qquad (23)$$

and this inclusion permits encompassing the sticking modes as well. The reader may wish to check that the inclusion in (23) really represents Coulomb's law of friction, using the definition of the subdifferential of an indicator function (that is nothing else but the normal cone to $\mathcal{D}$ at $F_t$). It is also possible to rewrite Eq. (23) in terms of a dissipation function $\varphi(V_t) = f \cdot N \cdot |V_t|$ so that $-F_t \in \partial\varphi(V_t)$. The function $\varphi(\cdot)$ is called a superpotential or pseudopotential (by analogy with the classical definition of forces that derive from a potential $F=\nabla V$, where this time the equality is replaced by an inclusion). In generalized coordinates as used in the works described in Section 7.4, the maximum dissipation principle is formulated as $F_q = \arg \max_{z \in \mathcal{C}(q)} (-z^T \dot{q})$, where $F_q$ is the generalized contact reaction impulse and $\mathcal{C}(q)$ is the generalized

friction cone. The velocity $V^{(j)} \in \mathbb{R}^3$ at a contact point $A_j$ and the generalized velocity are related as $V^{(j)} = J^{(j)}(q) \dot{q}$ whereas from the principle of virtual work $F_q = J^{(j)}(q) F^{(j)}$, where in our notations $F^{(j)} = (\frac{F_t^{(j)}}{N^{(j)}})$. $J^{(j)}(q)$ is a Jacobian matrix.

**Remark 5:**

* More general frictional characteristics can also be written in a form as in (23). However, when the considered set $K$ is convex one can rewrite Eq. (20) as a variational inequality as in (21). If $K$ is not convex as it is the case for more complex friction models, one has to resort to hemivariational inequalities, as introduced by Panagiotopoulos [105].

* The second set of conditions in (13) with $V_t=0$ is equivalent to $-\dot{V}_t \in \partial\psi_{\mathcal{D}}(F_t)$ [83,92] and $(15) \Leftrightarrow -\dot{V}_t \in \partial\psi_{\mathcal{D}}(F_t)$. The acceleration formulation is then useful to derive the dynamics of rolling (sticking) contacts, see Section 3.3.

* As we saw above, the disc $\mathcal{D}$ can be approximated by a convex polytope $\mathcal{D}$ [80,83]. This is an idea that is used in some algorithms to be presented later, see Section 7.4. However, the formulation in Eq. (23) remains valid, even with the modified friction cone. The advantage is that one deals with LCPs instead of NCPs. The case when $\mathcal{D}$ is a square (so that the cone is approximated by a pyramid) corresponds to the set $K$ in Fig. 1a. Actually basic results from convex analysis can be used to derive complementary slackness conditions (the best way to understand this is to look at [30] (Ch 23) corollary 23.8.1 and bottom of page 226, and [83]) from the polyhedral approximation of (23) or its acceleration formulation. Some event-driven and time-stepping algorithms that we will describe later hinge on such convex analysis tools to derive suitable (from a numerical analysis point of view) formulations of contact laws.

* The maximum dissipation principle is extended in [90] to more general friction laws. Some existence results (in the quasistatic case) are provided for $(q, N, F_t)$. As we shall see in Section 7.4, some authors formulate the problem in the configuration space, and $P(q, \dot{q}, \lambda)$ in (2) is written as $\mathcal{D}(q)\beta$, with $\beta$ satisfying $\psi(\beta) \le f \cdot \lambda$. The function $\psi(\cdot)$ (that defines $\mathcal{D}$) should be convex. Then the maximum dissipation principle reads [91]

$$\beta = \arg \min_{\psi(z) \le f \cdot \lambda} [\dot{q}(t^+)^T D(q) z]. \qquad (24)$$

This can be shown to be equivalent to

$$\begin{align} 0 &\in fD(q)^T \dot{q}(t^+) + \mu \partial\psi(\beta) \nonumber \\ \mu &\ge 0, \quad f \cdot \lambda - \psi(\beta) \ge 0, \quad \mu[f \cdot \lambda - \psi(\beta)] = 0 \tag{25} \end{align}$$

where $\mu$ is a Lagrange multiplier (or a slack variable), and $\lambda$ is as in (2).

* The impact rule in Eq. (12) satisfies also, a maximum dissipation principle, see (11). Its expression within the sweeping process framework has been recalled in Section 3.4.1. It can therefore also be transformed into a (convex)
quadratic program. Notice that both (11) and (12) can be written in terms either of the contact percussion $p_k$ or of $\dot{q}(t_k^+)$, using (3).

When combined with a velocity-impulse formulation, the principle of maximum dissipation allows one to avoid some inconsistency problems during the numerical integration that are encountered in rigid body dynamics with friction, see Section 7.4. Let us note that (23) is often formulated with $V_t(t^+)$ (the right limit of $V_t(t)$), a choice motivated by numerical implementation [94].

### 3.4.3 Second-order formulations

The formulations presented below concern frictionless rigid body dynamics, ie, $P(q, \dot{q}, \lambda) = 0$ in (2). The free motion of the system is described in generalized coordinates by the ODE

$$M(q)\ddot{q} = \bar{Q}(q,p,t),$$

where $p \triangleq M(q)\dot{q}$ is the generalized momentum and $\bar{Q}(q,p,t) \triangleq Q(q, M^{-1}(q)p, t)$. It should be noted that this equation allows dissipative terms. Let us assume that $\phi(\cdot)$ is a smooth function, ie, $\partial K$ (see (1)) is a smooth codimension one surface. From (2), one gets the equation

$$M(q)\ddot{q} = \bar{Q}(q,p,t) + \mu \quad (26)$$

where $\mu$ is a measure such that

$$\mu = \nabla \phi(q) \lambda, \quad \lambda \ge 0, \quad \text{almost everywhere} \quad (27)$$

and

$$\text{supp}(\mu) \subset \{t: \phi(q(t)) = 0\}. \quad (28)$$

Roughly speaking, relation (28) means that the contact forces are switched on only when the constraints are saturated and is therefore equivalent to the complementarity conditions in (2). The measure differential equation (26) does not provide a complete description of the motion and we must add an impact law in order to describe the transmission of kinetic energy at impacts. We assume a Newton's law for generalized momenta. More precisely, whenever the constraints are saturated, the tangential part of $p$ is conserved while the normal part is reversed and multiplied by a restitution coefficient $e \in [0,1]$. Here the normal and tangential part of $p$ are taken in the sense of the local metric for the space of generalized momenta, which is defined by $M^{-1}(q)$. Let us denote by $\Pi_q p$ the projection of $p$ on the tangent plane $\Pi$ to $\partial K$ at $q$ with respect to the local metric, ie, $\Pi_q p = p - \nabla \phi^T M^{-1} p / \nabla \phi^T M^{-1} \nabla \phi / \nabla \phi$. The complete formulation is finally given by Eqs. (26), (27), (28), and

$$p(t_k^+) = \Pi_{q(t)} p(t_k^-) - e(1 - \Pi_{q(t)}) p(t_k^-) \quad \text{if } \phi(q(t)) = 0 \quad (29)$$

and

$$\text{supp}(\mu) \subset \{t: \phi(q(t)) = 0\}, \\ \mu = \lambda \nabla \phi(q), \quad \lambda \ge 0, \quad \text{almost everywhere.} \quad (30)$$

One notes that the impact law in (29) is equivalent to the ones presented in Section 3.2 ([18] claim 6.1). Let us denote

by $T(t_k^+)$ the kinetic energy at an impact time $t_k^+$ and define $|p|_{q(t)}^* = p^T M^{-1}(q)p$, the norm of $p$ in the local metric of generalized momenta. One has:

$$\begin{aligned} T(t_k^+) &= \frac{1}{2} |p(t_k^+)|_{q(t)}^{*2} \\ &= \frac{1}{2} \{ |\Pi_{q(t)} p(t+0)|_{q(t)}^{*2} + |(1 - \Pi_{q(t)}) p(t_k^+)|_{q(t)}^{*2} \} \\ &= \frac{1}{2} \{ |\Pi_{q(t)} p(t_k^-)|_{q(t)}^{*2} + e^2 |(1 - \Pi_{q(t)}) p(t_k^-)|_{q(t)}^{*2} \} \\ &\le T(t_k^-) \end{aligned} \quad (31)$$

and the equality holds if and only if $e=1$. Hence this model is energetically consistent.

In the particular case of a constant mass matrix, we can consider the generalized coordinates defined by $u = M^{1/2}q$, which makes sense since $M = M^T > 0$. Then $p = M^{1/2}\dot{u}$ and with this new coordinates the impact law reduces to

$$\begin{aligned} \dot{u}(t_k^+) &= \dot{u}_T(t_k^-) - e\dot{u}_N(t_k^-) \\ &= \text{whenever the constraints are saturated} \end{aligned} \quad (32)$$

where $\dot{u}_T$ and $\dot{u}_N$ are the tangential and normal parts of $\dot{u}$, with respect to the Euclidean metric. The set of admissible positions is then described in terms of these new coordinates by $K = \{u: \phi(M^{-1/2}u) \ge 0\}$. If $K$ is convex with respect to $u$, we can replace the relations (26)–(28) by the following differential inclusion:

$$-\ddot{u} + f(t, u, \dot{u}) \in \partial \psi_K(u) \quad (33)$$

where $\partial \psi_K$ is the subdifferential of the indicator function of $K$ and can be interpreted similarly as the right-hand-side of (22), replacing $V(q(t))$ by $K$ and the right-velocity by the position, see Fig. 1a. Therefore, Eq. (33) is a second-order differential inclusion, and the total dynamics is given by (32) and (33), or equivalently (29) and (27), (28). Notice from (33) that the inertial forces are therefore a subgradient of the indicator function $\psi_K(u)$. From (22), one sees that Moreau's idea has been to consider them as subgradients of $\psi_{V(q(t))}(v(t^+))$ instead.

# 4 RIGID BODY VERSUS COMPLIANT MODELS

## 4.1 Rigid body models: Some properties

First of all it is necessary to make it very clear what is meant by *rigid* and by *flexible*, see [24] for a discussion. By rigid one may mean any model that uses velocity-impulse formulation for the collisions (like the Darboux-Keller or Routh approaches). However, this does not necessarily imply an algebraic form of the shock dynamics as in Eq. (3). In other words, what is called the rigid body approach does not preclude the incorporation of deformation effects, though the way they appear in the dynamics may differ from the usual spring-dashpot formulations. This is the case of the Darboux-Keller's shock equations [18]. Clearly, in general, a model may mix rigid body modelling features (complementarity conditions, restitution coefficients) and compliance (defor-mation of the bodies that contact). In this section compliant models are to be understood essentially as those models that give rise to ODEs with switching vector field and constant
dimension, but not to MDIs as in Eq. (2). Studies based on continuous mechanics for the bodies modelling and contact complementarity conditions [60] therefore do not fall into the type of compliance discussed here.

**Hyperstatic systems:** In the case of hyperstatic systems (95% of all real systems) the Lagrange multipliers cannot be calculated uniquely, as is well known. However, the accelerations are uniquely determined [43,70,71]. The reader may think of a chair on four legs standing on rigid ground. Clearly, for reasonable applied external forces- ie, no strange forces that create right- accumulations of impacts like in Bressan's counter-example [106]—the acceleration of the chair is unique since it detaches if the normal force compensates for gravity, or it stays at rest if it does not. This shows that if one is not interested in computing exactly the contact forces, such inconsistency is not bothering.

**Energetical behavior:** The energetical coherence can be in some cases guaranteed by the dissipativity properties of the schemes. At least this holds for certain choices of the restitution coefficients ($e=0$ in case of sweeping process with friction [96], or same value of normal $e$ and tangential $e_t$ restitution coefficients for Anitescu et al [7,84] in case of multiple impacts). In [55], an estimation of the (discrete-time) energy is also made for all $e \in [0,1]$.

**Impact detection:** Some schemes do not require an explicit calculation of the impact times: convergence of the discretized solutions towards a solution of the continuous problem is assured so one does not need a special module of detection. Decreasing the (constant) integration step automatically guarantees the improvement of the numerical result. This is the case for the time-stepping schemes described in Sections 7.2, 7.3, and 7.4. In other words, such schemes do not require (nor provide) the accurate determination of the times $t_k$ and a specific (possibly time consuming) procedure of local computation. They will therefore be suitable when one does not desire a very accurate knowledge of the dynamics at the contact points and/or the time of contact, but instead is content with a global picture of the system (eg, detection of particular orbits or attractors).

**Restitution laws:** The restitution coefficients need not at all to be constant. They can incorporate dependence with respect to initial relative velocity, material properties, shapes, dimensions etc, which take into account the influence of global and/or local deformations and sources of kinetic energy loss during the shock. The derivation of such coefficients has been and is still the object of research works. In particular, vibrational effects and transmission of energy in multiple impacts are quite fundamental issues. Some recent results are very promising [107–109].

**LCPs:** The resolution of LCPs may create some problems, however this is a topic that supersedes mechanics [110,111] and which is the object of many research works, see Sections 3.1 and 6.6.

**Global motion analysis:** It is known that nonsmooth complementarity modeling allows one to simplify the analysis of motion and to permit the study of first-return 1D Poincaré maps, which would not have been possible with compliant models of contact [112].

One major and still largely open modelling problem is that of multiple impacts with or without friction.

## 4.2 Compliant models

**Spring-dashpot model:** A compliant model that is often used is a spring-dashpot model. The linear spring-dashpot model possesses some strange properties which are often ignored or tolerated in the literature. Let us consider the following dynamics:

$$m\ddot{q} = \begin{cases} -kq - d\dot{q} & \text{if } q<0 \\ 0 & \text{if } q \ge 0. \end{cases} \qquad (34)$$

Let us consider $m=1$, $d=3$, $k=1$, $q(0)=0$, $\dot{q}(0)<0$. The values of the parameters are chosen such that the mass $m$ crushes the spring-dashpot and there is no subsequent detachment: this is a plastic impact. Numerical calculations yield $F=-kq-d\dot{q}=6.86 \exp(-2.62t) - 0.14 \exp(-0.38t)$ and this contact force can take negative values which would correspond to some “gluing” effects (which obviously are not the goal of this model). However the impulse of $F(t)$ is positive, which explains why this model provides a limit (as $k$ and $d$ diverge to $+\infty$ in a certain way) that is compatible with rigid body collision modeling [18,53].

**Compliance location:** The first point is about which sort of compliance is to be modeled: local deformation at the contact point? Or global deformations distributed in the bodies or in the kinematic chain (finite-element-like models [60,109,113–116], or elasticity at joints [117] for instance) to model the vibrational effects? Or both? A discussion on the relationships between compliant models and rigid body models can be found in [24]. In [118], the domains of application of various rigid body (also called stereocontact) models and compliant contact/impact models (Newton's, Poisson's restitution, Hertz, Saint Venant's impacts) are discussed from their respective time-scales. It is argued that multi-rate/multi-method integration codes provide better results. A simulator has been implemented that uses Gear's multi-rate integration [119]. This type of arguments is consistent with [120] who proposes to mix multi-rigid body with finite elements methods, and also with the general message at the end of this survey.

**Impact detection:** Notice that the impact detection also exists in compliant models since one has to compute the position at which contact is made so as to switch to the new vector field. In case of stiff equations and too large integration steps the energetical behavior has to be carefully checked, because long run simulations may lead to wrong results (energy loss or energy gains).

**Integration step:** The integration step $h$ has to be chosen satisfying $h \ll O(1/\sqrt{k})$ where $k$ is the contact stiffness since the time of penetration (say, of deformation) is itself $O(1/\sqrt{k})$ [53]. If one wants, for instance, to calculate 100 points during the shock and with $k=10^{10}$ N/m (a physical value in many instances like gears, pinions, systems with joint clearance [121]), then $h=10^{-7}$ s is required. If some real-time applications are required or if $n$ is very large (it can be $\ge 10^4$ in granular matter, think of a mere sandpile), then this may be an obstacle. Notice that some authors recom-
mend the calculation of at least 1000 points for each colli-
sion to assure a good accuracy of long run simulations [10].
Clearly, real-time applications—or even fast enough numeri-
cal tests—are impossible with such integration steps, see ex-
ample 1 below for numerical values. Another important point
is the fact that, in practice, impact times may be quite close
one to another [107], [[54], Ch 4], due for instance to micro-
collisions. If one desires to approximate correctly such high
frequency impacts, one is led to choose a sufficiently large k
in order to numerically detect enough impacts. Such a choice
for k is often difficult.

**Example 1:** Numerical results are reported in [121] who chooses a Hertz contact model with stiffness $k=6.58 \times 10^{10}$ N/m$^{1.5}$ (let us recall that such contact model supposes low approach velocities and is valid for central impacts only). A slider crank mechanism with joint clearance is simulated. A predictor-corrector routine is used based on Adams-PECE formulas. The simulation is performed over 30 s and its duration is 3.2 hours. A Baumgarte technique (see Section 5.2) is used to stabilize the kinematic constraints. The results are not compared to any experimental data, and they do not concern the long-term dynamics: how does the scheme behave on longer simulation times? Is such an approach realistic for long run prediction and for more complicated systems with several clearances and other unilateral contacts, both from the simulation length and the outcome accuracy? (Let us recall that certain mechanical systems like simple circuit breakers may include from 15 to 30 unilateral contacts, and from 7 to 20 bars: the identification of contact parameters and the sensitivity of a long run simulation outcome with respect to their numerical value may be a big problem). The simulation duration may be a serious obstacle in certain applications. If one wants to simulate the system on several minutes—and this may be needed in some instances—then the simulation length becomes a real problem.

Implicit algorithms may be preferred to explicit ones as it
is known that they often provide better accuracy and stability
for stiff problems. When an explicit code encounters stiff-
ness, the integration step needs to be decreased to keep sta-
bility ([25] page 21). Implicit schemes have the tendency to
filter out the high frequencies and therefore treat stiff ODEs
as DAEs. Larger integration steps are allowed for the same
accuracy. However, real-time applications may require fast
enough algorithms that preclude the use of implicit discreti-
zations (consequently of iterative algorithms).

**Physical parameters identification:** In practice, the identification of the physical parameters (damping, stiffness, or any other coefficient that appears in the contact model) may be quite difficult. Furthermore it should be recalled that some well known models like Hertz hinge on a particular stress law and are valid only for central impacts, for a simple shock between two bodies (ie, the line joining centers of gravity coincides with the normal direction at the contact). Collisions may be far from central in some applications, like colliding kinematic chains. Moreover, the influence of the parameters numerical values on the long run simulation outcome may be important. In certain simple cases of multiple

shocks (such as Newton's cradle), it can be shown that the
collision outcome is quite sensitive to such parameters val-
ues [18].

**Remark 6 (Coulomb friction):** The same type of discussion could be done about the friction model. All the algorithms described in this paper and which incorporate frictional effects use the Coulomb friction model. On one hand it is already difficult enough to properly identify and discretize this model, so that there is no need, from this point of view, to complicate things more. And, as we pointed out in Section 3.2, the main problem that we face is not the complication of the dry friction model for one contact, but its formulation and time-discretization for multiple contacts. On the other hand, this model, though simple, incorporates already a lot of physical phenomena and proves to be sufficient in many cases. In view of the large spectrum of mechanical systems that have been simulated with the simplest model, ranging from granular matter to buildings made of blocks and deep drawing processes, there seems to be little need to use a more complex model. In particular, the sticking modes are quite important in many applications (... the physics points to real discontinuities, and there is little advantage numerically in smoothing the discontinuity. The discontinuity is here to stay [91]). The simulation of a pile of rigid blocks will usually not require anything much more complex than Coulomb's law, because the phenomena other than sticking and sliding do not play a major role in the motion. However, it will necessarily involve complementary slackness conditions for both normal and tangential directions, and consequently the need for complementarity problem solvers.

Any more sophisticated model than Coulomb's law that is
to be incorporated in the framework of the developments in
this paper, should at least satisfy the principle of maximum
dissipation, both for physical reasons (what is the physical
validity of a model that is not dissipative with respect to the
supply rate $V_t^T F_t$?), and for mathematical programming rea-
sons (underlying complementary slackness conditions and
variational formulations). Models of this sort have been stud-
ied in [122,123]. More discussions on friction models in non-
smooth multibody problems can be found in [91].

4.3 Conclusions

The principal drawbacks of the rigid body approach are in-
determinate and inconsistent configurations, plus the need
for discovery of multiple impact rules. First of all, since the
rigid body model has to be seen as a limit model (in a sense
it is the model that contains the least physical information on
the process) these peculiarities are not surprising: the less
information, the worse prediction. However, the other ap-
proaches are not perfect, either, and possess their own draw-
backs. In [57,124], calculations show that elastic impact of a
particle in an angle with a penalization, implies essentially
unpredictable results after two reflections. Moreover, physi-
cal details can be introduced in a rigid body model, through
suitable restitution coefficients that may incorporate vibra-
tional effects (global deformations) like in [107], local prop-
erties at the contact point, etc. See [18] and references
therein. A lot depends on one's goals and on the domain of
application which most often requires finding a compromise
between the model complexity and its tractability: eg, de-
signing a stable feedback control algorithm or simulating the
long-run motion of a complex kinematic chain using sophis-
ticated compliant contact models seems unrealistic in most
cases. On the other hand, fine calculation of contact forces is
impossible with rigid body approaches since they rely on a
two-time scale of the dynamics and do not care about forces
during the shock but only about the impulses of forces and
their effect on the velocity. Notice, however, that for some of
the reasons listed above (physical parameters identification,
unpredictable motion for multiple shocks, etc), using compli-
ant models may not at all solve this problem.

In summary, some reasons that may motivate one to use
rigid body modeling approaches are:

* Stiff ODEs are to be avoided for real-time applications or all tasks which require fast simulation (because stiff ODEs usually require implicit algorithms). The principal source of stiffness is the contact model. Hence, rigid body contact/impact modeling allows one to avoid stiff differential equations.

* Estimation of contact parameters for compliant models (stiffness, damping, etc) can be quite difficult in practice, and the collision outcome may be highly sensitive to such physical parameter values.

* Simple compliant models of contact may essentially be useful for mathematical aims (study of convergence of penalized problems<fnref n="4" />) and justification of rigid body restitution laws, [18,51,53,74,75,125], but not for numerical simulation due to some fundamental drawbacks [126].

* Rigid models serve as a convenient model for control design purpose (backlash, bipedal locomotion, manipulators during complete robotic tasks, liquid slosh phenomenon, etc) especially when the impact cannot be controlled because the input values should then be chosen too large (consequently the collisions are autonomous).

* Provide good predictions for long term motion simulation, especially if one is not interested in details of local collision (or contact) behavior that involve complex material characteristics, but rather in the effect of collision on a global scale, for example the manner in which vibrations progress from end-effector to base in a flexible robotic manipulator due to impact forces or impulses.

* It can provide the basis for a tool that is easier to use by design engineers.

* When *n* is too large (like in granular matter) penalization models with high stiffness become impossible to simulate due to too long simulation time (to say nothing about first and second items).

To end this part, let us cite Baraff in [127]: *Although the penalty method is useful in some contexts (namely largely static environments) it has become increasingly apparent that the performance of spring-and-damper systems for*

*simulating rigid body motion is inefficient and has unpredict-
able accuracy in dynamic settings...* Notice that compliant
and rigid body approaches may also be mixed, see, eg, [120]
where rigid body and finite element approaches are mixed
for collision detection.

5 NUMERICAL ANALYSIS OF BILATERALLY
CONSTRAINED MECHANICAL SYSTEMS

**5.1 Introduction**

From a general point of view, there are two main classes of
impacting systems: the ones whose orbits undergo collisions
at separated instants (ie, there exists $\delta>0$ such that $t_{k+1} \ge t_k + \delta$ for all $k \ge 0$), and the ones for which this is not true
(see also Section 6).

In the first case, there exist time intervals of strictly posi-
tive measure on which the dynamics are ODEs or DAEs.
One also finds in this category, eg, simple systems with pe-
riodic impacting orbits (vibro-impact systems [6], bipedal
locomotion systems [15], etc). However, even for very
simple nonsmooth systems one cannot always simulate the
motion by combining contact detection with DAE or ODE
solvers. Consider a rigid ball falling to the rigid ground,
where one models the elastic collision with Newton's colli-
sion law $0<e<1$. Then there exist an infinite number of
collisions in a finite time-interval. In a variation to Zeno's
argument, one can then argue that such a situation is impos-
sible, since it requires an object to pass through an infinite
number of collisions in a finite amount of time.

Obviously, the second class is the most general one, and
the one which creates the most interesting challenges in
terms of mathematical analysis and modeling. We shall come
back later to which ingredients one needs to construct a good
simulation software. But it is useful in this setting to recall
first the particular features of DAE simulation.

**5.2 Simulation of DAEs: General methods**

For unconstrained dynamical systems, there are many dis-
cretization methods available. Discretization formulas for
linear systems can be found in [128]. For nonlinear systems,
there are many numerical methods available to solve Ordin-
ary Differential Equations (ODEs), see, eg, [129,130]. This
is not the case for Differential Algebraic Equations (DAEs):
establishing solvers for DAEs is still a very active research
area. Overviews of the theory can be found in [25,29].

Numerical simulation of the conventional formulation of
a constrained mechanical system exhibits severe stability
problems, already, for simple systems. Since simulation of
multibody systems is an activity with a long history, several
computational procedures have been proposed to overcome
the stability problems. These include techniques where a dis-
tinction is made between dependent and independent vari-
ables (a solution is sought through singular-value-
decomposition), equilibrium correction strategies [131],
penalty formulations [132,133], coordinate partitioning
methods [62], predictor/corrector algorithms [134], a differ-
ential algebraic approach [66,67], and projection methods
[135].

<footnote n="4">⁴A fine examination of the planar impact of a particle in an angle [74,75] shows that a penalized problem (damped linear spring-dashpot) yields Moreau's rule Eqs, (11) and (12) for generalized soft impacts in almost all the cases (in a measure theoretical sense).</footnote>
In the literature on numerical integration, DAEs are often characterized by their index [25,29,136]. Roughly speaking, the index equals the number of times the constraints must be differentiated to arrive at a set of ODEs. The index can be viewed as a measure of how far a DAE is from being an ODE. DASSL, designed by Petzold in the early eighties, is capable of solving DAEs which have a low index [137]. Constrained mechanical systems often have an index equal to three [29]. It is well known that index three systems can not be solved directly by standard ODE solvers [135,138]. The dynamics in Eq. (7) have been used in [67] to reduce the index. A multistep (BDF) method with variable step-size, combined with a Newton iteration root-finding algorithm, has been shown to converge. For constrained mechanical systems, stable numerical algorithms are few, and are usually available as research code only. A summary of DAEs simulation algorithms applied to bilaterally constrained mechanical systems can be found in [26] (Ch 7.2).

In engineering practice, the constraint stabilization technique presented by Baumgarte [131] is often applied because it is conceptually simple and easy to implement. Differentiating the position constraint twice gives $\dot{\phi}(q)=0$. It is well known that the numerical solution of this equation can be unstable, and can lead to values of $\phi(q)$ and $\dot{\phi}(q)$ that are far from the desired value zero. The modified acceleration equation: $\ddot{\phi}(q)+2\alpha\dot{\phi}(q)+\beta^2\phi(q)=0$, is (asymptotically) stable for $\alpha>0$. The additional terms in the latter equation can be seen to act as a proportional/derivative control with gains equal to $2\alpha$ and $\beta^2$. Baumgarte also presented the proportional/integral counterpart, for the asymptotic stabilization of holonomic constraints [139].

One problem can readily be seen from the formulation of the stabilization technique: how to choose the coefficients $\alpha$ and $\beta$? Since the stabilization term can be interpreted as a proportional/derivative control law, it is noted that the use of the stabilization term shifts the poles of the system and alters its dynamic behavior. The choice of $\alpha$ and $\beta$ is merely a matter of how fast we want to damp out the constraint violations. Large values of $\alpha$ and $\beta$ lead to high-gain feedback laws. Note that the choice $\alpha=\beta$ yields a critically damped system. It is this choice that is commonly used when Baumgarte's method is applied. In [140], the gains $2\alpha$ and $\beta^2$ are related to the step size that is applied in the numerical algorithm. There it is remarked that their choice of gains tends to damp out constraint violations faster than any other choice, but accumulation of (integration) errors cannot be prevented. Furthermore, decreasing the step size results in larger gains. As a result, the damping terms dominate the numerical solution process: they make the system become numerically stiff. A further analysis of Baumgarte's method can be found in [141] and [142]. In spite of these drawbacks, the constraint stabilization technique of Baumgarte is often applied since it avoids iterative solution of algebraic constraints. This is in contrast to, for instance, a predictor/corrector algorithm and some of the other methods. For instance, the projection method proposed in [135] uses a combination of the numerical solvers known as Backward Difference Formulas (BDFs) and a Gauss-Newton projection method. These algorithms

require iteration processes to obtain values within a certain predefined error level: a number of corrector steps must be applied.

## 5.3 Application to mechanical systems

During simulation studies of unilaterally constrained mechanical systems an expression of the Lagrange multiplier is often used as a model for a force sensor, or simply as a nonlinear expression for the contact force, expressed in terms of the motion on the constraint manifold and the control input [143]. Let us start from Eq. (5) and add a control input $w$ in the right-hand-side. Under the assumption that $\nabla\phi^T(q)$ has full-row rank, and omitting friction terms, one can obtain:

$$ \lambda = -[\nabla\phi^T(q)M^{-1}(q)\nabla\phi(q)]^{-1}(\nabla\phi^T(q)M^{-1}(q) \\ \times(w+Q(q,\dot{q}))+\overbrace{\nabla\phi^T(q)\dot{q}}^{\cdot}). \quad (35) $$

In [66], it is noted that, although Eq. (35) provides an analytical expression for the Lagrange multiplier, it is not always a good starting point to obtain a numerical solution. The presence of numerical errors leads to violation of the constraints and eventually leads to a drift-off from the constraint manifold. As a consequence, one obtains physically meaningless solutions. Usually, for consistent initial conditions, it will take a longer period of time for the drift-off to become noticeable. And if the drift-off remains small, then the corresponding approximate solution may well be acceptable. But, generally speaking, a growing drift-off can not be accepted. Striking examples of drift-off can be found in [66] and [135].

For ease of notation, let us examine drift-off in a first-order formulation (which can be obtained, for example, starting from (6) by setting $x^T=(q^T, \dot{q}^T)^T$, and redefinition of terms).

$$ \begin{aligned} \dot{x} &= f(x) + g(x)w + g(x)C^T(x)\lambda \\ p(x) &= 0. \end{aligned} \quad (36) $$

Instead of finding a discrete-time expression for $\lambda$ directly from (36) another sequence of steps is advocated in [66]. This sequence of steps can be described as *discretize first-substitute next-combine later*. This approach to simulation has been applied to restricted ODEs in [66] and to boundary value problems of Partial Differential Equations (PDEs) in [144]. But the original idea can already be found in [145] where it is applied to index one systems with linear, stationary constraints in combination with the Forward-Euler integration method.

First the equations in (36) are discretized. We will use the Forward-Euler method only to illustrate the concepts, although in general it is not advisable to simulate mechanical systems using the Forward-Euler method. This gives

$$ \begin{cases} x_{i+1} = x_i + h[f(x_i) + g(x_i)w_i + g(x_i)C^T(x_i)\lambda_i^d] \\ p(x_{i+1}) = 0. \end{cases} \quad (37) $$
Note that the constraint is treated in an implicit manner. The notation $\lambda_i^d$ is used to distinguish the discrete Lagrange multiplier from the analytical one in eg, (35). The idea now is to obtain a discrete formula for $\lambda_i^d$ directly from Eq. (37). Let $\epsilon_i$ denote the numerical error at step $i$, and let us assume here that $\epsilon_i$ is a constant $\epsilon$. First assume that the matrices $C$ and $H = \partial p / \partial x$ are constant. The discrete generalized Lagrange multiplier $\lambda_i^d$ is given by

$$ \lambda_i^d = -(CC^T)^{-1} \left( H[f(x_i) + g(x_i)w_i] + \frac{p(x_i)}{h} \right). \quad (38) $$

We will refer to $p(x_i)/h$ as a compensation term in the remainder. One can now obtain [66]:

$$ p(x_{i+1}) = hCC^T\epsilon. \quad (39) $$

No error accumulation can take place. Indeed, if $h \to 0$ one has $p(x_{i+1}) \to 0$, as desired.

In contrast, if one used the analytical expression of the Lagrange multiplier as in Eq. (35) it is shown in [66] (see also [65], Ch 9) that one obtains, with neglect of the higher-order terms:

$$ p(x_{i+1}) = p(x_0) + h \sum_{j=0}^{i} C(x_j) C^T(x_j) \epsilon_j. \quad (40) $$

From the latter equation, it follows that once an error is made the solution is not on the constraint manifold. Now again consider the special case where $\epsilon_i = \epsilon$, $\forall i \in \mathbb{N}$, with $\epsilon$ constant, and that $C$ is a constant matrix. This gives

$$ p(x_{i+1}) = p(x_0) + t_i CC^T \epsilon. \quad (41) $$

Note that it makes no sense to let $h \to 0$ since $h$ does not even appear in Eq. (41). And if $t_i \to \infty$, for instance because we are interested in an equilibrium solution, one even has that $p(x_{i+1}) \to \infty$! Even if the initial conditions are consistent with the equality constraints, error amplification is inevitable due to the presence of numerical errors. Each error source will contribute to the drift-off. This is one of the reasons why many simulation codes project the state on a time $t_{i+1}$ back to the constraint surface for example using Newton-Raphson, before proceeding the time simulation.

The expression for the discrete generalized Lagrange multiplier $\lambda_i^d$ is useful also in combination with other explicit and implicit integration routines, and for nonlinear constraints as well. For this one treats the compensation term as a constant on an interval $[t_i, t_{i+1})$ and evaluates all other functions in the points needed by the numerical method that is applied. For application to mechanical systems, the expression for the Lagrange multiplier on intervals $[t_i, t_{i+1})$ now reads:

$$ \lambda_i^d = -[\nabla \phi^T(q) M^{-1}(q) \nabla \phi(q)]^{-1} \left( \nabla \phi^T(q) M^{-1}(q) \right. \\ \left. \times (w + Q(q, \dot{q})) + \frac{\dot{\nabla} \phi^T(q) \dot{q} + \nabla \phi(q(t_i)) \dot{q}(t_i)}{h} \right). \quad (42) $$

Using Eq. (7) one also obtains:

$$ \mu_i^d = -[\nabla \phi^T(q) \nabla \phi(q)]^{-1} \left( \nabla \phi^T(q) \dot{q} + \frac{\phi(q(t_i))}{h} \right). \quad (43) $$

Observe that analytically the constraints are strictly satisfied, and one has $\lambda_i^d = \lambda$ and $\mu_i^d = \mu$.

Using the “discretize first—substitute next—combine later” sequence of steps, a numerical method can be obtained that has the property that it is robust with respect to errors in the initial conditions, and stable with respect to errors made during numerical integration. Moreover, the use of the discrete generalized Lagrange multipliers $\lambda_i^d$ and $\mu_i^d$ does not yield numerically stiff equations when the time-step is reduced: the term $h$ in the denominator is canceled again. Since the compensation term is fixed on each time-interval, it is easy to combine the method with variable stepsize methods (see [146,147] for a variable stepsize method). The numerical method can be extended to cover constraints $\phi(q,t)=0$, ie, constraints that depend explicitly on time. These constraints arise, for instance, in dynamic path planning of robotic manipulators, and also in flight-path management of airplanes. Most of the theory remains valid, but in the case of, for example, Runge-Kutta-4 there are additional requirements on higher derivatives of the constraints that need to be satisfied. This still leads to stable numerical simulation of bilaterally constrained dynamical systems.

# 6 NUMERICAL ANALYSIS OF UNILATERALLY CONSTRAINED MECHANICAL SYSTEMS

## 6.1 General motivations for new specific schemes

The main problem in simulating DAEs is the stabilization of the constraints. For MDIs, additional difficulties occur:

* Complementarity conditions are an essential part of the dynamics and have to be carefully incorporated into the numerical scheme. In particular, Baumgarte's stabilization technique is not suitable since the resulting multipliers' signs are physically meaningless, and cannot be used (even in the codimension one constraint case) to detect release.

* Detection of contact instants and re-initialization of the state. In the case where one is able to detect analytically what the precise impact point and time are, there are still problems in the numerical approximation of these instances. First, one must be sure that an impact is not missed (see Section 6.5). Next, since machine zero is not identical to zero, there will almost always be a small mismatch between, for instance, actual impact location and simulated impact location. It must then be decided whether a collision will take place or that the bodies remain in contact with each other. Clearly, a priori determination of all possible contact points, necessarily including the velocity information as done in [148] and [149], is an asset. Last but not least, one must decide whether or not, after constrained motion a release takes place. Since zero force level can not be maintained numerically and since DAE solutions do not yield exact satisfaction of a constraint, also in this case a priori determination of all possible release points, including the velocity and force information
is a benefit [149]. In Section 6.4 we will come back to these difficulties in relation to event-driven software schemes.

* Accumulation point of the sequence {$t_k$}$_{k \ge 0}$ occur frequently. One has to impose a threshold or derive local analytical forms of the impact Poincaré map [150]—which are similar to the bouncing ball dynamics. For $n=m=1$, the transition to the permanently constrained mode can be decided via a threshold. For $n \ge 2$ and/or $m \ge 2$, things complicate because generally the coordinates are dynamically coupled. If one constraint is decided to be saturated because numerically one is unable to detect rebounds any longer (ie, computer accuracy may be attained), then the overall motion may be drastically modified. In general, the methods based on impact detection will provide good results for:

—Periodic motions with finite number of impacts per period, for all $n, m$.

—Periodic motions with finite or infinite number of impacts per period, for $n=m=1$ [151–153].

The multiplicity<fnref n="5" /> of the shock is not easy to calculate, see works on Newton’s cradle [154] and on impact of a particle with an angle [73,155].

## 6.2 MDIs are not concatenation of ODEs, DAEs, MDEs

It is important to realize that an MDI as in Eq. (2) is quite different from ODEs, DAEs, DIs $\dot{a}$ la Filipov, MDEs<fnref n="6" /> of the form

$$ \dot{x} = f(x,t) + g(x,t)\dot{w}, \quad x(0^-) = x_0 \quad (44) $$

where $x \in \mathbb{R}^n$, and the control input $w \in \mathbb{R}^m$ is of bounded variation, and even MDEs as in [156], which are defined as:

$$ \begin{cases} \dot{x} = f(x) & \text{if } t \neq t_k \\ x(t_k^+) = x(t_k^-) + a(x(t_k^-)) & \text{if } t = t_k \\ x(0^-) = x_0 \end{cases} \quad (45) $$

for a sequence {$t_k$}$_{k \ge 0}$, where $t_k \to +\infty$ as $k \to +\infty$, $t_k$ may or may not depend on $x^7$. The following fundamental discrepancies between the last two differential equations and our MDIs are [18,87]:

* It is clear that the overall dynamics with complementary-slackness mechanical systems is far from being only a differential equation with impulsive inputs. The way the set of indices for active constraints and sticking/slipping contacts are refreshed is a fundamental part of the (hybrid) system [157]. Only the very simplest dynamics of systems as in Eq. (2) seem to resemble Eq. (45). But they are of different nature, as explained below.

* The commutativity property of the vector fields $g_i(x) \in \mathbb{R}^n$, $i \in 1, \dots, m$ allows one to decouple the MDE $\dot{x} = f(x) + g(x)\dot{w}$ into simpler systems for which existence and uniqueness of solutions are easier to prove. These commutativity properties are always satisfied for mechanical systems, but they do not at all solve the problem of multiple impacts [18].

* If $f(x)$ and $a(x)$ in Eq. (45) are linear then the whole nonsmooth system is linear as well. On the contrary, linearity of the continuous vector field in a mechanical system subject to a unilateral constraint does not at all imply linearity of the total system (see [18], Ch 1).

* The system in Eq. (45) and the developments in [156] do not allow for finite accumulation points of discontinuities.

* The product $g(x)\dot{w}$ is quite meaningful for unilaterally constrained mechanical systems because $g(x)$ contains only positions and only the velocities are discontinuous [18,87]. Otherwise specific tools have to be developed to render Eq. (44) meaningful. The formulation in Eq. (45) is a way to overcome such problems.

* The abundant literature on (Lyapunov) stability properties of impulsive differential equations as in Eq. (44) deals only with the stability of the fixed point of $\dot{x} = f(x)$, where the remaining terms are considered as disturbances. Stability studies of unilaterally constrained mechanical systems require other tools. Mainly because [18–20] i) they are hybrid dynamical systems with different modes (in other words there is a natural Discrete Event part whose states correspond to the modes associated to the index sets of sticking and sliding contacts), ii) the fixed point of the overall system may not at all correspond to the one of the continuous dynamics (think of the bouncing ball example in which the vector field $\ddot{q} = -g$ does not even possess any fixed point!), and iii) the stabilization goals may be quite different depending on the task (stabilizing jugglers or hopping robots [14] is totally different from stabilizing a manipulator along a time-varying trajectory including free-motion, constrained motion and impacting phases [19,20]).

* Nonsmooth mechanical systems as in Eqs. (2) and (3) have solutions which are generally discontinuous with respect to initial conditions [51,39]. This is not the case for MDEs as in Eq. (44) as can be easily seen taking, for instance, $g(x,t) = g(t)$ [158].

* A lot of work has been dedicated to simulation of bilaterally constrained mechanical systems, since they are index 3 DAEs (see [25] for further references and see Section 5). But as we said above for hyperstatic systems, LCPs have to be solved at each step and integrating switching DAEs is far from being sufficient. Baumgarte’s method is not suitable since it does not guarantee that the sign of the multiplier has any physical meaning.

In summary, the MDIs in Eqs. (4), (22), (26)–(28), and (32)–(33) are quite different in nature from the other types of MDEs. They are not a simple concatenation (except in very simple cases) of ODEs and/or DAEs separated by reinitializations of the state (the velocity). In order to reinforce

<footnote n="5">⁵The multiplicity is to be understood here as the number of shocks that the system undergoes when it strikes in a neighborhood of the singularity of $\partial K$. It should perhaps better be called the order $O$ of the collision, which can be simple or multiple. In this sense, the ball falling on a ground undergoes, in finite time, a simple shock of order infinity. A 2D particle striking an angle undergoes a 2-shock of finite or infinite order depending on the angle, initial data, restitution coefficients. Actually, these two notions require more accurate definitions: if $O < +\infty$ then one just faces a series of simple shocks. So the only relevant case is $O = +\infty$.</footnote>

<footnote n="6">⁶Measure Differential Equations.</footnote>

<footnote n="7">⁷Evidently $k$ here is an index and has no relationship with the stiffness in Section 4.2.</footnote>
these ideas, let us consider the following third order system ([39], Ch 7.3; [159]), which belongs to the class of Linear Complementarity Systems [42,160]:

$$ egin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = x_3 \\ \dot{x}_3 = \lambda \\ \lambda \ge 0, \quad x_1 \ge 0, \quad \lambda x_1 = 0 \\ x^T(0^-) = (0, -1, 0). \end{cases} \tag{46} $$

The state has to jump initially to the value $x^T(0^+) = (0,0,0)$, and with $\lambda(0)=0$. Then the system stays at rest keeping the same state for all $t \ge 0$ (note that the DAE corresponding to Eq. (46) reduces to the trivial system $x \equiv 0$, $\lambda \equiv 0$). Let us apply the following backward Euler (implicit) scheme to the system in (46):

$$ \begin{cases} \frac{x_{i+1} - x_i}{h} = A x_{i+1} + B \lambda_{i+1} \\ x_{1,i+1} \ge 0, \quad \lambda_{i+1} \ge 0, \quad x_{1,i+1} \lambda_{i+1} = 0. \end{cases} \tag{47} $$

The matrices A and B are easily identifiable from (46) and $x_i$ denotes the discretized value of $x$ at step $i$. Roughly, one constructs a LCP whose unknown is $\lambda_{i+1}$ by inserting $x_{1,i+1}$ into the complementarity conditions (which requires that $(I - hA)$ be invertible, which is true for $h$ small enough). This way of doing is classical in time-discretization procedures, see Sections 7.1 and 7.4. Under certain conditions (satisfied for this example) the constructed LCP($\lambda_i$) in Eq. (47) possesses a unique solution for sufficiently small $h$. The initial state jump is computed by solving the LCP at the first step, and one finds $(x_{1,0}, \lambda_0) = (0,1/h^2)$. Then the solution should converge to the null solution. The application of the algorithm in Eq. (47) yields [39]: $(x_{1,i}, \lambda_i) = (i(i+1)/2h, 0)$ for $i \ge 1$. Moreover it can be shown that the nonimpulsive part of $x_1$ satisfies $\|x_1\| \ge N_h(N_h - 1)/2h^{3/2}$, where $N_h$ is the integer part of $T/h$ for some $T>0$ and $\|\cdot\|$ is the Euclidean norm. Therefore $x_1$ diverges as $h \to 0$!

Following [72], let us consider linear DAEs. If the initial values are inconsistent (ie, they do not satisfy the constraints), then fixed-step k-step backward differentiation formulas (BDF) produce solutions which may diverge during the first steps as $h \to 0$, for $k \le 6$. However, they are $O(h^k)$ after a maximum of $(\bar{m}-1)k+1$ steps, where $\bar{m}$ is the DAE index. Therefore, the example in Eqs. (46) and (47) shows that one should be very cautious when applying an implicit Euler discretization to a nonsmooth system as in (46): the results of convergence which are available for ODEs or DAEs may absolutely fail when considering complementarity systems. Indeed it is really the complementarity conditions (and consequently the LCP($\lambda_i$)) which create the divergence phenomenon. The initial jump drives the state away from the constraint, but convergence is not recovered afterwards. Fortunately, this sort of problem is shown not to occur for linear mechanical systems in [39,159], and more generally for dissipative systems (see [117] for details on dissipative systems theory). The example in Eqs. (46) and (47) demonstrates that the existence of solutions to the (discrete-time)

LCP($\lambda$) does not imply at all the stability of the overall integration scheme nor convergence of its solution as $h \to 0$. This fact is quite similar to what occurs for the continuous-time dynamics, [18] Remark 5.18. Another example that demonstrates the importance of the (discretized) complementarity conditions can be found in [91]. It concerns a bouncing ball, whose dynamics is in Hamiltonian form with state $(q,p)$, discretized by an implicit Euler midpoint rule. If one uses $0 \le \lambda_i \perp (q_{i+1} + q_i)/2 \ge 0$, then the effective restitution is not elastic and the energetical behavior is desastrous. If one uses $0 \le \lambda_i \perp (p_{i+1} + p_i) \ge 0$, the energetical behavior is perfectly elastic (this had also been noticed in [28,96]).

It has been stated [161] that “Differential/Algebraic equations are not ODEs.” We conclude that: *Measure Differential Inclusions that represent complementary-slackness mechanical systems are not concatenations of ODEs and DAEs, neither MDEs.*

## 6.3 Simulation algorithms

It is difficult to classify the existing methods of numerical integration since there are many different criteria that can be used to fulfill such a task: acceleration-force/velocity-impulse; computation of the contact forces or not; resolution of LCPs/NCPs for all times, at certain times, never; frictionless constraints; Coulomb friction, regularized Coulomb; Signorini conditions; polyhedral approximation; revolution cone of friction; 2D/3D order of the discretization, contact point; finite-element methods, quasistatic/dynamic, inclusion of complementarity conditions in the discretization; calculation of impact instants or not; proof of convergence towards a solution of the rigid body dynamics or not; inclusion of multiple shocks; implicit or explicit form of the numerical scheme etc. Following [101], we choose the following classification in this paper (another classification is chosen in [115]):

* Event-driven schemes: one integrates the DAEs between state re-initializations. One needs some basic ingredients that we will describe below.

* Time-stepping schemes: one considers a time-discretization of the nonsmooth dynamics (MDI) including complementarity conditions and impact rules. The whole set of discretized equations/conditions is used at each step to calculate the future state.

* Penalized-constraint schemes: the unilateral constraints are replaced by some compliant model of the contact.

Consequently, we will first describe the main ingredients that any good event-driven software should incorporate, and then we shall describe the most advanced methods that have been proposed in the recent years by several authors and which incorporate (implicitly or explicitly) such modules. Evidently we do not deal with penalized schemes and stiff ODEs here, as announced before. The problem of discretization of the continuous dynamics (ODEs or index 3 DAEs) has been discussed in Section 5, as well as the problem of constraint stabilization.
## 6.4 Overview of modules for event-driven schemes

It is clear that there is a need for simulation tools that provide routines that can be used to simulate MDIs. In the remainder of this section we will discuss some of the ingredients that should be incorporated in event-driven simulation softwares. It must be stated that even for bilaterally constrained mechanical systems, commercial tools are not abundant and in many cases the numerical code that implements the transition rules and discontinuities due to uncontrolled collisions must be supplied by the user. In this section an overview of different modules is given.

In order to simulate unilaterally-constrained dynamical systems, one needs to know whether or not there will be a jump in the state or its derivative (the acceleration), and if so, how this jump is made. A simulation of contact must use models for the collision maps, the transition rules and the dynamics. For affine nonlinear systems, [65,149] give algorithms that compute all possible contact and release sets off-line. Using these sets, explicit expressions for the collision maps can be made off-line as well. During simulation, once the contact point is known, a simple check followed by one function evaluation then suffices to obtain the desired results. Due to discretization however, there remain a number of problems with respect to simulating uncontrolled and controlled contact. The problem with (un)controlled contact is the approximation of the time that contact with the boundary set is made. The problem with controlled contact is the validation of the transition rules for contact and release. Since real-time simulation often involves trade-off between accuracy and computation time, in [65] (Ch 9), parameters are introduced to aid the user of a simulation program to customize simulation to his/her specific needs.

The introduction of the Lagrange multiplier makes detection of contact and release during simulation studies nontrivial, as pointed out in Section 3.1. For instance the decision that a release takes place cannot be based on a check whether or not $\phi(q)>0$, since an active Lagrange multiplier will ensure that $\phi(q)=0$. It must be decided when to (de)activate the Lagrange multiplier $\lambda$. The transition from free motion to constrained motion, and vice versa, or plastic collision detection can be based on the following analytical rules [43,149] (for codimension one constraints).

**Contact rule:** if $\phi(q)=0$ and

$$\nabla \phi^T(q) \dot{q} = 0 \quad \text{then activate } \lambda \quad (48)$$

**Release rule:** if $\lambda<0$ then deactivate $\lambda$.

First we discuss the problem of determining the time that a trajectory makes contact with the boundary set. Assume that at time $t_i$ there holds $\phi(q(t_i))>0$. If contact is made in the interval $(t_i, t_{i+1})$ then this can be detected only if $\phi(q(t_i)) \cdot \phi(q(t_{i+1})) \le 0$. This also implies that if a controlled contact and a controlled release take place in the interval $(t_i, t_{i+1})$, (and the motion stays on a boundary for a small period of time), this can be detected numerically only if the time-step is small enough. We assume that this is the case.

Depending on the characteristics of the numerical solution, it may well be that $\phi(q(t_{i+1}))$ is significantly smaller than zero. This implies that even with a fixed step-size, so-called step back is necessary. Step back means that the time-step $h$ is adjusted and the simulation is started again from time $t_i$. This means that the same dynamic equations are solved again. And if the new estimate of $q(t_{i+1})$ is not satisfactory, again step back is necessary. Clearly such an iterative procedure may increase computation time significantly.

The method of false position or regula falsi can be used to establish the time-instant $t^*$ where a trajectory makes contact with the boundary set. Since we assume stable numerical integration, the numerical approximated trajectory and the analytical trajectory are related through the accuracy of the numerical ODE solver that is used, and to the approximation error on the initial data (that is needed after each impact). This means that it makes sense to search for a time instant $\hat{t}$ that is close, but not necessarily equal to $t^*$. The step-size is adapted to a value $\hat{h} \in (0,h)$. This gives $\phi(q(t_i+\hat{h})) \approx 0$. This approximation of the contact point may already be acceptable to the user of the simulation program. And using this value of $\hat{h}$, a linear (or nonlinear) interpolation with the numerical value of the velocities at times $t_i$ and $t_{i+1}$ may give an acceptable value of the velocity component at contact. But we can go on by restarting the simulation from time $t_i$ by using $\hat{h}$ as the new time-step. Now the obtained numerical approximation is checked again to see if $\phi(q(t_i+\hat{h})) \le 0$. If this holds then the procedure outlined above can be repeated until a point is reached for which the position constraint is not satisfied. The approximation of the time of contact, $\hat{t}$, is now set equal to $t_i+\hat{h}$, where $\hat{h}$ is the last update such that $\phi(q(t_i+\hat{h})) \approx 0$. Next, the time $t_i+\hat{h}$ is also used to approximate the velocity vector at the time of contact. Subsequently we set the time-step to its old value $h$ again as the time-step $\hat{h}$ may be too small to maintain real-time simulation. A parameter $\gamma$ should be introduced to set an a priori bound on the number of iterations. The trade off between accuracy and computation time is then put in the hands of the user. In [162,163], $h$ is divided by two until $\phi(q(t_k))$ is smaller than a given tolerance. Another adaptive step-size method is proposed in [164]. Various methods similar to what is described above to compute $t^*$ are discussed in [151–153] as well as their influence on the scheme order when combined with Newmark, RK24 and Dormand-Prince RK (the benchmark is a one-degree-of-freedom system with a single constraint).

Even though the procedure outlined above is simple, it still uses an iteration process. If the time needed to execute the procedure violates real-time simulation, the initial linear interpolation step may be replaced by a higher-order interpolation method using information at velocity level. Another promising approach is to use in the iteration process itself only a part of the dynamics equation: $M(q)\ddot{q} = w$. This choice is motivated by the fact that the contact set is independent of the system matrices $M$ and $Q$ and the importance of the matrix $M$ in calculating the impulse [149]. This is still an active research area for real-time simulations of unilaterally constrained robotic manipulators (see also Section 6.5).
Next it must be determined whether one is dealing with a collision or not. In case plastic collisions are considered, after detection of contact, the Lagrange multipliers are activated and the motion remains on the surface, see Eq. (48). When elastic collisions are possible, there is another numerical problem since $d/dt\{\phi(q)\}(t*)$ will in general not be zero numerically. This problem is similar to the problem in bilaterally constrained mechanical systems when the reaction force is approximately zero. A heuristic approach is the following. If at a contact point one has that $v \le d\phi/dt(q(t*)) \le 0$, for a user defined value of the parameter $v$, then the motion is assumed to proceed on the boundary. If $d/dt\{\phi(q)\}(t*) < v \le 0$ then the motion is to proceed by the use of an uncontrolled collision map. To detect release from the surface a check is made whether or not $\lambda_i^d < 0$, where $\lambda_i^d$ is in Eq. (38). If this inequality holds true, then the Lagrange multiplier is deactivated. Since the compensation terms are usually much smaller then the value of $\lambda_i^d$ itself, these compensation terms have no influence on this decision [65].

There remains a problem when the desired (simulated) contact force during motion on the boundary is very small. But then controller design is also difficult: the slightest deviation from the desired path will mean that release takes place, although control is aimed at maintaining contact. This is the reason that in practice a certain amount of (constant) normal force on the surface is chosen.

## 6.5 Collision detection between bodies

This section and the next one concern a very important module: the management of contacts status.

Approximation of the shapes and approximate calculation of impact times are generally CPU-time intensive tasks [120]. Many works have been dedicated to collision detection, eg, [148,165–167] to cite a few. Roughly this module requires to calculate, explicitly or implicitly, the expressions for $\phi(q)$ and solve $\phi(q)=0$ (ie, determine the points that are going to touch—which are not necessarily the ones which are the closest at the instant of the computation, so several pairs of points have to be watched simultaneously). Even in very simple cases such as one degree-of-freedom systems, various numerical methods may be used to calculate the times $t_k$ such that $\phi(q(t_k))=0$ [151–153]. Their influence on the algorithm properties (consistency, order) may be significant.

The main problem is that an exact analytical description of the objects shapes, even when this is possible, is quite time consuming. Secondly one has to calculate with a suitable numerical routine the times $t_k$. In case of accumulation of impacts and for multiple contacts, the problem is harder because the influence of deciding the end of the series $\{t_k\}_{k \ge 0}$ according to the machine accuracy, is not easy to quantify. Micro-collisions phenomenon [107] prove that it is possible in some cases that there is a large quantity of re-bounds, but finite number of collisions, and an escape out of $\partial K$ after a finite time. Things even complicate for multiple impacts. What is the influence on the long run motion if one decides instead that one constraint becomes active?

As pointed out in Section 2.2, it may not be possible to define all the constraints $\phi_j(q) \ge 0$: there would be too many! Hence one usually employs procedures that eliminate useless constraints, ie, those bodies which are too far one from each other to be likely to collide in the next future steps of integration. Consequently one implements rough tests that select the bodies which may collide, and fine tests to compute the collision times [120]. Rough tests usually consist of surrounding the bodies by simple volumes (spheres, boxes) and watching whether they overlap or not. Concerning the finest tests, the main approaches are (see [166] for a review):

* Classification of typical contacts and geometries [163,168–172]. In other words, process the real surface of the objects and the type of contact (circle-circle, circle-line, angle-line etc). These methods are essentially studied in the Mechanical Engineering literature. They are restricted to certain types of geometries contained in the available library developed for the software. If the bodies surfaces are simple enough to be described by analytical curves, one gets an explicit function $\phi(q(t))$ (see for instance the developments in [83] for the derivation of $\phi(q)$ from a local frame at the contact point). The next step is to solve numerically $\phi(q(t))=0$-which can be done with a Newton-Raphson method or a polynomial root finding routine, since in case of several roots Netwon-Raphson may compute the wrong zero and there is penetration before the algorithm decides that contact has occurred. Others [162,163] use a time step halving process until $\phi(q(t_k)) = 0$ is satisfied within a specified tolerance. We have already outlined this problem in Section 6.4. For instance, for two bodies with parametric surfaces $\phi_i(u_i, v_i, t)=0$, $i=1,2$, one has $\phi(q,t) = \phi_1(u_1, v_1, t) - \phi_2(u_2, v_2, t)$ and one faces a nonlinear 5-dimensional root-finding problem [148]. These methods are however less fast and more complex to implement than the 2-dimensional ones [127].

* For 2D systems, one can approximate the bodies $B_i$, $i \in \{1, \dots, N\}$ by polygons made of edges and nodes $N_i$. Two main methods are used [120]: the node-in-polygon test (NIPT) and the ray-crossing approach (RCA), see Fig. 2a,b, respectively. Let $n_i$ be the number of nodes $N_i$ in polygon $i$, and $\alpha_{ii}$ the angle $(N_j N_i, N_j N_{i+1})$. Then if $\sum_{\text{nodes of } B_i} \alpha_{ii} = 0$, the node $N_j \notin B_i$, if $\sum_{\text{nodes of } B_i} \alpha_{ii} = 2\pi$, then $N_j \in B_i$: the bodies intersect. The RCA consists of looking at the number $\tilde{n}$ of intersections of a straight halpline (a ray) emanating from $N_j$, with the polygon $\partial B_i$. If $\tilde{n}$ is odd then $N_j \in B_i$, if $\tilde{n}$ is even then $N_j \notin B_i$. The RCA is more robust than the NIPT. Both methods are $O(n_i n_j)$ for two bodies $B_i$ and $B_j$. However, their generalization to 3D systems is not easy [120].

* Approximation of the objects surfaces and of the impact times by bounding boxes methods [166] are more efficient for 3D systems. These methods are essentially studied in the Computer Science literature. If the bodies are convex and subject to gravity (or more generally to any vector field that is integrable) it is possible to approximate the distance $\phi(q)$ and to calculate a lower-bound on the impact time [165]. The approximation can be refined as much as the constraints (desired accuracy, speed of computation)
<figure><img src="image_3.png" /><figcaption>Fig. 2 Collision detection methods</figcaption></figure>

permit to do it. In [148], it is pointed out that just watching
positions to determine collision times cannot work since
contact may occur between two sampling instants $t_i$ and
$t_{i+1}$ while $\phi(q(t_i)) > 0$ and $\phi(q(t_{i+1})) > 0$. So including
the velocity information in the algorithm is mandatory.
Adaptive subdivision of the bodies into simple volumes
(polygons or polyhedra [120], spheres [165,166], rectangu-
lar prisms [148]) and incorporation of a Lipschtiz bound-
edness condition on $\phi(\cdot)$ allows one to approximate the
collision times [148,173]. This method is called bounding
box schemes: each object is surrounded by bounding
boxes. When these boxes overlap, the objects must be
close one to each other. Then a more accurate collision test
is made once more. Bounding box schemes allow one to
avoid testing all possible contacts ($=O(N^2)$ for $N$ bodies),
but to focus on objects in close proximity only. Roughly
speaking, the Lipschtiz bounds permit to approximate the
next step motion of each simple volume (or surface) and to
determine if a collision has occurred. A refinement of the
mesh can be used to increase the accuracy of the collision
time computation, in an adaptive way. These methods ap-
ply well to convex bodies. Nonconvex bodies can be de-
composed into convex parts to be treated. Voronoi regions
for polytopes [167,174] are used to maintain a list of clos-

est distances during the simulation<fnref n="8" />. The change in
Voronoi cells from one step to the next one is usually
small, facilitating the calculations. An implementation of
the Lin-Canny algorithm with a running time linear in $N$,
can be found in [175]. Baraff [176] proposes a coherence
based bounding box test that is $O(N)$.

## 6.6 LCP and NCP solvers

**LCPs:** As we saw in Sections 3.1, 3.2, and 3.3, the dynamics in Eq. (2) with or without friction, can be transformed in a form involving LCPs or NCPs which is useful to cope with detachment, sticking-sliding, impacts. Notice from Eq. (8) that if $B \ge 0$ (componentwise), then $\lambda = 0$ is a solution to the LCP. Problems arise when $B_i < 0$ for some $i$. The most popular algorithm used to solve LCPs, Lemke's algorithm, is motivated by this observation [110,177,178]. Basically this is a pivoting method that converges in a finite number of pivots provided the LCP matrix $A$ in (8) is a copositive matrix. In the worst case, Lemke's algorithm is exponential but is expected to be polynomial in $m$. Lemke's algorithm is initial-

<footnote n="8">⁸A Voronoi cell associated to an object consists of the set of points whose distance to this object is the smallest. The object can be a node, an edge, a face. See Fig. 2c.</footnote>
ized using a so-called covering vector [89], or supplementary variable [178], or initial ray [177], that is chosen arbitrarily or even randomly [89]. Lemke's algorithm is formulated in [177] as finding the zero of a piecewise linear function, using a homotopy method. Even in the nice case where A is copositive, Lemke's method is known to fail when the dimension is too large (because of accumulated errors of pivoting), or even for medium-sized problems if the diagonal of A has null entries [115]. Further interesting properties of Lemke's algorithm have been proved in [179]. Other methods to solve LCPs exist, such as Mangasarian which consists of using a nonlinear (complementarity) function whose zero is the LCP solution [180], then a root-finding algorithm. Another homotopy method called Katzenelson's algorithm exist [181]. Mangasarian's method has been used in [182]. As we pointed out above, LCPs and Quadratic Problems under constraints are equivalent when A is copositive (then the QP is convex and its solution satisfies the so-called Karush-Kuhn-Tucker conditions which are complementarity conditions between N and a Lagrange multiplier—a slack variable). In [183], Interior Point algorithms are described that solve LCPs as QPs under constraints. In [89], a comparison between Lemke's and a new specific Interior Point algorithm is proposed. 35 different problems are tested, randomly generating 328 data sets, with LCPs of dimensions ranging from 2 to 170. The IP seems to supersede Lemke in terms of the number of problems it can solve when Lemke's covering vector is chosen arbitrarily. See [115], [178], and [184] for a comparison between Lemke, Gauss Seidel and augmented Lagrangian methods, using distinct-element models of various mechanical systems: Lemke's algorithm is shown to provide the smallest CPU computation time.

Baraff developed a Dantzig algorithm [185]. Let there be $m$ frictionless contacts with a set of complementarity conditions between $\dot{V}_N^{(j)}$ and $N^{(j)}$ at each contact. Initially the contacts are classified into 2 sets: $\dot{V}_N^{(j)}=0$ and $\dot{V}_N^{(j)}>0$, whereas one sets $N^{(j)}(0)=0, \forall j \in \{1, \dots, m\}$. Dantzig algorithm is a pivoting method of the indices between the 2 sets. One starts by computing suitable values of $\dot{V}_N^{(1)}$ and $N^{(1)}$, then one passes to contact 2 and adjusts both contacts 1 and 2. Indeed modifying the data $N^{(2)}$ and $\dot{V}_N^{(2)}$ must in general modify $N^{(1)}$ and $\dot{V}_N^{(1)}$, see Eq. (8). Then contact 3 is adjusted, taking into account 1 and 2 as well. And so on until contact $m$ is attained. Then all contacts satisfy the complementarity conditions and the LCP($\lambda$) in (8). The algorithm terminates if $A$ is copositive. Extensions to frictional contacts are possible, in 2 or 3 dimensions. It is noted in [185] that in many practical cases $\nabla\phi$ is constant, consequently $B$ in (8) is in the column space of $A$: this property may be used for singular $A$ (eg, hyperstatic systems). Moreau uses a Gauss-Seidel method [78], see Section 7.2.2.

In general, users may be tempted to directly connect their software with available codes. However the use of packages has the following drawbacks [185]:

* Interfacing the numerical software package with the simulation software may not be easy.

* The available codes are often implemented as research codes, especially for QPs.

* The tuning of the package adjustable parameters (numerical tolerances, iteration limits etc) may not be easy.

* The packages are often *black boxes* whose codes are not accessible to the user. The origin and remedy of numerical problems may therefore be difficult to isolate.

Algorithms for 2D frictional contact problems in the quasistatic case have been presented and compared in [186]. It is shown that this reduces to finding the zero of a function in a closed domain. Newton's method seems to supersede interior point ones.

**NCPs:** The NCP in Eq. (16) is equivalent to solving $\min[y,g(y)]=0$. Such nonsmooth equations can be treated with nonsmooth Newton methods [187]. As a generalization of the variational expression in Eq. (21), (16) is also equivalent to the variational inequality [111]: find $y \in \mathbb{R}^+$ such that $(z-y)^T g(y) \ge 0$ for all $z \in \mathbb{R}^+$. From a theoretical point of view NCPs are still an active research area, (see eg, [188] and [189] and references therein).

## 6.7 Event-driven algorithms

Basically these are algorithms which integrate the motion between events (shocks or stick-slip transitions) and use some event detection procedure. They are of the force-acceleration type and require the computation of contact forces. The works therefore focus on the improvement of the modules in Sections 3.2, 3.3, 6.5, and 6.6. The time-discretization problems are not treated explicitly, but the authors generally merely use available DAEs or stiff ODEs solvers. Roughly, the main discrepancy between what we call event-driven and time-stepping schemes is that the former are closer to the continuous-time description of the dynamics in Eq. (2) with a two-time-scale dynamics (which is close in spirit to the hybrid dynamical systems point of view [42]). The latter are a real difference equation approach, and therefore, better lend themselves to convergence analysis, a crucial property. On the contrary event-driven algorithms do not seem amenable for convergence analysis. This may have important consequences for the simulation of complex dynamics. The first three groups of algorithms in Sections 6.7.1, 6.7.2, and 6.7.3 clearly supersede the others in the way they formulate the nonsmooth dynamics. They are briefly presented here since most of the specific material they contain has been presented elsewhere in this paper.

### 6.7.1 Pfeiffer and Glocker's formulations

These authors proposed in [7] a general formulation of the dynamics in Eq. (2), at acceleration-force level. The basic convex analysis tools described in Section 3.4 are used to derive complementarity conditions and write friction and impact rules with LCPs or NCPs [83]. Two-dimensional friction (planar systems) are treated in [7], Lemke's algorithm is used. Extensions to 3D systems is studied in [83], using polytope approximations of the friction cone and convex analysis tools which allow one to express such approximate Coulomb's law with complementary slackness conditions. It encompasses collisions with friction and tangential restitution, see Section 3.3. The work in [92] uses the second formulation in Eq. (21) to express the complementarity conditions in Eq. (10) and dry friction in dimension 3 (see (23) and its acceleration counterpart). The dynamics in Eq. (2) is
transformed in an implicit equation for the acceleration and N, $F_t$, to be solved at each step by a root finding algorithm (a subroutine NEWT from [190] is used, and the dynamics are integrated with a RK-Fehlberg with automatic control of the integration step). Again convex analysis tools are at the core of the proposed method. The numerical results are compared to experiments (a planar mass on an oscillating oblique plane). In another work [182] the same authors use an augmented Lagrangian approach to solve the same problem, with a linearized friction cone. The numerical predictions fit well with experiments.

In most of these schemes, the integration of the continuous dynamics is done with a RK scheme. The existence of solutions of the obtained LCPs (Painlevé paradoxes) is discussed but not treated. Some of these results have been implemented in a commercial software [191,192].

### 6.7.2 Modified Moreau's scheme

The algorithms presented in [4] and [5] may be seen as an adaptation of Moreau's schemes (in particular the Gauss-Seidel method to solve the multi-contact problem with friction). They provide the user with the event occurences and the contact forces values, which are of primary importance for engineers doing virtual prototyping. Due to these industrial constraints, an event-driven force-acceleration formulation is more suitable than a time stepping velocity-impulse one. Also low-order algorithms (Euler) may not be very accurate when applied to systems with sparse events, a drawback that is well-known for free-motion systems [26]. The algorithms have been thoroughly tested on the company Schneider Electric circuit breakers (low and average tension) and have proved to supply the design engineers with much more reliable results than the available penalty-based softwares, see Section 9. In particular, a significant advantage is that a change of topology in the mechanisms does not necessitate a re-estimation of the contact parameters as it is the case when spring-dashpot contact models are used.

### 6.7.3 Baraff's algorithms

Baraff, motivated by problems in computer graphics and animation, essentially focused on the calculation of contact forces and development of specific methods to increase the speed of calculation and the resolution of LCPs or NCPs. Many different aspects of rigid body simulations are discussed in [179], like Painlevé paradoxes, the formulation of friction in 3 dimensions, as well as quite interesting developments on Lemke's algorithm. Lemke's algorithm is used in the early works, but Baraff developed a Dantzig algorithm in [185], See Section 6.6.

### 6.7.4 Other schemes

The works in [64], [162], [163], [168–172], [191], and [193–198] essentially focus on the dynamics formulation (Section 2.1), the collision detection (Section 6.5), the impact rule. The type of integrator that is used for the continuous phases of motion is generally not provided (except [169,196] RK 4, [162] DADS 2D, [198] Newmark). The problems raised in Sections 3.1, 3.2, and 3.3 are usually ignored, except [171] and [193] who analyze multiple impacts.

It is for instance not clear how the capture and finite accumulations of impact times $t_k$ are treated numerically. An extension of DASSL (called DASSLRT) has been used in [199] to simulate several discontinuous phenomena that affect robot motion. Basically these authors consider that the dynamics in Eq. (2) is a time concatenation of DAEs (no mention is made of inequality constraints and related notions). A root finding algorithm is used to compute transition instants. DASSLRT and the ADAMS routine supplied in SIMULINK with a coordinate reduction, see Section 2.2—are compared on an example (a double pendulum with a singular configuration, ie, $\nabla \phi_b(q)$ in Eq. (2) is not full column rank). Other routines failed. DASSLRT seems to be quite fast, $O(0.1 \text{ s})$, while ADAMS took 20 min, due to the singular configuration, in the neighborhood of which acceleration diverges<fnref n="9" />. No numerical result incorporating detachment and collisions is presented. The specific multiple contacts modelling problems and complementarity conditions are almost totally ignored in most of these studies. The study in [151–153] may be a first step towards a better understanding of interconnections between the integrators of ODEs or DAEs and impact detection algorithms, despite it is limited to a very simple one degree-of-freedom system (with however possible finite accumulations of impacts).

## 7 TIME-STEPPING ALGORITHMS

We have described above some various (and classical) manners to discretize DAEs and ODEs. The question one may ask to oneself is: how does the addition of complementarity conditions and state re-initialization modify the properties of such schemes? In other words, many schemes are known to be consistent (ie, the discrete piecewise constant solution converges—in a certain sense—towards the solution of the real dynamics). As we shall see further, many time-stepping schemes, though not all, consist of discretizing simultaneously the continuous dynamics and the complementarity conditions, consequently forming a LCP or a NCP to be solved at each step. The simulation is then easily advanced in time by solving the LCP, using the available solvers, see Section 6.6. Additionally, state re-initializations are needed. The interest of time-stepping methods over event-driven ones, is that they aim at providing a difference-equation approach to the simulation problem, which is suitable for convergence proofs. This may be quite important in case there are many impacts, because it remains to be shown that event-driven schemes are robust with respect to the accumulation of initialization errors (after each impact). This fact combined with the problem of sensitivity with respect to initial data, may render their use quite delicate.

Actually the methods presented in this section all belong to the time-stepping schemes family. But as we shall explain,

<footnote n="9">⁹Notice that this is an artifact. Indeed consider Eq. (35). Clearly, if one tries to compute the multiplier using this formula, difficulties arise whenever the inverted matrix loses its rank. However, as shown in [43] and [71], if $\nabla \phi \in \mathbb{R}^{n \times m}$ has full column rank, then the multiplier vector and the acceleration are unique. If the matrix loses its column rank, then the acceleration remains unique, but not the multiplier vector. Loss of column rank may occur because of hyperstatism or at a singular configuration. Special techniques can be developed to cope with this problem [4,7]. Some authors [83] recommend to avoid inverting matrices of the type $W^T M^{-1} W$.</footnote>
there are significant discrepancies between the schemes within this family. For instance, the ones in Section 7.2 and Section 7.3 do not explicitly require the calculation of the contact forces, contrary to the ones in Sections 7.4 and 7.1. Moreover the schemes in Section 7.2 and Sections 7.3, 7.4 do not rely on an accurate determination of the shock instants: they work with constant h. These schemes are, therefore, true difference equations of the MDI in Eq. (2), with however, possible iterations within a step due to their implicit form. The schemes in Section 7.4 work with h constant, but the shock equations (in some of them) are treated apart when a velocity reinitialization is needed. The schemes in Section 7.1 are based on an accurate calculation of the times of contact with a linear interpolation, similarly to event-driven algorithms. The algorithms in Section 7.2 and Section 7.4 are of the velocity-impulse form. This has some consequences on the integration when friction is present, see Section 8.1.

## 7.1 Lötstedt's algorithms

In [27], Lötstedt introduced time-discretization procedures, at the acceleration/force level. We choose to classify these algorithms into the time-stepping section, because Lötstedt explicitly discretized the dynamics and the complementarity conditions (on the velocity) to form a LCP or a NCP whose unknown is the multiplier $\lambda_i$ at each step $i$. However the discontinuity instants (on the velocity—when there is a shock— and the acceleration—when there is a transition between stick and slip phases—) are computed by an inverse linear interpolation, similarly to an event-driven algorithm. The detection of these instants is made by monitoring the impulse (considered to be zero under a certain threshold as in Section 6.4). Lötstedt's work can with no doubt be considered as an important pioneering work in the field of time-discretization of nonsmooth mechanical systems.

### 7.1.1 The frictionless case

Let us consider first the frictionless case. The following numerical scheme is proposed to compute the state at step $i$:

$$ egin{cases} q_i = \frac{1}{\alpha_0^1} b_i^1 \\ \dot{q}_i = \frac{1}{\alpha_0^2} [h \beta_0^2 M^{-1} (Q_i + \nabla \phi_i \lambda_i) + b_i^2] \\ \nabla \phi_i^T \dot{q}_i \ge 0, \quad \lambda_i \ge 0, \quad \lambda_i^T \nabla \phi_i^T \dot{q}_i = 0 \end{cases} \quad (49) $$

where

$$ b_i^1 = h \sum_{k=0}^{r} \beta_k^1 \dot{q}_{i-k} - \sum_{k=1}^{r} \alpha_k^1 q_{i-k} \quad (50) $$

$$ b_i^2 = h \sum_{k=1}^{r} \beta_k^2 M^{-1} (Q_{i-k} + \nabla \phi_{i-k} \lambda_{i-k}) - \sum_{k=1}^{r} \alpha_k^2 \dot{q}_{i-k}. $$

Clearly, the complementarity relations in Eq. (49) correspond to the active constraints at step $i$, see Section 3.1. They encompass the persistent contact as well as plastic impacts phases. The formulas in Eqs. (49) and (50) correspond to two linear $r$-step methods [25] Section VII.3. $\nabla \phi_i = \nabla \phi(q_i)$ and

$\phi_i = \phi(q_i)$. The coefficients $\alpha_k^1$ and $\beta_k^1$ are determined from an Adams-Bashforth family of explicit formulas ([26], page 250), denoted as AB-r, $\beta_0^2=1$, $\beta_k^2=0$ for all $k=1,2,\dots,r$. The second equation in (49) is a backward difference formula, denoted BDF-r. Notice that the mass matrix $M$ is assumed to be constant (hence the Coriolis and centrifugal torques are zero), which restricts the application to simple mechanical systems with Euclidean configuration space (like collections of particles). It is however argued that this is just a matter of convenience to allow for an easy factorization of $M$, and that the extension towards $M(q)$ is possible ([27], Section 6). The torque $Q_i = Q(q_i, \dot{q}_i, t_i)$ therefore contains gravity, viscous friction, and external actions (like control inputs). The integration step is chosen constant, equal to $h$.

When $Q=Q(t,q)$, it is shown in [27] that a LCP whose unknown is $\lambda_i$ can be formulated from (49). As we already pointed out, this LCP can be rephrased as a quadratic program:

$$ \min_{\lambda_i \ge 0} \frac{1}{2} \lambda_i^T \nabla \phi_i^T M^{-1} \nabla \phi_i \lambda_i + h^{-1} \lambda_i^T \nabla \phi_i^T [b_i^2 + h M^{-1} Q_i]. \quad (51) $$

Consequently, the set of equations in (49) allows one to advance the solution in time from $i-1$ to $i$. The methods AB-1 (forward Euler)-BDF-1 and AB-2-BDF-2 are chosen in [27], where it is recalled that it is useless to use methods of order $\ge 3$ (linear multistep A-stable methods have an accuracy of order $\le 2$, ie, at most $O(h^2)$ [26] pages 250–251). After discontinuities in $\dot{q}_i$ or $\ddot{q}_i$ (which are detected from the value of the impulse on one step-with a threshold under which it is considered to be zero), the AB-1-BDF-1 algorithm is used during two steps to restart the simulation (it is known that multistep methods are not self-starting and require the help of a single-step algorithm initially).

When $Q=Q(q,q,t)$, then the LCP formulation is lost. However Lötstedt proves that provided the matrix

$$ A(\dot{q}) = M - \frac{h}{\alpha_0^2} \frac{\partial Q}{\partial \dot{q}}(q, \dot{q}, t) \quad (52) $$

is full rank and $\nabla \phi_i^T A^{-1} \nabla \phi_i$ is positive definite, then (49) still possesses a unique solution so that the algorithm can be used to safely advance the solution in time. However this time $\lambda_i$ is generally the solution of a NCP (a quick look at the second equation in (49) allows one to realize this). The condition in (52) can be used with the implicit function theorem to express $\dot{q}_i = g_i(\lambda_i)$ for some function $g_i(\cdot)$. The second condition is used to prove the existence of a solution to the NCP. A way to solve the NCP is proposed, based on functional iteration. Certainly this could be improved using new tools, see Sections 3.1 and 6.6. In summary, Lötstedt's algorithm is given as follows [27]:

* Compute $q_i$ using AB-1 or AB-2, with $h$ such that the local error in $q_i$ is smaller than $h \epsilon$ for a prescribed tolerance $\epsilon$ (ways to estimate such a $h$ are given in [27]).

* Calculate $\nabla \phi_i$ to a prescribed accuracy and calculate $\ddot{q}_i = M^{-1}(Q_i + \nabla \phi_i)$ and $\dot{q}_i$ by BDF-1 or BDF-2.
* Test whether velocities and accelerations are discontinuous between $t_{i-1}$ and $t_i$, either due to an impact (detected from a nonzero value of the impulse) or to the activation of a new constraint ($\phi_k(t_{i-1})>0$ and $\phi_k(t_i)\le 0$ for some $k$), or to the deactivation of a constraint. The time of such jumps is calculated by inverse linear interpolation. After a shock a new velocity $\dot{q}_{i+1}$ is computed by a collision rule (the rule in Eq. (12) is used). Then restart the algorithm at the first step with AB-1 and the new set of active constraints.

* Test the detachment conditions by checking whether one entry of the vector $\lambda_i$ passes through zero, and whether the corresponding entry in the normal velocity $\nabla \phi_i^T \dot{q}_i$ is positive. Then refresh the set of inactive constraints if needed.

* End.

**Remark 7:**

* The algorithm is based on velocity constraints to reduce the index and form an LCP. The drawback is evidently a possible drift away from $\partial K$ during permanent contact phases (DAE simulation). The stabilization on the constraints $\partial K$ during permanent contact phases needs special attention [22,66], see Section 5.

* Lötstedt also shows that the LCP($\lambda_i$) can be reformulated as the minimization problem

$$ \min \left\| \sum_{k=0}^{r} \alpha_k^2 \dot{q}_{i-k} - h M^{-1} Q_i \right\|_M, \quad \nabla \phi_i^T \dot{q}_i \ge 0. \quad (53) $$

* The velocity can be calculated directly from Eq. (53). However it is argued in [27] that it is better to compute the multiplier and then insert it in the dynamics because this simplifies the initialization of the algorithm (a feasible $\lambda_0$ is easy to find), the QP in Eq. (51) is a particular case of the QP constructed for the case with friction (see below with comments), and it may be interesting to know the value of $\lambda_i$ explicitly, since $\lambda_i$ can be used as a measure of force. This is of interest for control design of real-world systems [21].

* Dissipativity of Eqs. (49), (50), and (51) plus the impact rule and convergence of the algorithm are not proved.

### 7.1.2 Constraints with 2D friction

As recalled briefly in the introduction and Section 8, the dynamics in Eq. (2) is much more complicated when friction is considered. Friction may create some unexpected phenomena [18,96,200,201] as the divergence of the contact force (but with bounded impulses) or so-called *Impacts without Collisions*. These phenomena are not to be considered as artefacts due to the model deficiencies. Although they disappear when rigidity is relaxed, or when the Coulomb model is replaced by some regularized law, they really represent physical phenomena which occur in real systems. More comments are given in Section 8. Since Lötstedt's algorithms are acceleration/force schemes, these phenomena should be taken into account. Lötstedt was perfectly aware of such problems [202], and therefore, proposed a particular numerical procedure to avoid them.

The same algorithms AB-1-BDF-1 or AB-2-BDF-2 are used as in the frictionless case. When 2D friction is incorpo-

rated in the algorithm, one has to add the tangential contribution of the contact force in the right-hand-side of the second equation in (49). The contact force is split into two parts, see Section 3.3: $\nabla \phi \lambda + P(q, \dot{q}, \lambda) = [G_1(q) + H(q, \dot{q})]\lambda_1$. Roughly $G_1(q)\lambda_1$ contains the normal generalized force and the contribution of the sticking contacts, whereas $H(q, \dot{q})\lambda_1$ accounts for the sliding contacts. The vector $\lambda_1$ contains the normal multipliers $\lambda_{n,j}$ and the tangential ones $\lambda_{t,j} = \pm f\lambda_{n,j}$ for the $j$th contact. There are two features in the algorithm proposed in [27]. The first one is the approximation of $\lambda_{1,i}$, the second one is the calculation of the impulses at the shock instants. Let us denote the $j$th component of $\lambda_1$ by $\lambda_1^j$ and its $i$th iteration by $\lambda_{1,i}^j$. Then the approximated value is $\bar{\lambda}_{1,i}^j = \lambda_{1,i-1}^j + h_i(\lambda_{1,i-1}^j - \lambda_{1,i-2}^j)/h_{i-1}$, for a variable step of integration $h_i$. A QP is constructed that allows the computation of the term $G_1(q)\lambda_1$. As noted in [185] (Section 9.2.1), it possesses the advantage of assuring that the tangential force is opposite the tangential acceleration. But it has the strong drawback that sliding generally implies the QP matrix to be non-symmetric, rendering the problem harder to solve. It is clear that the introduction of $\bar{\lambda}_{1,i}^j$ in the dynamical equations modifies the subsequent calculations in a non-physical manner right after the first step$^{10}$, and should be avoided. Special procedures are also used after a shock and a discontinuity in the acceleration. The error introduced in $\dot{q}_i$ by the use of $\bar{\lambda}_{1,i}^j$ in a permanent contact phase are shown to be $O(h^3)$ when $h_i=h$, a constant. They are $O(h)$ after a reinitialization of the velocity or of the acceleration.

The second point (calculation of the impulse at a shock instant), is formulated as follows. Taking frictional effects at impacts into account, let us denote the right-hand-side of Eq. (3) as $P_k=G(q)\Lambda_k$, where $\Lambda_k$ is a vector of normal and tangential percussions. Then Lötstedt proposes to calculate the impulse from the QP in Eq. (17). The dual version of this QP is:

$$ \begin{align} & \min \frac{1}{2} [\dot{q}(t_k^+) - \dot{q}(t_k^-)]^T M [\dot{q}(t_k^+) - \dot{q}(t_k^-)] \nonumber \\ & W^T v = G^T \dot{q}(t_k^+), \quad v \ge 0, \quad v^T W \Lambda_k = \Lambda_k^T G^T \dot{q}(t_k^+) = 0 \nonumber \\ & W = \begin{pmatrix} I & 0 \\ f_{imp} I & -I \\ f_{imp} I & I \end{pmatrix} \tag{54} \end{align} $$

where $I$ is the identity matrix with dimension equal to the number of active constraints, and $f_{imp}$ can be considered as an impulse ratio [93]. The main problem with the calculation in Eq. (54) is that although it looks like the plastic impact rule Eq. (12), it is not like (12): there may be rebounds. In addition, if there is a tangential velocity reversal during the shock (ie, the post and pre-impact tangential velocities have opposite signs), then there may be a kinetic energy gain at the shock instant (this phenomenon is well-known in the literature, see eg, [18,93]).

<footnote>$^{10}$Indeed the state at step 2 is calculated with $\bar{\lambda}_1^j$, so the next $\lambda_1^j$ is calculated from wrong positions and velocities.</footnote>
**Remark 8:**

• Lötstedt discusses the issue of Section 6.5 for collections of polyhedral objects. He points out the need for a selection procedure for the collision detection, but does not provide many details on this part of his algorithm.

• The algorithm in [164] is close in spirit to Lötstedt's one (time-stepping with accurate detection of contacting times). It uses a trapezoidal discretization of the continuous frictionless dynamics (implicit one-step scheme, solved by a Newton method with an initial guess from a Euler's discretization), and an adaptive step size procedure. Several simple examples show that $h$ may decrease to very small values as $10^{-12}$ s during the simulation. Lemke's algorithm is used to solve the contact force LCP and the impact percussion LCP (see Eqs. (11) and (12)).

## 7.2 Discretization of Moreau's sweeping process

### 7.2.1 Frictionless sweeping process

We describe in this section and the next one the time-discretization of the general MDI presented in Sections 3.4.1 and 3.4.2. It has been named the *NonSmooth Contact Dynamics* (NSCD) method by Moreau and Jean, and can be considered as one of the results of the research led by Moreau in Montpellier on Convex Analysis and Nonsmooth systems since the early 60s (see eg, [70,95,104,203,204]). The presented methods constitute the first attempt [28] of time-discretization of MDIs as in (22), ie, simulation of multibody systems without regularization of either the normal or the tangential friction laws of contact/impact. The Contact Dynamics method provides a very general and powerful framework for the simulation of various nonsmooth mechanisms, including granular matter [12,78,205], buildings made of blocks and monuments [13], deep drawing process [94,206,207], robotic systems [208] and kinematic chains [5]. It so happens that the time-discretization of the MDI in Eq. (22) yields an intrinsic implicit formulation. But it can be transformed into an explicit scheme using basic convex analysis [101,209]. We first concentrate on the discretization of (22), then we indicate how friction is treated, and finally we focus on the general NSCD method. Following [28,96] choose in Eq. (22) at step $i$:

$$ egin{cases} t_{m,i} = t_i + \frac{1}{2}h \\ q_{m,i} = q_i + \frac{1}{2}hv_i \\ M_{m,i} = M(q_{m,i}) \\ Q_{m,i} = Q(t_{m,i}, q_{m,i}, v_i) \\ v_{l,i} = v_i + hM_{m,i}^{-1}Q_{m,i} \\ -M_{m,i}(v_{i+1} - v_{l,i}) \in \partial\psi_{V(q_{m,i})}(v_{i+1}) \\ q_{i+1} = q_{m,i} + \frac{1}{2}hv_{i+1} \end{cases} \quad (55) $$

where $h$ is the integration step (that might be varying from one step to the next). One sees that the computation of the values at step $i+1$ depend on intermediate calculations at the midpoint $t_{m,i}$. The midpoint velocity $v_{l,i}$ is equal to the approximated velocity that the system would have at step $i$

+1, if there was no collision on $(t_i, t_{i+1})$. This is easily seen from the equivalence between the inclusion in Eq. (55) and

$$ v_{i+1} = \text{prox}_{M_{m,i}}[V(q_{m,i}), v_{l,i}] \quad (56) $$

where prox means the proximation operation in the kinetic metric (with an underlying projection in the same metric). The tangent cone $V(q(t_i))$ is defined also outside the admissible domain $K$ in order to cope with possible violation of the constraints in the course of the simulation, see Eq. (18). Indeed if there is a contact event between $i$ and $i+1$, which is detected by checking the negative sign of $\phi(q_{m,i})$, one uses Eq. (56) to compute $v_{i+1}$. This is a quadratic program in the kinetic metric $M_{m,i}$. The reader may recognize again that the formulation in the last line of Eq. (55) encompasses the whole dynamics (continuous motion and shocks). This discrete-time inclusion is a discretization of the MDI in Eq. (22) with $dv \approx v_{i+1} - v_{l,i}$, $v(t^+) = v_{i+1}$ and $v_{l,i}$ can be interpreted as a left-velocity at $t_{i+1}$.

The general case $0 \le e \le 1$ can also be handled, as we observed in Section 3.4.1 [35,96]. Let $u = v(t^+) + ev(t^-)/1 + e$ in Eq. (22) (with $v(t^+)$ in the right-hand-side of (22) also replaced by $u$, see also Remark 4). Then Eq. (22) can be discretized as follows [35]:

$$ \frac{h}{1+e} Q(t_{i+1}, q_{i+1}, v_{i+1}) - M(q_{i+1}) \frac{v_{i+1} - v_i}{1+e} \\ \in \partial\psi_{V(q_{i+1})} \left( \frac{v_{i+1} + ev_i}{1+e} \right) \quad (57) $$

which is clearly an implicit formulation but is equivalent to:

$$ \frac{v_{i+1} + ev_i}{1 + e} = \text{prox} \left( V(q_{i+1}), v_i + \frac{h}{1 + e} Q(t_{i+1}, q_{i+1}, v_{i+1}) \right). \quad (58) $$

One notes that there are some differences between Eqs. (57), (58), and (55). This last formulation is used in [32] and [35] for the study of convergence of the algorithm. This indicates that mathematical convergence proofs and real implementation of an algorithm may lead to different discretizations. Equation (58) is an explicit form as long as $Q = Q(t, q)$ and provided one chooses $q_{i+1} = q_i + hv_i$. Indeed in this case introducing $q_{i+1}$ inside the left hand side of Eq. (58) one gets a direct way to calculate $v_{i+1}$ by a proximation tool (ie a quadratic programming approach).

**Remark 9:**

• **Convergence:** Proofs of convergence of the discretized Moreau's sweeping process can be found in [32], using a discretization as in Eq. (57) and $e=0$. The scheme in Eqs. (57) and (58) is proved to converge globally for $0 \le e \le 1$ in [35] using techniques inspired from [32], with $\partial K \in C^{1,\beta}$, $\beta > \frac{1}{2}$ and codimension one constraints ($m=1$).

• **Impact calculation:** The usefulness of the midpoint calculation in Eq. (55) lies for instance in the fact that, as one easily calculates, $q_{i+1} = q_i + \frac{1}{2}h(v_i + v_{i+1})$. Contact is detected if $\phi(q_{m,i}) < 0$, while $\phi(q_i) > 0$. In the case of an elastic collision (think for instance of the simple one-
degree-of-freedom bouncing ball), one gets $v_{i+1} = -v_i$, since the velocity is simply reversed as can be checked from Eqs. (55) and (56): one has $\text{prox}[R^+, v_i + h/2g] = 0$ for $h$ small enough since $v_i < 0$ at a shock instant; modify (56) in the same lines as Eq. (58) to get the result. Thus $q_{i+1} = q_i$, hence the ball rebounds exactly (at the machine accuracy) to the same height, whatever the length of the simulation. This is also pointed out in [91] who discusses symplectic integrators and energy conservation problems on the bouncing ball example. As we shall see some schemes do not possess this property, and may yield energetic inconsistencies. One sees that in the time-stepping scheme in Eqs. (55)-(56), $v_{i,l}$ plays the role of the pre-impact velocity, while $v_{i+1}$ plays the role of the post-impact velocity. This is on contrast with some event-driven schemes, where the impact rule is applied at step $i$. The midpoint Euler scheme is further justified in [28] (Section 6), by the fact that it assures a much more accurate estimate of $q_i$ than a Euler algorithm. This fact combined with the need of jumps detection led Moreau to choose a discretization as in Eq. (55). Notice that (55) is a one-step method, but is not a second-order explicit RK ([26], page 247). Further sophistications of the algorithm are presented below.

*   **Constraint drift:** When there is only one contact that is made, the proximation in Eq. (56) is equivalent to $\nabla \phi_{m,i}^T v_{i+1} = 0$. This will generally result in a bad stabilization of the constraints during persistent contact phases. The approximation of the tangent cone $V(q_{m,i})$ can be replaced by $V(q_i + h v_i)$. This may help in satisfying the constraint better [96]. In addition if the constraints are not respected, then Moreau [28] proposes to use a projection of $q_{i+1}$, denoted as $q_{i+1}^c$, onto $\partial K$, computed in the kinetic metric as:

$$q_{i+1}^c = q_{i+1} - \phi(q_{i+1})(\nabla\phi M^{-1}\nabla\phi^T)M^{-1}\nabla\phi \quad (59)$$

where the last term is evaluated at $q_{i+1}$.

*   **Transition phases and capture:** As we explain in the Impact rules paragraph of Section 7.4, the transition between free-motion and persistent-contact phases—ie, the problem of capture—, via a shock or a sequence of shocks, is automatically treated by the algorithm in Eqs. (55) and (56) or (57) and (58). This feature is actually shared by other time-stepping schemes as the ones in Section 7.3, and constitute a serious advantage of time stepped algorithms (which are a truly difference equation approach to simulation, contrary to event-driven schemes which are in a sense closer to the continuous-time description for impact detection and collision effects). Note however that in event-driven schemes contact/rebound happens at one time instance $t_i$, whereas with time-stepping schemes there is always at least a time-delay $h$ between contact and release. Incidentally, notice that replacing the left-hand-side in Eq. (56) by $v_{i+1} + ev_i/1 + e$ does not change much the capture problem once $v_i \approx 0$. Here again one sees that $v_{i+1}$ plays the role of $v(t^+)$ whereas $v_i$ plays the role of $v(t^-)$ when an impact is detected at $t_{m,i}$.

*   Actually as shown by Moreau [100,101] the first-order

sweeping process formulation $-dv/dt \in \partial\psi_{C(t)}(v(t))$ with $C(t)$ a moving convex set, cannot be discretized via any Euler explicit scheme of the form $v_i - v_{i+1} \in \partial\psi_{C(t_i)}(v_i)$ since this is equivalent to $v_i = \text{prox}(C(t_i), v_{i+1})$ which is nonsense. On the other hand, $v_i - v_{i+1} \in \partial\psi_{C(t_{i+1})}(v_{i+1})$ makes sense since this is equivalent to $v_{i+1} = \text{prox}(C(t_{i+1}), v_i)$ which provides an explicit way to calculate $v_{i+1}$. One notices that the integration step does not appear in these expressions: this is due to the fact that since the right hand side is a cone, multiplying it by any positive constant does not change the inclusion. In the second-order sweeping process formulation, $h$ appears explicitly at the velocity level, see Eq. (57).

### 7.2.2 Constraints with Coulomb friction

At each contact $A^{(j)}$ the Coulomb friction law is formulated as [94,104]:

$$F_{t,i+1}^{(j)} \in \mathcal{D}_{i+1}^{(j)}, \quad \forall s \in \mathcal{D}_{i+1}^{(j)}, \quad [s - F_{t,i+1}^{(j)}]^T V_{t,i+1}^{(j)} \ge 0 \quad (60)$$

where the notations are the same as in Section 3.4.2. In particular, $\mathcal{D}_i^{(j)} = \mathcal{D}(f|N_i^{(j)})$ has an obvious meaning from the definition of $\mathcal{D}$. The discrete form in Eq. (60) is equivalent to $F_{t,i+1}^{(j)} = \arg\max_{z \in \mathcal{D}_{i+1}^{(j)}} (-z^T V_{t,i+1}^{(j)})$, and to $F_{i+1}^{(j)} = \arg\max_{z \in \mathcal{C}_{i+1}^{(j)}} (-z^T V_{i+1}^{(j)})$. In the 2D case, Eq. (60) exactly represents the Coulomb graph between $V_{t,i}^{(j)}$ and $F_{t,i}^{(j)}/N_i^{(j)}$. In the 3D case, one recovers a friction cone.

The Gauss-Seidel method [78,101] roughly consists of the following Signorini Coulomb loop [94]. Suppose that for contact $l$ at step $i$ the data $F_{t,i}^{(j)}, N_i^{(j)}, V_{t,i}^{(j)}, V_{n,i}^{(j)}$ are known for all $j > l$ from the previous iteration, and from the current iteration for all $j < l$. Then compute the status of contact $l$ by solving the Signorini Coulomb conditions (which monitor the transitions at $A^{(l)}$: sticking/sliding/detachment). Iterate until the last contact. Then apply a convergence test, that may be on the contact forces, ie, $|F_i^{k+1} - F_i^k|/|F_i^k| \le \epsilon$, where $k$ is the index of the iteration of the Signorini Coulomb loop at step $i$, or on the distance to the Signorini graph. If it is not satisfied redo the calculations for all contacts. One may also choose to stop arbitrarily after a certain number of iterations, since in some practical cases convergence is not at all guaranteed and the algorithm could stuck at one step. The fact that the contacts are treated independently one after each other, guarantees at least that Coulomb friction law is satisfied by $V_{t,i+1}^{(j)}$ and $F_{t,i+1}^{(j)}$ at each contact at the end of the iteration process. This is not necessarily the case for other formulations, see Section 7.4. Such a cyclic procedure is similar to a nonlinear block Gauss Seidel algorithm (which is for instance used to solve QPs [5]). Some convergence results can be found in [113].

**Remark 10:**

*   The dynamics is written from the Lagrange equations as in (2) whereas the dry friction law involves local quan-
tities at the contact points $A_j$ as in (60). Consequently one has to calculate the Jacobians $J^{(j)}(q)$ to relate both (see Section 3.4.2). It is argued in [94] that $J^{(j)}(q)$ may be evaluated at various intermediate values of $q$, eg, $q_i$, or $q_i + h/2\dot{q}_i$, or even a value obtained from an iteration process to approximate $q_{i+1}$ (in case of an implicit formulation of the dynamics one needs to implement a Newton-Raphson like algorithm to compute $q_{i+1}$ and $\dot{q}_{i+1}$ at step $i$). This is however related to the curvature of the contacting surfaces, which regulates the change in $\nabla\phi$ from $i$ to $i+1$.

* A discretization of the sweeping process with friction that generalizes that in [32] is proposed in [210], for $m > 1$. The inclusion (23) is merely written at each contact point.

* The sweeping process (second order) discretization does not resort to any calculation of the acceleration and the contact forces. Only velocities (and implicitly the contact impulses) are calculated at each time step. This allows one to avoid the problems related to Painlevé's paradoxes when friction is present (see Sections 7.2.2 and 8.1), ie, the computation of unbounded contact forces (for certain friction coefficients and configurations of the system). Indeed as shown in [200,201] on a simple example, the force impulse remains bounded when the orbit crosses so-called critical points. The so-called impacts without collisions, or tangential impacts (which are of a different nature from the contact force unboundedness problem) are handled via the principle of maximal dissipation: at each time step, the tangential impulse is calculated so that it maximizes the dissipated frictional energy.

### 7.2.3 Simulation results

Extensive numerical tests have been performed by Moreau on granular matter [78,205]. Simulation of granular matter (sand piles, planetary rings) is a difficult subject, essentially due to the very large number of degrees of freedom. It is also difficult to make experiments and to compare them with simulations (try to follow the motion of a sand grain in a sand pile!). Only macroscopic phenomena may be checked (resulting pressures, average stress tensors, distribution of reaction forces), and may be expected to be robust with respect to numerics [94]. It is even difficult to make numerical simulations, and compliant models may simply be impossible to use, see Section 4.2. On the other hand, this is a topic with major applications in industry and in theoretical physics. In [78], a 2D vertically shaken cylindrical vessel is simulated. It contains 3999 beads with diameter 0.2 cm and one bead with diameter 0.5 cm (ie, $n=4000$ in (2)). In [205], another test is made with 2000 beads with diameter 0.2 cm and 200 beads with diameter 0.02 cm. A 3-parameter impact law is chosen ($e_n, e_t, f$), where $e_t$ is a tangential restitution coefficient [93]. The numerical integration allows one to test the influence of the physical parameters on the global behavior, and to verify if some phenomena like clusterization, fluidization, bulk segregation, convection effects, occur and under which conditions.

**Remark 11:**

The dynamical effects can be quite important in the overall

motion even for almost static systems, because they permit to simulate possibly sparse events that have a crucial influence on the final configuration.

### 7.2.4 Related algorithms

**NSCD methods:** The work in [116], that is part of the NSCD method, is dedicated to the simulation of a continuous medium, approximated by a mesh of nodes. This is called the *distinct elements method*, in order not to confuse with the classical Finite Element Method (FEM). The nodes collide with the obstacle while the whole structure deforms. The inertia matrix is constant. Various discretization procedures based on $\theta$- and $\theta$-Euler methods are compared. Let us recall that the $\theta$-method yields $\int_{t_i}^{t_{i+1}} f(t)dt \approx h[\theta f(t_{i+1}) + (1 - \theta)f(t_i)]$. The $\theta$-Euler method consists of discretizing the Lagrange equation in (2) by a $\theta$-method, and the position by a Euler method. A modified $\theta$-method is also tested which consists of using a $\theta$-method, then replacing $q_{i+1}$ by $q_{i+1} + h(1 - \theta)\dot{q}_{i+1}$ in the contact relations. In particular the influence of the discretization of the complementarity conditions (second line in Eq. (2)) is discussed in [116]. The term $\phi^{(j)}(q(t))$, which represents the local distance between node $j$ and the constraint, is discretized at each node as (see Remark 12):

$$ \begin{aligned} \phi_{i+1}^{(j)} = & \phi_i^{(j)} + h\theta[\nabla\phi_{i+1}^{(j)}]^T\dot{q}_{i+1} + h(1-\theta) \\ & \times [\nabla\phi_i^{(j)}]^T\dot{q}_i, \quad (\theta\text{-method}) \\ \phi_{i+1}^{(j)} = & \phi_i^{(j)} + h[\nabla\phi_{i+1}^{(j)}]^T\dot{q}_{i+1}, \end{aligned} \quad (\theta\text{-Euler and modified } \theta\text{-methods}). \quad (61) $$

The reader will see that contrary to Eqs. (49), (50), and (69), (70) which base on the analytical form of $[\nabla\phi^{(j)}(q)]^T\dot{q}$ and on the calculation of this expression using $q_i$ or $q_{i+1}$ and $\dot{q}_{i+1}$ or $\dot{q}_i$, the expressions in Eq. (61) are a time-discretization of $\phi^{(j)}(t)$. As shown in [155] this may have important consequences on the numerical results (eg, concerning the calculation of Lyapunov exponents). For the $\theta$ and $\theta$-Euler methods, the complementarity conditions are at node $j$: $\phi_{i+1}^{(j)} \ge 0$, $N_{i+1}^{(j)} \ge 0$, $\phi_{i+1}^{(j)} N_{i+1}^{(j)} = 0$. For the modified $\theta$-method, they are formulated with $\bar{\phi}_{i+1}^{(j)} = \phi_{i+1}^{(j)} + h(1 - \theta)\dot{\phi}_{i+1}^{(j)}$. It is noted that for the $\theta$-method, $\phi_i^{(j)} = 0$ does not imply that $[\nabla\phi_{i+1}^{(j)}]^T\dot{q}_{i+1} = 0$, but this is the case for the $\theta$-Euler and modified $\theta$ methods. In [94], it is pointed out that the position complementarity conditions with $\bar{\phi}_{i+1}^{(j)}$ are not recommended for large collections of bodies (like granular matter with $n \ge 1000$), because the correcting impulse that is calculated after a penetration to send back the system in $\text{Int}(K)$ may be a non-negligible numerical artifact. Then a velocity complementarity formulation is preferable.

**Remark 12:** The idea of consistency of the gap approximation with unilateral condition is introduced in [94]. It means that $\bar{\phi}_i^{(j)} = 0$ and $\bar{\phi}_{i+1}^{(j)} = 0 \Rightarrow (\nabla\phi_{i+1}^{(j)})^T\dot{q}_{i+1} = 0$ is satisfied. For instance if one chooses simply $\bar{\phi} = \phi$, this consistency is satisfied only when $\theta = 1$. Other choices are possible
as the ones above. When this consistency property is not satisfied, some artificial numerical oscillations may appear [116]. The choice of the gap approximation is not important, provided it has the consistency property.

For the 2D case, dry friction is transformed in [116] via velocity slack variables into a set of two complementarity relations per contact, see Eq. (14). Then the whole stuff is assembled as for the time-stepping schemes in Section 7.4 (a technique already used in [114]). Lemke's algorithm is used to solve the LCP in [116]. Inelastic impacts with friction are treated in [94] by combining the resolution of complementarity conditions between the normal quantities $\bar{\Phi}_i^{(j)}/h + V_{n,i}^{(j)}$ and $N_{i+1}^{(j)}$, and the resolution of Eq. (60). In a compact form, they might be written as $\text{SignCoul}(i, V_{i+1}^{(j)}, F_{i+1}^{(j)})$. Using the second formulation in Eq. (21) it may be shown that the mapping $(V, F) \mapsto \text{SignCoul}(i, V, F)$ is piecewise affine in 2D, and continuous raywise in 3D.

**Simulation results:** Due to its implicit formulation and the Signorini Coulomb loop, the NSCD method consists of two nested iteration loops at each time step $i$ [94]. Although large integration steps can be used, it is therefore time consuming, which renders its use for real-time applications less easy. The numerical results obtained by the three discretization methods above, are compared to the exact solution in the case of the impact of two identical elastic bars moving on a line. They are further compared to a Newmark method with $\gamma=0.5$ and $\beta=0.25$ (these values are the smallest ones that guarantee A-stability [26] Ch 7.1.5). The Newmark and $\theta$-method with $\theta=0.5$ generate oscillations during the contact phase. Taking $\theta=0.55$ damps out the oscillations (increasing $\theta$ renders the algorithm more implicit, and it is known [94] that implicit schemes behave nicely when applied to nonsmooth problems, although their implementation is more complex). The transition contact-non contact is too slow with the $\theta$-Euler method. So the modified $\theta$-method is the best one for this 1D case. For more complex case studies the $\theta$-Euler method is said to be a good compromise between the $\theta$-method (too smooth velocity) and the modified $\theta$-method (oscillations). A 2D example from [211] (a disc colliding a rigid ground) is simulated with the $\theta$-method ($\theta=0.55$). The mesh has 99 nodes. The oscillations obtained with a Newmark scheme and an impact detection procedure as in [211], and which cannot be explained by acoustic waves propagation, are damped out. Evidently in these cases only a comparison with experimental results would really be telling (although the validation of a code cannot be made by comparison with experimental results [115], but by comparing various numerical schemes tested on benchmarks, or by comparing the discrete solution with the analytical one in simple cases). In [13], buildings made of blocks are simulated (an arch under various loads). Each block is approximated as in Fig. 3. Comparisons with experimental results with wood blocks led in the Ecole Supérieure de Mécanique de Marseille (F) are encouraging. It is important in this setting to recall that Coulomb friction is the main effect that precludes the existence of a unique equilibrium point (despite its dissipativity, see [18] pages 207-208), so that the obtained state depends on the history of loading. Physical

situations are like this. The simulation of a Couette granular flow with 2400, 4000 and 16000 polydisperse disks can be found in [94]. An implicit Euler scheme has been used ($\theta=1$ in Eq. (61)). The NSCD method seems to be particularly suitable for granular matter simulation and has been also used to study stress transmission and granular pressure [12]. Different quasistatic examples are presented in [115]: a dove-tail assembly (4345 nodes, 49 contact nodes), a block sliding on a plane (4193 nodes, 65 contact nodes), and a pressure vessel (674 nodes, 54 contact nodes), a high pressure screw press (11933 nodes, 250 contact nodes), the extrusion of an aluminum cylinder (105 nodes, 21 contact nodes). Emphasis is put on the efficiency of various LCP solvers in terms of CPU times. Some preliminary results have been obtained for kinematic chains [208]. However such systems differ a lot from granular matter and distinct element systems, so that other algorithms (more accurate, with explicit contact forces and events calculation) may be preferred, see Section 6.7.2. The Contact Dynamics method has been implemented in a software called Simem3 [206].

**3-parameter impact law:** A numerical scheme inspired from [78] and [96] has been proposed in [212]. It applies to a 2D rigid body hitting a wall. Its focus is on the proper use of LCPs for impacts with friction and tangential restitution, so that 3 parameters are used (the 3-parameter impact law $(f, e_n, e_t)$ is one of the most used in the literature [18] for impacts of spheres). A contact with Coulomb friction is treated with two LCPs, see Eq. (14). The LCPs are solved by a pivot algorithm. Experimental results of a rubber ball thrown in a box with wooden walls and steel plate are reported. The real motion of the ball between two impacts is recorded by a high-speed camera. The coefficients $f$ and $e_n$ were measured off-line and $e_t$ fitted afterwards to get a good matching between experiments and simulations. The choice $e_n=e_t$ (which assures some energetical consistency [7,84]) is made in most experiments, although $e_t$ is sometimes varied to get a better result. Further comparisons are made with experimental results available in the literature. The conclusions are mitigated and it is pointed out that a more complex impact model including moment impulse [93] could improve the results. The work is extended to 3D problems in [213], relying on an impact rule proposed by Moreau [78]. Roughly, it uses a variational formulation of this impact rule and of Coulomb friction (similar to the variational inequality in (21) and to (60)) to express the dynamics as a nondifferentiable equation, inspired by the works for quasi-static frictional contact in [214]. The problem is solved using a Newton-type algorithm specially devised for such nondifferentiable equations [215]. We note that such a procedure to formulate the complementarity relations from Eq. (21) is quite similar to what is used in [92,182] for the acceleration calculation. Similar experiments as in [212] are made and compared to the numerical results, showing good agreement when $e_t=e_n$.
## 7.3 Discretized second-order MDI

A numerical scheme especially suited to the second-order formulation presented in Section 3.4.3 is proposed by Paoli and Schatzman in [54-56,216,217].

### 7.3.1 Description of the scheme

The scheme is written in terms of positions only. So the relevant local metric is defined locally by the mass matrix of the system. Let us use the same notations as in Section 3.4.3. If we assume that $\partial K$ and $q \to M(q)$ are smooth (the second property being satisfied in most cases), it is possible to define locally a projection on $K$ relatively to the Riemannian metric defined by $M(q)$.

Let $F$ be a continuous function, consistent with $Q$, ie,

$$F(t,q,q,v,0)=M^{-1}(q)\bar{Q}(q,M(q)v,t) \quad \forall (q,v,t).$$

The scheme is given by:

$$q_{i+1} = -eq_{i-1} + (1+e)P_K \left( \frac{2q_i - (1-e)q_{i-1} + h^2 F_i}{1+e} \right) \quad (62)$$

$$F_i = F \left( t_i, q_i, q_{i-1}, \frac{q_{i+1} - q_{i-1}}{2h}, h \right)$$

provided that $q_{i+1}$ is uniquely defined in a neighborhood of $K$. The projection $P_K$ is done in the kinetic metric.

In order to understand how this scheme approximates the system (26)-(29), we must say a few words about its construction. Let us consider the simplest formulation which corresponds to the case of convex constraints with a constant mass matrix. The dynamics in Eq. (33) is approximated by the implicit algorithm

$$\frac{q_{i+1} - 2q_i + q_{i-1}}{h^2} + \partial \psi_K \left( \frac{q_{i+1} + eq_{i-1}}{1+e} \right) \ni F_i. \quad (63)$$

**Remark 13:**

The reader may notice the similarity between Eq. (63) and the last equation in (55). However this time the positions are involved, not the velocities. Notice that $F_i$ is premultiplied by $h^2$ in Eq. (62) whereas the same term is multiplied by $h$ in Eq. (55).

We can transform Eq. (63) by using a result of convex analysis. It can be shown that if $K$ is a nonempty closed subset of $\mathbb{R}^d$ ($d \ge 1$), then with $y \in \mathbb{R}^d$ and $\lambda > 0$ given, the equation

$$x + \lambda \partial \psi_K(x) \ni y$$

admits a unique solution $x = P_K(y)$. Taking $x = q_{i+1} + eq_{i-1}/1+e$ in Eq. (63) one gets:

$$q_{i+1} = -eq_{i-1} + (1+e)P_K \left( \frac{2q_i - (1-e)q_{i-1} + h^2 F_i}{1+e} \right) \quad (64)$$

which is exactly the first equality in Eq. (62).

**Remark 14:**

Notice that the proximations in Eqs. (56) and (58) can also be formulated as projections in the kinetic metric but the set of projection (ie, $V(q_{m,i})$ in (56) and $V(q_{i+1})$ in (58)) depends on the step $i$. Here the set of projection is the same at each step: one always projects on $K$. Another common point between the schemes (55) and (62) is that they do not require the calculation of contact efforts.

It should be noted that the constraints are satisfied at each time step by the average position $\bar{q}_i = q_{i+1} + eq_{i-1}/1+e$. Moreover, if $\bar{q}_i$ belongs to $\text{Int}(K)$, equation (64) reduces to

$$q_{i+1} = 2q_i - q_{i-1} + h^2 F_i \quad (65)$$

which is a classical second-order approximation for the equation of the free motion of the system. Furthermore, the constraints on positions and the impact law are taken into account at the same time by using the average position $\bar{q}_i$. Due to the choice for $F_i$ in Eq. (62), Eq. (64) is an implicit equation as soon as $\bar{Q}$ depends on $p = M(q)\dot{q}$. In such a case some iteration (Newton-Raphson like algorithm) has to be used at each step to compute $q_{i+1}$ (see eg [109]). Notice that the formulation (57) and (58) used in the convergence proof of the discretized Moreau's sweeping process is also an implicit formulation. For ODEs or DAEs implicit methods are known to possess larger domains of conditional stability ([218], page 239). What about MDIs?

As pointed out above the scheme in Eq. (62) does not require the systematic detection of impact times and does not need to refine the time step when the discrete positions are close to the boundary of $K$. As an example let us consider the 1D model problem described by $K = \mathbb{R}^+$, $\bar{Q}(q,p,t) = 0$ and the initial data $q(0) = 1$, $\dot{q}(0) = -1$. The motion is given by

$$q(t) = 1 - t \quad \text{if } t \in [0,1], \quad q(t) = e(t-1) \quad \text{if } t \ge 1. \quad (66)$$

The algorithm (62) yields

$$q_{i+1} = -eq_{i-1} + \max[2q_i - (1-e)q_{i-1}, 0].$$

Let us choose $q_0 = 1$ and $q_1 = 1 - h$. For $i \ge 1$, let $w_i = 2q_i - (1-e)q_{i-1}$. For $i = 1$, we have

$$w_1 = 2q_1 - (1 - e)q_0 = (1 + e) - 2h$$

and $w_i > 0$ if and only if $h < (1+e)/2$. From now on, let us assume that $h < (1+e)/2$. We define $n = \inf\{k \in \mathbb{N}^*: w_k < 0\}$. We have $n \ge 2$ and for all $i \in \{2, \dots, n\}$ we get

$$q_i - q_{i-1} = q_{i-1} - q_{i-2} = -h$$

<figure><img src="image_4.png" /><figcaption>Fig. 3 Elementary block</figcaption></figure>
hence $q_i=1-ih$.

One can also observe that $q_{n-1}>0$ since $w_{n-1}=(1+e)q_{n-1}-(1-e)h\ge 0$. Since $w_n<0$, we have $q_{n+1}=-eq_{n-1}<0$. Let us compute $w_{n+1}$:

$$ 
\begin{aligned}
w_{n+1} &= 2q_{n+1} - (1-e)q_n = -2eq_{n-1} \\
&\quad - (1-e)(-eq_{n-2} + w_{n-1}) \\
&= -e(2q_{n-1} - (1-e)q_{n-2}) - (1-e)w_{n-1} \\
&= -w_{n-1} \le 0.
\end{aligned} 
\quad (67) $$

It follows that $q_{n+2}=-eq_n$ and

$$ w_{n+2} = 2q_{n+2} - (1-e)q_{n+1} = -2w_n \ge 0 \quad (68) $$

$$ q_{n+3} = -eq_{n+1} + w_{n+2} = e^2 q_{n-1} + w_{n+2} \ge 0. $$

Moreover, a straightforward induction gives $q_i=q_{n+2}+eh(i-n-2)$ for all $i \ge n+2$.

The approximate positions $q_i$ do not satisfy the constraints at each time step, similarly to the discretized sweeping process. Of course, the average positions $\bar{q}_i=q_{i+1}+eq_{i-1}/1+e$ belong to $K$ for all $i \ge 1$ but not necessarily $q_{i+1}$ and $q_{i-1}$. In this 1D example, at least one and at most two approximate positions are outside $K$.

It can be proved that the penetration tends to 0 uniformly with respect to the time step since

$$ \|q_i - \bar{q}_i\| \le h \|\dot{q}_i\|_{L^2} = O(h). $$

Moreover, modification has been proposed in [155] that assures penetration on one step only. It follows that the scheme is at most of order 1. This fact is not surprising since we approximate a second-order differential equation in terms of positions only. Such a choice means that we prefer to propose a fast scheme than a very accurate one. We should remind that we have to deal with highly nonlinear problems that are often very sensitive to initial data. Thus the accuracy of the scheme may be less important. But its convergence is always a crucial property.

### 7.3.2 Convergence

The convergence of this scheme is proved. The result has been established first in the case of convex constraints with a trivial mass matrix in [54] and [216], then an extension to non convex constraints (but still trivial mass matrix) is proposed in [217], finally the general case is studied in [55] and [56]. The proof follows the sketch described in Section 7.5 and is based on a rather natural geometrical idea: with an appropriate choice of local coordinates we can describe the set of admissible positions $K$ by a half-space and obtain a simpler expression of the constraints and the projection on $K$. Nevertheless, the change of coordinates introduces other quadratic terms in the algorithm due to curvature effects. These new terms interact with the constraints and create serious difficulties in the study of the scheme.

### 7.3.3 Implementation

The scheme could seem to be difficult to implement since it requires to solve at each time step an implicit equation involving a projection on $K$ for a Riemannian metric. But, in

most of the cases, it is possible to describe the system with a choice of generalized coordinates such that the projection term is easy to compute exactly or to approximate. This scheme has been implemented on different examples. In [54], [219], and [220], a 1D model of tight joints is considered; in [54] and [220], the discretization of a guided beam is examined; and in [109], the motion of a slender bar is studied. For the two first cases the mass matrix of the system is constant, but not for the last one since a model in large deformations is considered. The results are compared to other numerical results obtained by systematic detection of impact times for the first and last examples and to experimental results for the last two examples. For instance in [109] free-motion dynamics is discretized with a Newmark scheme (implicit algorithm) with $\gamma=0.5$ and $\beta=0.25$ (such values guarantee A-stability [26], Ch 7.1.5), and impacts are detected by a root-finding algorithm (a Newton scheme initialized by the values computed at the preceding step). In all the cases, the performances of the scheme in Eq. (62) were satisfactory: the scheme is substantially faster than the detection method (almost 40% faster in [109], three times faster in [219]) and gives good approximation of the motion even with rather large time steps and/or on long time intervals (see [219]). This scheme has also been implemented in [151-153]. Therein they consider a spring-dashpot system with one degree-of-freedom. They compute its motion for two sets of data leading to periodic motions, by using the time-stepping scheme (62) and different event-driven schemes. They compare the numerical results to the analytical solution. In the two cases the scheme (62) is of order 1 and is faster than the event-driven schemes for a given time-step. On the other hand, due to its low order, the scheme (62) requires a larger computing time than the event-driven schemes, in order to approximate the solution to a given precision.

## 7.4 Velocity-impulse formulations

These works have been performed by a group composed of Stewart [33], Trinkle [84], Pang [85], Anitescu [86], and Potra [87]: different formulations have been proposed by Stewart et al following the works of Moreau and Lötstedt. They may be seen as variants of the semi-implicit Euler method for DAEs (which are very attractive for systems as in Eq. (6) [25], page 524). Some details have been given concerning friction formulation, see Eqs. (24) and (25). In Eq. (70), a polyhedral approximation of the friction cone is used, so that the conditions in Eq. (25) are modified. The algorithms have the general form:

$$ 
\left\{ 
\begin{array}{l}
M(q_{i+1})(\dot{q}_{i+1} - \dot{q}_i) = \\
\qquad \nabla \phi(q_i) N_{i+1} + D(q_i) \beta_{i+1} + hQ(q_i, \dot{q}_i) \\
q_{i+1} = q_i + h\dot{q}_{i+1}
\end{array}
\right. 
\quad (69) $$
$$
\left\{
\begin{array}{l}
\mathbf{N}_{i+1} \ge 0, \quad \nabla \phi(q_i)^T \dot{q}_{i+1} \ge 0, \quad \mathbf{N}_{i+1} \nabla \phi(q_i)^T \dot{q}_{i+1} = 0 \\
\beta_{i+1} \ge 0, \quad \lambda_{i+1} \mathbf{e} + D(q_i)^T \dot{q}_{i+1} \ge 0, \\
\beta_{i+1}^T [\lambda_{i+1} \mathbf{e} + D(q_i)^T \dot{q}_{i+1}] = 0 \\
\lambda_{i+1} \ge 0, \quad f \mathbf{N}_{i+1} - \mathbf{e}^T \beta_{i+1} \ge 0, \\
\lambda_{i+1} [f \mathbf{N}_{i+1} - \mathbf{e}^T \beta_{i+1}] = 0.
\end{array}
\right.
\qquad (70)
$$

It is clear that the complementarity relations in Eq. (70) apply when the constraint is active only (including impact times), ie, if $\phi(q_i + h\dot{q}_i) < 0$. If $\mathbf{N}_{i+1} = 0$, then $\beta_{i+1} = 0$ as well. In Eqs. (69) and (70), $\mathbf{e}^T = [1, 1, \dots, 1] \in \mathbb{R}^e$ where $e$ is the number of edges of the polyhedral approximation of the friction cone, and $f$ denotes here the friction coefficient. Hence $\beta \in \mathbb{R}^e$ as well. Indeed the friction cone is approximated by the polyhedral set $\hat{C}(q) = \{\mathbf{N}\nabla\phi + D(q)\beta, \mathbf{N} \ge 0, \beta \ge 0, \mathbf{e}^T\beta \le f\mathbf{N}\}$. The columns $D^l$ of the matrix $D(q)$ are vectors that span the tangent subspace at the contact point. It is also assumed that there is always $l$ and $j$ such that $D^l = -D^j$. For instance in the case of a planar point-mass system $D^1 = \mathbf{t}$ and $D^2 = -\mathbf{t}$, and the tangential part of the contact reaction is given by $D(q)\beta = D^1\beta_1 + D^2\beta_2$, where $\mathbf{t}$ is the tangent direction at the contact point. Let $q^T = (x, y)$ with $x$ the tangential coordinate. The two sets of complementary slackness conditions in Eq. (70) then become (in continuous time, so drop the indices):

$$
\begin{align}
0 \le (\lambda + \dot{x}) \perp \beta_1 \ge 0 \nonumber \\
0 \le (\lambda - \dot{x}) \perp \beta_2 \ge 0 \tag{71} \\
0 \le (f \mathbf{N} - \beta_1 - \beta_2) \perp \lambda \ge 0. \nonumber
\end{align}
$$

Assume for instance that there is a sliding motion with $\dot{x} > 0$. Then since $\lambda \ge 0$, one has $\beta_1 = 0$. Now necessarily $\lambda > 0$ since $\lambda \ge \dot{x} > 0$. Thus $\beta_2 = f\mathbf{N}$. If $\dot{x} < 0$ one would find $\beta_1 = f\mathbf{N}$ and $\beta_2 = 0$. Consider now $\dot{x} = 0$. One finds that if $\lambda \ne 0$, then $f\mathbf{N} = \beta_1 + \beta_2$. Since $\lambda > 0$ the third relation implies $\beta_1 + \beta_2 \le f\mathbf{N}$: the contact reaction is inside the friction cone. Such a reasoning generalizes to 3D cases and the pyramid cone formulation. The conditions in Eq. (25), when discretized, yield a highly nonlinear complementarity problem [91].

It is assumed that $Q(q, \dot{q})$ derives from a potential energy. It can be shown that if $M(q) \equiv M$ then the set of discretized equations in (70) can be transformed in an LCP whose unknown is the vector $(\mathbf{N}_i, \lambda_i, \beta_i)$, and this LCP possesses a solution. Thus the algorithm can be advanced in time. Notice that $\lambda_i$ and $\beta_i$ are to be considered as impulses since they are proportional to forces times $h$. The last two sets of complementarity conditions in (70) represent an approximation of Coulomb model, where the friction cone is replaced by a polyhedral set (a pyramid if $e=4$). They are the Kuhn-Tucker conditions for the maximum dissipation principle in Eq. (72). This is what allows one to express this model as a LCP, using tools from convex analysis, as pointed out in Remark 5. Otherwise in dimension 3, one would end up with a NCP to formulate dry friction. The tangential impulse is represented by $\beta$, while $\lambda$ is a slack variable.

In relationship with Section 3.4.2, let us note that the complementarity conditions in Eq. (70) are equivalent to ($\Lambda_q$ is the generalized impulse in coordinates $q$)

$$
\Lambda_{q,i+1} = \arg \max_{z \in \hat{C}(q_i)} \{-z^T \dot{q}_{i+1}\}. \qquad (72)
$$

If they were formulated with $\nabla\phi(q_{i+1})$ and $D(q_{i+1})$ then one would get $\Lambda_{q,i+1} = \arg\max_{z\in\hat{C}(q_{i+1})}\{-z^T\dot{q}_{i+1}\}$, yielding however a NCP instead of a LCP to be solved at each step.

**Remark 15:** The polyhedral approximation of $\mathcal{C}$ hampers to have the tangential velocity in a direction exact opposite to the tangential force, since the latter belongs to one of the corners of the polygon [89]. If the number of faces of the polyhedral set is increased, this effect is decreased, but the number of LCPs to be solved at each step is increased as well. Hence the interest of looking for solutions that do not approximate the cone $\mathcal{C}$ and minimize the number of equations to be used [83].

In case of multiple $m$ contacts, the (generalized) friction cone (see eg, [18,28,201]) is taken as the sum of the friction cones at each contact, ie, $\mathcal{C} = \sum_{i}^{m} C_i$ (in other words an element $R \in \mathcal{C}$ can be written as $R = \sum_{i}^{m} R_i$ with each $R_i \in C_i$). The vectors $\nabla\phi$, $\lambda$, $\beta$, $\mathbf{N}$, are then simply constructed as the concatenation of the vectors for each contact. In the multi-contact case (say 2), it may therefore occur that $\dot{q}_{i+1}$ and $\Lambda_{q,i+1}$ are computed from Eq. (72) inside $\hat{C}_1 + \hat{C}_2$, whereas $\nabla\phi^{(1)T}(q_i)\dot{q}_{i+1} > 0$, ie, there is detachment from constraint 1 at $t_{i+1}$. However as we pointed out above this implies in turn that $\mathbf{N}_{i+1}^{(1)} = 0$ and $\beta_{i+1}^{(1)} = 0$ too. Thus the impulse $\Lambda_{q,i+1}$ is calculated inside $\hat{C}_2(q_i)$. It is not clear at this stage whether this formulation would allow to get satisfaction of Coulomb law as $h \to 0$ or not. It seems that the main obstacle towards such a result is the possibility of Painlevé-like effects, which have not yet been understood for several contacts. Another problem that might occur as $h \to 0$ is that there may be some cyclic situations in which $\hat{C}$ switches infinitely fast between several values, eg, $\hat{C}_1 + \hat{C}_2$, $\hat{C}_1$ and $\hat{C}_2$. What happens in the limit then? This problem is—at least at the philosophical level—of the same nature as that of multiple impacts. Interestingly enough, the conclusions in [125] about the relationship between the limit of a penalized problem with Coulomb friction and the LCP (as constructed in Section 3.3.3), hold also only for the one-contact sliding mode. Singular perturbation analysis is used and the stability of the so-called boundary layer is no longer equivalent to the existence and uniqueness of solution to the LCP($\lambda$) when $m_s \ge 2$.

**Impact rules:** The impact rules are chosen according to Moreau's *maximal dissipation principle* [104] in [33] and [34], ie, Eqs. (69) and (70) represent a plastic generalized impact as in Eq. (12). In [84,86], the collision rule of Pfeiffer and Glocker [7] is chosen, ie, one solves a LCP to compute the velocity at the end of the compression phase, and another LCP for the velocity at the end of the expansion phase (recall that such an approach relies on Darboux-Keller's model for collisions [18]). Poisson coefficients of restitution are chosen. Energetical consistency holds when all coefficients are chosen equal to a unique value $e_p$. Roughly one considers