Alexandre Patrick Ballard
email: patrick.ballard@dalembert.upmc.fr

Alexandre Charles
email: alexandre.charles@safrangroup.com

An overview of the formulation, existence and uniqueness issues for the initial value problem raised by the dynamics of discrete systems with unilateral contact and dry friction

Keywords: dynamics, rigid solid, discrete system, impact, friction, contact, existence, uniqueness' dynamique, solide rigide, système discret, impact, frottement, contact, existence, unicité'

An overview of the formulation, existence and uniqueness issues for the initial value problem raised by the dynamics of discrete systems with unilateral contact and dry friction

Introduction

The dynamics of rigid solids with contact and dry friction conditions is a venerable subject, which was developed mainly in the second half of the nineteenth century and the beginning of the twentieth century to answer some questions raised by engineering. Then, the attention of engineers began to be driven towards elasticity and continuum mechanics and less attention was paid to frictional contact multibody dynamics. In the seventies, a renewal of interest occurred, mainly driven by the development of numerical modelling in granular dynamics and the control issues associated with robotics. For both concerns, it turned out that the foundations of the venerable theory were not firm enough and that they should be reconsidered in the light of more recent mathematics.

A fundamental impulse was given simultaneously by Michelle Schatzman [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF] and Jean Jacques Moreau [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF] who first considered an evolution problem for the configuration q : [0, T ] → R d in the framework of functions whose second derivative (in the distributional sense) is a Radon measure. At that time, the antique point of view of different systems of equations applying to the different phases of motion (without or with active contact), was still prevailing. Their new point of view permitted to formulate for the first time, a mathematical evolution problem associated with multibody contact dynamics. It also paved the road for the design of efficient strategies for numerical computing and enabled the first investigations on the stability and control issues that are crucial in the analysis of the problems facing frictional contact events in robotics.

The seminal work of Michelle Schatzman and Jean Jacques Moreau also initiated a series of contributions on the general formulation and the mathematical analysis of the initial value problem associated with multibody contact dynamics. A brief sketch of the history follows.

-The first studies about the formulation were restricted to the model problem of the dynamics of a particle evolving in an admissible region of R d bounded by an obstacle. In [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF], Michelle Schatzman formulated a consistent evolution problem in the frictionless case, under the additional restriction that the admissible region is convex. She was able to successfully implement a penalty method to prove the existence of a solution for the initial value problem. Her original work was restricted to impacts preserving the kinetic energy (the so-called elastic impact law) and an external force depending only on time. This result was generalized later by Paoli in her PhD thesis [START_REF] Paoli | Mouvement à un nombre fini de degrés de liberté avec contraintes unilatérales : cas avec perte d'énergie[END_REF] to the case of an arbitrary impact law and an external force possibly also depending on current position and velocity. In parallel, an alternative strategy for proving the existence of solution was designed by Monteiro Marques [START_REF] Monteiro Marques | Differential inclusions in nonsmooth mechanical problems -shocks and dry friction[END_REF]. He introduced a time-stepping approximation and proved the convergence (of a subsequence) towards a solution. It was restricted to the completely inelastic impact law (zero restitution coefficient), but he was able to relax the convexity assumption of Schatzman on the admissible region. More importantly, he was able to generalize to the case where the contact with the obstacle obeys the Coulomb law of dry friction with given friction coefficient µ (the frictionless case is recovered by taking µ = 0). One benefit of this new strategy is that it directly suggested an algorithm for numerical computations. The time-stepping approach was further developed by Paoli [START_REF] Paoli | Time-stepping approximation of rigid-body dynamics with perfect unilateral constraints. I and II[END_REF] and [START_REF] Paoli | A proximal-like algorithm for vibro-impact problems with a non smooth set of constraints[END_REF], who, in particular, extended it to the case of an arbitrary restitution coefficient. Her work, however, is up to now concerned only with the frictionless case.

-It was recognized very early by Michelle Schatzman that issues should be expected with the uniqueness of the solution for the initial value evolution problem. In particular, she exhibited in [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF] a striking example of multiple solutions for the unilateral dynamics of the one degree-of-freedom particle submitted to an external force that is a C ∞ function of time. This issue was further considered by Percivale [START_REF] Percivale | Uniqueness in the elastic bounce problem, i[END_REF][START_REF] Percivale | Uniqueness in the elastic bounce problem, ii[END_REF] who noticed that uniqueness of the solution could be recovered in the one degreeof-freedom problem considered by Schatzman, provided that the given external force was assumed to be not only a C ∞ but an analytic function of time. His work was suggesting that uniqueness could be expected in general, provided a regularity assumption of analyticity on the data. This was proved in full generality for the frictionless problem by Ballard in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF] and [START_REF] Ballard | Dynamics of rigid bodies systems with unilateral or frictional constraints[END_REF]. Local uniqueness in the analytic framework was also exploited in this work to design a third alternative strategy (in addition to penalty and time-stepping methods) to prove the existence of the solution. This new strategy turned out to yield more general (except for the additional assumption of analyticity) existence results than those which were available at that time from the penalty and time-stepping methods, encompassing the multi-constraint case with an arbitrary impact law. This strategy of proof was also adapted in [START_REF] Ballard | Existence and uniqueness for dynamical unilateral contact with coulomb friction: a model problem[END_REF] and [START_REF] Charles | Existence and uniqueness of solutions to dynamical unilateral contact problems with Coulomb friction: the case of a collection of points[END_REF] to the dynamics of a point particle with contact conditions and Coulomb friction. It yields a slightly more general (except for the analyticity assumption) existence result than that of Monteiro Marques, and provides in addition, the uniqueness of the solution from a given initial value.

-In the eighties and the nineties, most of the articles that appeared on the subject of the mathematical formulation of the initial value problem and the issues of existence and uniqueness of solutions were restricted to the model problem of a point particle evolving in R d , or rather in an admissible region of R d . There is one noteworthy exception: the seminal article [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF] by Jean Jacques Moreau.

In this article, Jean Jacques Moreau addressed the formulation of the dynamics of a collection of rigid bodies submitted to frictionless unilateral constraints such as the ones arising from the noninterpenetration conditions. The framework is from the beginning that of Lagrange in which the motion is represented as a curve q(t) in the configuration space, identified with a subset of R d . In that framework and in accordance with the ideas of Lagrange, the reaction force that appear in the formulation is a generalized reaction force. In particular, the detailed distribution of reaction forces in the real world (meaning forces in R 3 from one body onto another) is generally undefined. The existence and uniqueness result of Ballard in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF] and [START_REF] Ballard | Dynamics of rigid bodies systems with unilateral or frictional constraints[END_REF] applies to this general framework for frictionless unilateral multibody dynamics, under the assumption of analyticity of the data. It yields existence and uniqueness for the motion q(t) and for the generalized reaction force, but the detailed distribution of reaction forces in the real world remains undetermined in general, and there may be several such distributions that are compatible with the generalized reaction force that is associated with the solution. The use of generalized forces, originating in Lagrange's idea, has now been made systematic all over continuum mechanics under the name 'Principle of Virtual Power'. It conveys the idea that the representation of forces within a mechanical theory must be made consistent with the kinematics: forces must be taken in the dual space of the velocity space, the duality product between forces and velocities being nothing but the power. Surprisingly, coming to frictional unilateral multibody dynamics, the use (originating from the nineteenth century) was to invoke the Coulomb law of dry friction applying to reaction forces in the real world, contradicting Lagrange's point of view and the Principle of Virtual Power that stipulate that the dynamics should be formulated in terms of appropriate generalized forces. However, in some simple cases involving a small number of contacting rigid bodies, a mathematical evolution problem can still be formulated based on this historical point of view [START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of painlevé's problem[END_REF], although no modern formulation of the general evolution problem has ever been obtained in that setting. Unsurprisingly, adopting the historical point of view of real world forces (when it is possible) raises some inconsistencies, known as Painlevé and Kane paradoxes [START_REF] Painlevé | Sur les lois du frottement de glissement[END_REF] and [START_REF] Kane | A dynamics puzzle[END_REF][START_REF] Kane | Dynamics : theory and applications[END_REF]. Because of this lack of consistency within the usual view about the formulation of frictional unilateral multibody dynamics, no precise formulation of the evolution problem was obtained and the only mathematical investigations about existence and uniqueness of solutions for frictional unilateral dynamics have been so far restricted to the case of a finite collection of point particles. That was only recently that such a consistent general formulation in the line of Lagrange's ideas about the use of generalized forces was derived by Charles [START_REF] Charles | The formulation of dynamical contact problems with friction in the case of systems of rigid bodies and general discrete mechanical systems-Painlevé and Kane paradoxes revisited[END_REF] in his PhD thesis.

The one-degree-of-freedom problem

Consider a point particle, of unit mass, that is constrained to move along a line, the location of which is represented by the abscissa q ∈ R. We assume that an obstacle is located at the origin so that the particle is constrained to remain in the half-line defined by q ≥ 0. To enforce this constraint during an arbitrary motion q(t) of the particle, an unknown reaction r force must be added in the equation of motion:

q = f + r, (1)
where the external force f (t) is supposed to be a given (integrable) function of time only (for the sake of simplicity). The usual physical assumption is that the existence of an obstacle has no influence on the motion of the particle when contact is not active, and the reaction force must therefore be supported in those instants where the contact is active:

Supp r ⊂ t q(t) = 0 .
Elementary examples then shows that velocity jumps cannot be avoided in this framework and the acceleration that appears in the equation of motion (1) should be understood in the sense of (Schwartz's) distributions. Hence, the reaction force r should not be expected to be a function but rather a distribution. As it is usually assumed that the obstacle is only able to repel the particle, a nonnegativity assumption also has to be required on the unknown reaction force: r ≥ 0. In the extended framework where r is a distribution, this can only mean that the distribution r returns nonnegative real values when tested by means of a nonnegative C ∞ trial function with compact support. But, it is a classical (and easy) result that such a nonnegative distribution must actually be a nonnegative measure. Hence, given a bounded time interval [0, T ], the largest possible functional space in which the motion q(t) can be sought is the space MMA([0, T ]) (the acronym standing for 'Motions with Measure Acceleration') of those distributions on [0, T ] whose second derivative is a Radon measure (q ∈ M ([0, T ]). Distributions in that space MMA([0, T ]) are actually continuous functions, admitting left and right derivatives q-(t), q+ (t) (in the classical sense) at every instant t ∈ ]0, T [. The side derivatives q-(t) and q+ (t) are actually equal, except possibly at some instants belonging to a countable subset of [0, T ]: the impact instants. The two functions q-(t) and q+ (t) are functions with bounded variation. Finally, given an initial condition (q 0 , v 0 ) ∈ R + × R compatible with the obstacle (q 0 = 0 ⇒ v 0 ≥ 0), the evolution problem associated with the unilateral dynamics of the particle reads as follows.

Problem P 1 . Find q ∈ MMA([0, T ]) and r ∈ M ([0, T ]) such that:

q(0) = q 0 , q+ (0

) = v 0 , -q = f + r, -∀t ∈ [0, T ], q(t) ≥ 0, -Supp r ⊂ t ∈ [0, T ] q(t) = 0 , -r ≥ 0, -q(t) = 0 ⇒ q+ (t) = -e q-(t).
Here, the last line is an additional requirement with respect to the introductory discussion. If it were not stated, a particle free of external force impacting the obstacle could either subsequently remain stuck on the obstacle, or bouncing according to a sign-reversed velocity (infinitely many intermediate choices are possible), and all of these events would be compatible with both the equation of motion and the contact conditions. Such an indeterminacy was already noted by Newton who introduced the concept of a restitution coefficient to recover determinism. The indeterminacy originates in the fact that real bodies are always deformable and their bouncing is actually governed by the deformation waves that travel within the body during an impact. As the model of a point particle is too coarse to describe these deformation waves, the outcome of the impact remains undetermined and should be artificially added to the equations by means of an impact law, which has the same status as constitutive laws in continuum mechanics. The restitution coefficient e is supposed to be given in the interval [0, 1] (the choice of e = 1 is usually baptized the elastic impact law and that of e = 0, the completely inelastic impact law ).

Actually, it was only recently that it was discovered that the original proposition of Newton to add the impact law fails to yield uniqueness of the solution, even in the case where f is supposed to be a C ∞ function of time. The first discovery of this striking fact seems to be that of Bressan [START_REF] Bressan | Incompatibilità dei Teoremi di Esistenza e di Unicità del Moto per un Tipo molto Comune e Regolare di Sistemi Meccanici[END_REF], although he did not have such a precise statement of the evolution problem. The construction of such a C ∞ function f enabling multiple solutions q can be found in [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF] in the case of the elastic impact law e = 1 and in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF] in the case of the completely inelastic impact law e = 0. In both cases, the C ∞ functions f that are exhibited are highly oscillating (having in particular infinitely many zeroes in the bounded interval [0, T ]).

The following results about problem P 1 have been proved in various articles.

(i) Suppose that f ∈ L 1 (0, T ) and e ∈ ]0, 1]. For all ε > 0, the penalized initial value problem:

q ε (0) = q 0 , qε (0) = v 0 , qε (t) = f (t) - 2| log e| ε[π 2 + (log e) 2 ] qε (t) sgn -q ε (t) - q ε (t) ε sgn -q ε (t) , for a.a. t ∈ [0, T ],
(where sgn -is the function taking the value 1 on ]-∞, 0[ and 0 on [0, ∞[), has a unique solution q ε ∈ W 2,1 (0, T ). As ε goes to zero, one can extract a subsequence in q ε that converges strongly in W 1,1 (0, T ) towards some q ∈ MMA([0, T ]) that solves problem P 1 . This result was proved in [START_REF] Paoli | Mouvement à un nombre fini de degrés de liberté avec contraintes unilatérales : cas avec perte d'énergie[END_REF], whereas the simpler case e = 1 (where the damping term in the penalty differential equation reduces to 0) had been previously treated in [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF]. Hence, although the penalty approach fits naturally with the particular case of the elastic impact law e = 1, it can be extended to dissipative impact laws.

(ii) Suppose that f ∈ L 1 (0, T ) and e ∈ [0, 1]. Picking n ∈ N \ {0} and setting h = T /n, we define a sequence of approximants q n ∈ C 0 ([0, T ]) by the following induction.

-

Q 0 n = q 0 , V 0 n = v 0 , -∀i ∈ {1, 2, . . . , n}, F i n = 1 h ih (i-1)h f, V i n = V i-1 n + hF i n , if Q i-1 n > 0 -eV i-1 n + (1 + e)V i-1 n + hF i n + , if Q i-1 n ≤ 0, Q i n = Q i-1 n + hV i-1 n -q n (t) = Q i-1 n + t -(i -1)h V i-1 n
, ∀t ∈ (i -1)h, ih , where x + = max{x, 0} stands for the positive part function. Then, it was proved in [START_REF] Paoli | A proximal-like algorithm for vibro-impact problems with a non smooth set of constraints[END_REF] that a subsequence of q n (t) converges strongly in W 1,1 (0, T ) towards some q ∈ MMA([0, T ]) that solves problem P 1 . The simpler case e = 0 had been previously treated in [START_REF] Monteiro Marques | Differential inclusions in nonsmooth mechanical problems -shocks and dry friction[END_REF]. Hence, although the timestepping approach fits naturally with the particular case of the completely inelastic impact law e = 0, it can be extended to the general case of an arbitrary restitution coefficient e ∈ [0, 1].

(iii) Uniqueness for problem P 1 can be recovered in MMA([0, T ]), if the function f (t) is more regular than C ∞ . In the analytic case, uniqueness was first proved for the particular case e = 1 in [START_REF] Percivale | Uniqueness in the elastic bounce problem, i[END_REF], and for an arbitrary restitution coefficient e ∈ [0, 1] in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF]. Needless to say, these results apply in the case of piecewise analyticity.

(iv) If the function f (t) is analytic, or piecewise analytic, then the solution mapping:

   R 2 → C 0 ([0, T ]) (q 0 , v 0 ) → q(t)
is continuous. In other words, the unique solution of problem P 1 depends continuously on the initial data. A proof is to be found in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF].

Frictional unilateral dynamics of a point particle

We now consider the motion of a point-particle, of unit mass, in a region of R d defined by the unilateral constraint ϕ(q) ≥ 0, where ϕ is a smooth (of class C 1 at least) function such that:

ϕ(q) = 0 =⇒ dϕ q = 0,
where dϕ q denotes the differential of the real-valued function ϕ at point q. The equation ϕ(q) = 0 defines the geometry of the obstacle, and dϕ q defines the direction of the outward normal to the obstacle. When the particle is in contact with the obstacle, the reaction force exerted by the obstacle will be denoted by r ∈ R d . It can be classically split into normal and tangential parts:

r = r t + r n dϕ q |dϕ q | , with r n = r, dϕ q |dϕ q | ,
where q denotes the location of the particle at contact and •, • , | • | are the canonical scalar product and norm in R d . The case where dry friction between the particle and the obstacle can occur will be considered. The simplest law describing dry friction is the empirical law of Coulomb, which reads as:

r t ≤ µ r n and if r t < µ r n , then q+ t = 0, if r t = µ r n , then ∃λ ∈ R + , r t = -λ q+ t .
(

) 2
where q+ t denotes the tangential (with respect to the obstacle) component of the right-velocity, and µ ≥ 0 is a given friction coefficient. The Coulomb law above can be compactly and equivalently rewritten under the weak form:

∀v ∈ R d , r t , v t -qt + µ r n |v t | -| qt | ≥ 0,
where •, • stands for the canonical scalar (duality) product in R d , as seems to have been first pointed out independently in the sixties by Jean Jacques Moreau and Georges Duvaut. This weak form turns out to be the appropriate form to extend the classical Coulomb law to the case where the reaction force r is a Radon measure (with respect to time):

∀v ∈ C 0 ([0, T ]; R d ), [0,T ] r t , v t (t) -qt (t) + µ r n |v t (t)| -| qt (t)| ≥ 0. Note that if q ∈ MMA([0, T ]; R d ), then q+ : [0, T [ → R d is
a function with bounded variation which is therefore universally integrable (integrable with respect to any measure). Hence, the integral in the above weak form of the Coulomb law is well-defined for q ∈ MMA([0, T ]; R d ) and r ∈ M ([0, T ]; R d ). This weak form can be assumed globally in [0, T ], so that the Coulomb law will be enforced both during the impacts and the possible smooth phases of motion along the obstacle.

Given an initial condition (q 0 , v 0 ) ∈ R d × R d compatible with the obstacle:

ϕ(q 0 ) ≥ 0 and ϕ(q 0 ) = 0 ⇒ dϕ q0 , v 0 ≥ 0,
we can now formulate the initial value problem that governs the dynamics of the particle in frictional contact with the obstacle. The external force f (t; q, q-) is now allowed to depend on the current location and velocity (here, we can equivalently use either q-or q+ , with no influence on the solution).

Problem P 2 . Find q ∈ MMA([0, T ]; R d ) and r ∈ M ([0, T ]; R d ) such that: -q(0) = q 0 , q+ (0) = v 0 , -q = f (t; q, q-) + r, -∀t ∈ [0, T ], ϕ q(t) ≥ 0, -Supp r ⊂ t ∈ [0, T ] ϕ q(t) = 0 , -r n ≥ 0, -∀v ∈ C 0 ([0, T ]; R d ), [0,T ] r t , v t (t) -qt (t) + µ r n |v t (t)| -| qt (t)| ≥ 0, -ϕ q(t) = 0 ⇒ q+ n (t) = -e q- n (

t).

The following results about problem P 2 have been proved in various articles.

(i) In the frictionless case µ = 0, with given restitution coefficient e ∈ ]0, 1], in the case where f is a continuous function that is globally Lipschitz with respect to (q, q), then it is possible [START_REF] Paoli | Mouvement à un nombre fini de degrés de liberté avec contraintes unilatérales : cas avec perte d'énergie[END_REF] to extend the penalty method described in section 2 (i), at least when the admissible subset of R d defined by ϕ(q) ≥ 0 is convex. In the same way, there exists a subsequence of the penalty approximates q ε that converges strongly in W 1,1 (0, T ), as ε → 0+, towards some q ∈ MMA([0, T ]) that solves problem P 2 (a proof is to be found in [START_REF] Paoli | Mouvement à un nombre fini de degrés de liberté avec contraintes unilatérales : cas avec perte d'énergie[END_REF] in the case of a convex admissible subset of R d ).

(ii) In the frictionless case µ = 0, with given restitution coefficient e ∈ [0, 1], in the case where f is a continuous function that is globally Lipschitz with respect to (q, q), then it is possible [START_REF] Paoli | A proximal-like algorithm for vibro-impact problems with a non smooth set of constraints[END_REF] to extend the time-stepping method described in section 2 (ii). In the same way, there exists a subsequence of the time-stepping approximates q n that converges strongly in W 1,1 (0, T ), as n → +∞, towards some q ∈ MMA([0, T ]) that solves problem P 2 (a proof is to be found in [START_REF] Paoli | A proximal-like algorithm for vibro-impact problems with a non smooth set of constraints[END_REF]).

(iii) In the frictional case µ ≥ 0, and the completely inelastic impact law e = 0, in the case where f is a continuous bounded function of t and q only, Monteiro Marques [START_REF] Monteiro Marques | Differential inclusions in nonsmooth mechanical problems -shocks and dry friction[END_REF] defines a sequence of approximants q n ∈ C 0 ([0, T ]) (n ∈ N \ {0}) by the following induction.

-

Q 0 n = q 0 , V 0 n = v 0 , -∀i ∈ {1, 2, . . . , n}, Q i n = Q i-1 n + hV i-1 n , V i n = V i-1 n + hf (ih, Q i n ), V i n = V i n , if V i n ∈ V (Q i n ) Proj 0; V i n + C (Q i n ) ∩ T (Q i n ) , if V i n / ∈ V (Q i n ), -q n (t) = Q i-1 n + t -(i -1)h V i-1 n , ∀t ∈ (i -1)h, ih , where h = T /n, V (q)
is the set of admissible right-velocities at location q:

V (q) = v ∈ R d dϕ q , v ≥ 0 , if ϕ(q) ≤ 0, R d if ϕ(q) > 0,
(dϕ q is assumed to be nowhere 0), T (q) is the tangent hyperplane:

T (q) = v ∈ R d dϕ q , v = 0 ,
and C (q) is the friction cone:

C (q) = r ∈ R d r, dϕ q ≥ |r| dϕ q / 1 + µ 2 .
Then, there exists a subsequence in the sequence (q n ) that converges strongly in in W 1,1 (0, T ), as n → +∞, towards some q ∈ MMA([0, T ]) that solves problem P 2 (a proof is to be found in [START_REF] Monteiro Marques | Differential inclusions in nonsmooth mechanical problems -shocks and dry friction[END_REF]).

(iv) In the frictional case µ ≥ 0, with given restitution coefficient e ∈ [0, 1], in the case where ϕ is an analytic function and f is an analytic function that is globally Lipschitz with respect to (q, q), then there exists a unique solution to problem P 2 (a proof is to be found in [START_REF] Charles | Existence and uniqueness of solutions to dynamical unilateral contact problems with Coulomb friction: the case of a collection of points[END_REF]).

Frictionless unilateral multibody dynamics

In the preceding sections, only the dynamical evolution of a point particle has been considered. However, the applicability of the theory to situations of practical interest requires to extend it to the case of one or several rigid bodies. The framework must therefore be extended to that of Lagrange about discrete mechanical systems (that encompasses in particular any finite collection of rigid bodies, some of them being possibly connected by so-called perfect joints).

Virtual power and Lagrange equations

A discrete mechanical system is a mechanical system whose arbitrary configuration in the space can be described by a finite number d of independent real numbers: the generalized coordinates denoted by q 1 , q 2 , . . . , q d . Here, "independent" means that it is always possible to conceive a motion for which all generalized coordinates but an arbitrary one, remain fixed. The number d is called the number of degrees of freedom of the discrete mechanical system. The notation q = (q 1 , q 2 , . . . , q d ) will be used and q will be called the abstract configuration (or sometimes simply, the configuration) of the system. A motion of the system is simply a mapping q(t) defined on some time interval and taking abstract configuration values. Its derivative with respect to time is denoted by q(t) and the vector q(t) is called the generalized velocity at time t. The generalized velocity is a convenient mathematical representation of the whole velocity field over the body (or bodies). The kinetic energy K(q, q) is quadratic with respect to the generalized velocity and defines the so-called kinetic matrix M(q) :

K(q, q) = 1 2 q • M(q) • q.
The kinetic matrix is always symmetric and positive. It is actually positive definite, if the system does not involve any massless components, which will always be assumed in the sequel. More precisely, the set Q of all abstract configurations is endowed with the structure of a differentiable manifold. This point of view makes it possible to consider an abstract configuration q with no need of a prior definition of generalized coordinates (that is, a parameterization of the system) and is particularly useful when it is important to distinguish intrinsic quantities, that is, quantities not relying on a specific choice of generalized coordinates on the configuration manifold. In this context, the generalized velocities belong to tangent spaces to the configuration manifold. The kinetic matrix endows each tangent space with a scalar product, so that the configuration manifold is actually a Riemannian manifold. When emphasis is to be put on intrinsic quantities, the following alternative notations will be possibly encountered in the sequel:

q1 , q2 q = q1 • M(q) • q2 , q q = q • M(q) • q. In any case, the point of view of Riemannian manifold on the configuration space, although enlightening, is not strictly needed and there will be no harm in identifying in the sequel the configuration manifold with (an open subset of) R d and the tangent space T q Q at configuration q, with R d endowed with the scalar product M(q). This identification is nothing but considering a particular choice of generalized coordinates.

The generalized force f of Lagrange is defined by means of the virtual power that the internal and external forces develop in any virtual generalized velocity v. This virtual power reads as:

f, v = f i vi ,
(with the usual convention of summation on repeated indices) and shows that the generalized force belongs to the dual space T * q Q of the vector space T q Q of all generalized velocities at configuration q. In particular, the generalized force for a mechanical system with d degrees of freedom has d components. The virtual power point of view is usually used to compute the generalized force from a given distribution of forces f (x) in the real world (here x denotes the space variable in the three-dimensional space). First, the real world velocity ṽ(x) is computed in terms of the generalized velocity v:

ṽ(x) = l(q, x) • v,
where l(q, x) : T q Q → R 3 is a linear mapping depending in general on the current configuration q and x. Then, the virtual power paradigm:

∀v, f, v = f (x) • ṽ(x) dx = t l(q, x) • f (x) dx , v =⇒ f = t l(q, x) • f (x), (3)
provides the expression of the generalized force f in terms of the real world force distribution f (x) (here, in the case where f (x) consists in finitely many point forces, the above integral reduces to a finite sum). Let us point out once more than it is always possible to compute the generalized force from a real world force distribution but that the real world force distribution cannot be recovered in general from the generalized force.

The generalized acceleration γ γ γ of Lagrange is defined by means of the virtual power it develops in any virtual generalized velocity v.

γ γ γ, v = D q Dt , v q = d dt ∂K ∂ qi - ∂K ∂q i vi ,
where D/Dt stands for the covariant derivative along the motion. It satisfies:

dK dt = D q Dt , q q .
The fundamental principle of classical dynamics asserts the Lagrange equation of motion γ γ γ = f, which is equivalent to the principle of virtual power:

∀v ∈ T q Q, D q Dt , v q = f, v .

Frictionless unilateral constraints in Lagrange's setting

In the case of several rigid bodies, it must be expected, in general, that there will be several unilateral constraints: ϕ α (q) ≥ 0 (α = 1, 2, . . . , n), as the non-interpenetration conditions gives rise to at least one such unilateral constraint for each pair of rigid bodies. The set of all the indices for which the corresponding constraint is active in the configuration q will be denoted by J(q): J(q) = α | ϕ α (q) = 0 .

The unilateral constraints ϕ α are supposed to be such that, for any admissible configuration q, the dϕ α,q (α ∈ J(q)) are linearly independent. The realization of the unilateral constraints requires to complement the equation of motion with a (generalized) reaction force r. In the particular case of frictionless unilateral constraints, it takes the form:

r = n α=1 λ α dϕ α,q , with λ α ≥ 0 and Supp λ α ⊂ {t ∈ [0, T ] | ϕ α (q(t)) = 0}.
Although the sole restitution coefficient e is enough to convey the more general frictionless impact law in the case of a single unilateral constraint, a great deal of additional complexity of the impact law is permitted in the case of several active unilateral constraints at an impact. The impact law will be postulated under its more general possible form:

q+ = F q, q-,
where q is the configuration of the system at the impact time, q-the left (impacting) velocity and q+ the right (outgoing) velocity. The function F is a datum of the problem. In a practical situation, it has to be identified by using either experiments or a refined theory such as the mechanics of deformable bodies. Naturally, the function F cannot be arbitrary, it has to fulfil some compatibility conditions in order to be compatible with the equation of motion. More precisely, it must satisfy the three following conditions. (i) ∀q, q-, ∀α ∈ J(q), dϕ α,q , F q, q-≥ 0, (the post-impact velocity must not violate the unilateral constraints), (ii) ∀q, q-, M(q) • F q, q--q-∈ α∈J(q) R + dϕ α,q , (no friction: the generalized reaction force impulse is directed along the normal), (iii) ∀q, q-, F q, qq ≤ qq , (the kinetic energy cannot be increased by an impact).

There exist many functions F satisfying these requirements. One canonical example is the Moreau impact law [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF]. It is based on the decomposition of an arbitrary vector v of the Euclidean vector space T q Q R d , endowed with the scalar product M(q): v = Proj M(q) v; V (q) + Proj M(q) v; N (q) , on the two mutually polar cones:

V (q) = v ∈ T q Q ∀α ∈ J(q), dϕ α,q , v ≥ 0 , N (q) = α∈J(q) R -∇ϕ α,q = α∈J(q) R -M -1 (q) • dϕ α,q ,
where Proj M(q) is the orthogonal projection operator with respect to the scalar product M(q) of T q Q R d . Then, the Moreau impact law [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF] reads as:

F q, q-= Proj M(q) q-; V (q) -e Proj M(q) q-; N (q) (4) in which the restitution coefficient e can be chosen arbitrarily in [0, 1]. It is an easy matter to prove that it fulfils requirements (i), (ii) and (iii) above. Reciprocally, at a location q where only one unilateral constraint is active (card J(q) = 1), requirements (i), (ii) and (iii) above, entail that the most general impact law must take the form (4) of that of Moreau, for some e ∈ [0, 1].

The evolution problem associated with frictionless unilateral multibody dynamics

We are given an initial condition (q 0 , v 0 ) ∈ T Q R d ×R d assumed to be compatible with the unilateral constraints:

∀α, ϕ α (q 0 ) ≥ 0, and ∀α ∈ J(q 0 ), dϕ α,q0 , v 0 ≥ 0, and an impact law function F fulfilling the three requirements (i), (ii), (iii) above. The initial value problem that governs frictionless unilateral multibody now reads as follows.

Problem P 3 . Find q ∈ MMA([0, T ]; Q) and λ α ∈ M ([0, T ]; R) (α ∈ {1, 2, . . . , n}) such that: -q(0) = q 0 , q+ (0) = v 0 , (initial condition) -M(q) • D q Dt = f (t; q, q-) + n α=1 λ α dϕ α,q , (equation of motion) -∀α, ϕ α (q) ≥ 0, λ α ≥ 0, Supp λ α ⊂ t ∈ [0, T ] ϕ α (q(t)) = 0 , ( contact
conditions) -∀t, q+ = F q, q-, (impact constitutive law).

The following results about problem P 3 have been proved in various articles.

(i) If the functions ϕ α , M and f are analytic functions of their arguments, then problem P 3 has a unique maximal solution [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF][START_REF] Ballard | Dynamics of rigid bodies systems with unilateral or frictional constraints[END_REF]. Here, 'maximal solution' means that the solution is not necessarily defined all over [0, T ], but rather on a subinterval that cannot be extended, since a possible blow up at finite time can occur. This blow up at finite time is classically known to be possible in the case where there is no unilateral constraint, but is dismissed whenever f satisfies a global Lipschitz condition with respect to (q, q). The same is true for problem P 3 : if f is, in addition, Lipschitz with respect to (q, q), then the maximal solution is defined all over [0, T ] and problem P 3 has truly one and only one solution [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF][START_REF] Ballard | Dynamics of rigid bodies systems with unilateral or frictional constraints[END_REF].

(ii) The above general existence and uniqueness result of solution for problem P 3 should not put a mask on a major ill-posedness issue that is generally encountered by the solution: it does not depend continuously on the data, in general. This unpleasant feature can be observed in a simple example in R 2 , with f ≡ 0, two affine functions ϕ α (representing two intersecting straight obstacles) and an elastic impact law. On figure 1, two arbitrarily closed initial locations of a point particle are considered, associated with the same initial velocity, and they are seen to produce post-impact velocities that are orthogonal, so that the subsequent motions diverge. This pathology is very unpleasant because it precludes any rational attempt to compute an approximation of the unique solution of problem P 3 , in general. However, this issue is specific to the situation where the unique solution of problem P 3 experiments at least one multiple impact (that is, passing through a location q at which several unilateral constraints are simultaneously active: card J(q) ≥ 2). Hence, the difficulty is connected with multiple impacts. However, all the multiple impacts do not give rise to this pathology. In particular, the counter-example of figure 1 would be impossible to build in the case where the two straight obstacles would intersect with a right angle. This was generalized as follows by Ballard.

Proposition 1 [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF] Consider problem P 3 in which the data ϕ α , M and f are analytic functions of their arguments and where the impact law is chosen to be that of Moreau (4) for an arbitrary e ∈ [0, 1]. Suppose, in addition, that f is Lipschitz with respect to (q, q), which ensures that problem P 3 admits a unique solution for any admissible initial condition. Consider such an initial (q 0 , v 0 ) and the corresponding solution q ∈ C 0 ([0, T ]; R d ). Suppose that this solution q(t) experiments only multiple impacts that are right with respect to the scalar product M -1 (q):

∀t ∈ [0, T ], ∀α, β ∈ J(q(t)), dϕ α,q(t) • M -1 (q(t)) • dϕ β,q(t) = 0, (5)
then the mapping:

   T Q → C 0 ([0, T ]; Q) (q 0 , v 0 ) → q(t)
is continuous at (q 0 , v 0 ).

Physically, condition (5) (the angle between two arbitrary constraints at a multiple impact is a right angle) is satisfied for a multiple impact consisting of simultaneous collisions at different locations by two pairs of solids in a dynamical evolution of a collection of rigid bodies, whereas it is generally not satisfied in the case of three bodies colliding simultaneously. This issue of continuous dependence of the solution of problem P 3 with respect to the data was further investigated by Paoli [START_REF] Paoli | Continuous Dependence on Data for Vibro-impact Problems[END_REF] who proved that, in the particular case where e = 0 (completely inelastic impact law), the orthogonality condition in proposition 1 can be weakened into:

∀t ∈ [0, T ], ∀α, β ∈ J(q(t)), dϕ α,q(t) • M -1 q(t) • dϕ β,q(t) ≤ 0, that is, in the case of the completely inelastic impact law, the angles between the constraints at a multiple impact does not need to be right, but only acute, to ensure continuous dependence of the solution on the initial condition. She also proved that this condition of orthogonality of the constraints (or acute angles, in the case of the completely inelastic impact law) ensures not only continuous dependence of the solution on the initial condition, but also on the other data of the problem, such as M , f and ϕ α , in the sense of appropriate topologies. (iii) Assuming that the impact law in problem P 3 is that of Moreau (4) for an arbitrary e ∈ [0, 1].

Then, Paoli [START_REF] Paoli | A proximal-like algorithm for vibro-impact problems with a non smooth set of constraints[END_REF] defines a sequence of approximants q n ∈ C 0 ([0, T ]) (n ∈ N \ {0}) by the following induction.

-

Q 0 n = q 0 , V 0 n = v 0 , -∀i ∈ {1, 2, . . . , n}, Q i n = Q i-1 n + hV i-1 n , V i n = V i-1 n + h M -1 (Q i n ) • f (ih, Q i n , V i-1 n ), V i n = -e V i-1 n + Proj M(Q i n ) V i n + e V i-1 n ; V (Q i n ) , -q n (t) = Q i-1 n + t -(i -1)h V i-1 n , ∀t ∈ (i -1)h, ih ,
Now, suppose that the data ϕ α , M and f in problem P 3 , are analytic functions of their arguments and that the impact law is chosen to be that of Moreau (4) for an arbitrary e ∈ [0, 1]. Suppose, in addition, that f is Lipschitz with respect to (q, q) and:

∀q, ∀α ∈ J(q) dϕ α,q • M -1 (q) • dϕ β,q = 0, if e ∈ ]0, 1] , ∀q, ∀α ∈ J(q) dϕ α,q • M -1 (q) • dϕ β,q ≤ 0, if e = 0, (6)
then, Paoli proves in [START_REF] Paoli | A proximal-like algorithm for vibro-impact problems with a non smooth set of constraints[END_REF] that the sequence q n converges strongly in W 1,1 (0, T ) towards the unique solution q of problem P 3 . Actually, her proof seems to cover also the situation where the conditions [START_REF] Paoli | A proximal-like algorithm for vibro-impact problems with a non smooth set of constraints[END_REF] are not necessarily fulfilled for all q, but only for all the locations q(t) (t ∈ [0, T ]) associated with the unique solution q of problem P 3 .

An alternate time-stepping scheme is proposed in [START_REF] Paoli | Time-stepping approximation of rigid-body dynamics with perfect unilateral constraints. I and II[END_REF] with proof of the same results of convergence.

Frictional unilateral multibody dynamics

In most use of unilateral multibody dynamics in engineering problems, it is desirable to account for (dry) friction. Unfortunately, frictional unilateral multibody dynamics has not reached yet the same stage of completion as in the frictionless case.

The need for a general formulation

It is very surprising that the first attempt to obtain a general formulation of an evolution problem associated with frictional unilateral multibody dynamics was only that of Charles in his PhD thesis [START_REF] Charles | The formulation of dynamical contact problems with friction in the case of systems of rigid bodies and general discrete mechanical systems-Painlevé and Kane paradoxes revisited[END_REF].

Formerly, the general and rather imprecise idea which was prevailing is that one should use the Lagrange equation for dynamics, complemented with the Coulomb friction law applying to real world reaction forces and to the velocities of the material contact points. This paradigm has not led so far to any general formulation of frictional multibody dynamics. It has only resulted in two outcomes:

-In very specific cases, there is a one-to-one correspondence between generalized reaction forces and real world reaction forces, and this paradigm can then be turned into a sound mathematical evolution problem. One such example is the Painlevé example of a rigid bar above a rigid obstacle, so that there can only be, at most, one material point at which the contact takes place. Such a precise formulation of the evolution problem can be found in Stewart [START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of painlevé's problem[END_REF] in the case of the Painlevé example. The motivation of Stewart was to prove the existence of a dynamical solution for the Painlevé example (which he did), since suspicion of non-existence of solution had fed a one-centuryold controversy which could be entitled: 'is the (local) Coulomb law consistent with rigid body dynamics?'. -Numerical methods to compute discrete-time solutions in more complex cases. Examples of such methods are those of Jean [START_REF] Jean | The Non-Smooth Contact Dynamics Method[END_REF], which was inspired by the time-stepping technique of Monteiro Marques, and the method of Stewart in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of painlevé's problem[END_REF]. A common feature of all these methods is that they are derived on the basis of the paradigm of the Coulomb law applying to reaction forces in the real world. But, they are not derived on the basis of a general continuous-time evolution problem. This was recognized by Stewart in [START_REF] Stewart | Rigid-Body Dynamics with Friction and Impact[END_REF]: 'In many ways it is easier to write down a numerical method for rigid-body dynamics than it is to say exactly what the method is trying to compute.' As was already pointed out when formulating the evolution problem associated with frictionless multibody dynamics, the equation of motion involves the generalized reaction force which is therefore completely determined, when the motion is calculated, whereas the precise distribution of the real world reaction forces remains undetermined, in general. Besides, the modern point of view in mechanics stipulates, through the systematic use of the Principle of Virtual Power, that the appropriate representation of forces within any mechanical theory is the one that comes from duality with velocities. Within the dynamics of rigid bodies, the appropriate representation of forces is therefore that of Lagrange's generalized forces. Hence, it is not surprising that sticking to a formulation of the friction law in terms of the real world reaction forces has not yet come to a successful end, with respect to the issue of obtaining a general consistent formulation of the dynamics. In the specific cases where it is possible though, such a paradigm is known to raise the following paradoxes, once again unsurprisingly.

-Painlevé paradox. In the case of the planar dynamics of a rigid bar in a half-plane with possible frictional contact at one extremity of the bar, there is a one-to-one correspondence between generalized reaction force and real world reaction force. It is therefore possible to obtain a precise formulation of an evolution problem in that situation (it is to be found in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of painlevé's problem[END_REF]). However, it was recognized as early as the nineteenth century that multiple solutions can be encountered [START_REF] Painlevé | Sur les lois du frottement de glissement[END_REF]. The example of the rigid bar would therefore be different in nature of that of a point particle for which the uniqueness of the solution can be proved rigorously (see the discussion in section 3). At the time of Painlevé, there was also a suspicion of possible non-existence of a solution. Later, Lecornu [START_REF] Lecornu | Sur la loi de Coulomb[END_REF] pointed out that the existence of a solution could be recovered by allowing an impact (velocity jump and associated Dirac mass of the reaction measure), although the incoming normal velocity vanishes. This was further discussed by Moreau [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], in the framework of the motions with measure acceleration, under the name 'tangential impact' or 'frictional catastrophe'. Finally, Stewart proved in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of painlevé's problem[END_REF] an existence result for the Painlevé example with completely inelastic impact law, showing that the allowance for 'tangential impacts' was ruling out the suspicion of non-existence in the Painlevé example. Hence, the analysis of Stewart, which was restricted to the system with one material contact point only, solved the non-existence issue in the Painlevé paradox, but let open the non-uniqueness issue. -Kane paradox. The non-uniqueness part of the Painlevé paradox was sometimes claimed to be acceptable. However, another paradox connected to the use of a friction law applying to real world reaction forces arose more recently in [START_REF] Kane | A dynamics puzzle[END_REF][START_REF] Kane | Dynamics : theory and applications[END_REF], which was more serious. The example is that of a double pendulum whose extremity can experiment frictional contact with a straight obstacle. The Coulomb law and an impact law with restitution coefficient of 0.5 are adopted. In that example also, there is an exceptional one-to-one correspondence between generalized and real world reaction forces, making it possible to formulate an evolution problem associated with the Coulomb law applying to real world reaction forces. A motion solving that evolution problem is exhibited, in which the kinetic energy of the system is increased during an impact, due to friction. This is an energetic inconsistency raised by the use of the Coulomb law applying to real world reaction forces within rigid body dynamics. In the case of the Painlevé example with completely inelastic impact law, Stewart proves in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of painlevé's problem[END_REF] that the solution he constructs is dissipative, and therefore does not exhibit this energetic inconsistency. However, since no uniqueness is to be expected, there is no guarantee of energetic consistency of all the solutions.

The case of one unilateral constraint

In this section, we shall obtain the evolution problem governing the frictional dynamics of a discrete system obeying only one unilateral constraint of the form:

ϕ(q) ≥ 0.
According to the discussion in the preceding section, we are going to resist invoking any real world reaction force and therefore stick to generalized reaction forces and velocities.

Splitting into normal and tangential parts

At a configuration where the constraint is active, it is easy to see how to define what a tangential velocity is: it is simply a velocity that arises in a motion maintaining the contact active: dϕ q , q = 0.

However, the issue of defining what a normal velocity is, turns out to be trickier, as it requires the choice of a scalar product. Even if the configuration manifold can locally be identified with R d by making the choice of a parameterization, the use of the canonical scalar product of R d is not an option, for, in that case, the definition of what a normal velocity is, would depend on the choice of the parameterization. Actually, there is only one natural scalar product which is intrinsic, that is, does not depend on the choice of a parameterization: the one induced by the kinetic energy whose matrix in a given parameterization is the (positive definite) mass matrix M(q). We shall see later that this inevitable choice of a scalar product has the fallout of enforcing energetic consistency and therefore protecting from Kane paradox. A normal velocity v n is therefore one which is orthogonal to a tangential velocity according to the scalar product

M(q): λ ∈ R, v n = λ M -1 (q) • dϕ q = λ ∇ q ϕ,
where the scalar product M(q) enables, as usual, to make the link between the gradient and the differential of the function ϕ. Likewise, a normal generalized reaction force r n is one that develops zero virtual power in any tangential velocity:

λ ∈ R, r n = λ dϕ q ,
and a tangential reaction force r t is one that develops zero virtual power in any normal velocity:

r t • M -1 (q) • dϕ q = 0.
Some insight is gained from a(n elementary) geometric point of view: velocities live in the tangent space T q Q which is endowed with an Euclidean structure by the scalar product (mass matrix) M(q). Generalized forces are linear forms applying on velocities through the Virtual Power paradigm, and therefore live in the dual space T * q Q which is endowed with a natural Euclidean structure by the dual scalar product M -1 (q). Hence, the natural split into tangential and normal parts read as:

∀v ∈ T q Q, v = v t + v n ∇ q ϕ
|∇ q ϕ| q , with v t , ∇ q ϕ q = 0, and |∇ q ϕ| q = dϕ q • M -1 (q) • dϕ q , ∀r ∈ T * q Q, r = r t + r n dϕ q |dϕ q | * q , with r t , dϕ q * q = 0, and |dϕ q | * q = dϕ q • M -1 (q) • dϕ q .

The generalized friction law

We are now going to rely on the splitting of generalized reaction forces and velocities into normal and tangential parts to infer the general form of a friction law expressed in terms of these quantities.

Let us first recall some classical definitions of convex analysis. Let C be a nonempty closed convex subset of R d . Its indicatrix function I C is defined by: The pointwise formulation (2) of the Coulomb law is therefore equivalent to any one of the following equivalent statements:

I C (x) = 0 if x ∈ C +∞ if x / ∈ C It is
(i) -v t ∈ ∂I µrnB r t (ii) r t ∈ ∂S µrnB -v t (iii) ∀v ∈ R 2 , r t , v -v t + µr n |v| -|v t | ≥ 0,
where ∂f denotes the subdifferential of the proper, lower semi-continuous function f . Since ∂I C is the closed convex cone of all the outward normal vectors to C , statement (i) is also called the normality rule, when µr n is assumed to be given as in Coulomb's experiments (where r n was merely the opposite of the gravity force). The equivalent statement (ii), is often referred to as the maximum dissipation principle.

To generalize the Coulomb friction law (2) for the point particle to the case of an arbitrary discrete mechanical system, we are naturally led to postulate the existence of a set C * q in the space T * q ⊂ T * q Q of all tangential reaction force, such that all the admissible tangential components of the generalized reaction force are characterized by the condition:

r t ∈ r n C *
q , where the normal reaction force r n is required to be nonnegative and vanishes when the constraint is not active (that is, when ϕ(q) > 0). As for the Coulomb friction law (2), this condition entails that whenever the normal reaction force vanishes (for example in motion episodes where the constraint is not active), then the tangential reaction also vanishes.

In general, the set C * q will not be required to be a disk (case of isotropic friction (2) for a point particle), but only to be nonempty. It is also required to be closed and convex so that the normality rule can be adopted as a flow rule: this is nothing but assuming that for a given value of the reaction normal component, the flow rule should obey the maximum dissipation principle. With this assumption, the flow rule is formally expressed by any one of the following equivalent statements:

(i) -q+ t ∈ ∂I rnC * q r t (ii) r t ∈ ∂S rnC * q -q+ t (iii) ∀v ∈ T q , r t , v -q+ t + r n S C * q (-v) -S C * q (-q+ t ) ≥ 0,
where the identity S rnC * q = r n S C * q has been used. Here, the subdifferentials should be understood in the sense of the duality between the space T q of tangential generalized velocities and the space T * q of tangential generalized forces.

Actually, the set C * q can be viewed either as a subset of the space T * q of tangential generalized forces or as a subset of the space T * q Q of all generalized forces. In the latter case, one has: ∀v ∈ T q Q, S C * q (v) = S C * q (v t ), since C * q is contained in the space T * q of tangential generalized forces. Therefore, statement (iii) is equivalent to:

(iii ) ∀v ∈ T q Q, r t , v -q+ + r n S C * q (-v) -S C * q (-q+ ) ≥ 0, where we recall that T q Q stands for the space of all generalized velocities at the generalized configuration q (tangent space at q to the configuration manifold). This latter form of the (generalized) friction law is -∀v ∈ C 0 [0, T ]; T Q , with v(t) ∈ T q(t) Q [0,T ] r t , v -q+ + r n S C * q (-v) -S C * q (-q+ ) ≥ 0, ϕ q(t) = 0 ⇒ q+ n (t) = -e qn (t). The following results about problem P 4 have been proved in [START_REF] Charles | The formulation of dynamical contact problems with friction in the case of systems of rigid bodies and general discrete mechanical systems-Painlevé and Kane paradoxes revisited[END_REF].

-If e ∈ [0, 1], then any solution of problem P 4 is dissipative, that is, for all t 1 , t 2 ∈ [0, T ], such that t 1 ≤ t 2 , one has:

1 2 q+ (t 2 ) • M q(t 2 ) • q+ (t 2 ) ≤ 1 2
q-(t 1 ) • M q(t 1 ) • q-(t 1 ) + t2 t1 f t; q(t), q-(t) , q-(t) dt. [START_REF] Percivale | Uniqueness in the elastic bounce problem, ii[END_REF] Problem P 4 is the general evolution problem raised by the frictional dynamics of an arbitrary discrete system submitted to one unilateral constraint. The impact law governs the normal component of the velocity, as in the frictionless situation. The generalized friction law rules both the impacts and the smooth episodes of motion. The above result shows that any solution is energetically consistent. Let us just recall that the example of Kane in [START_REF] Kane | A dynamics puzzle[END_REF][START_REF] Kane | Dynamics : theory and applications[END_REF] was that of a discrete system (double pendulum), with one unilateral constraint. However, he did not know such a formulation as our problem P 4 , and used the usual Coulomb friction law applying to real world reaction forces. As a result, some solutions were possible in which an increase of kinetic energy during an impact was possible. The above result is therefore a demonstration of superiority of formulating the dynamics along problem P 4 , in which energetic consistency is built-in. -Existence and uniqueness of solution for problem P 4 are still to be proved (under the assumption that the data are analytic). Note, however, that this existence and uniqueness result has already been proved in the case of problem P 2 , which is nothing but a particular case of problem P 4 , for which the friction set C * q reduces to a disk. In addition, a discussion around the Painlevé example is developed in [START_REF] Charles | The formulation of dynamical contact problems with friction in the case of systems of rigid bodies and general discrete mechanical systems-Painlevé and Kane paradoxes revisited[END_REF]. When the dynamics of the Painlevé system is formulated on the basis of the usual Coulomb friction law applying to real world reaction forces (as in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of painlevé's problem[END_REF]), multiple solutions can be exhibited for large enough friction coefficients. However, when the dynamics of the same system is formulated along problem P 4 , among these explicit multiple possible dynamical evolutions, only one of them remains a solution of the corresponding problem P 4 . This is a clue that formulation of the dynamics along problem P 4 escape not only from energetic inconsistencies as in Kane paradox, but also from indeterminacies and multiple solutions as in Painlevé paradox.

The general case

There remains to derive a systematic formulation of frictional dynamics of discrete systems, in the case where several unilateral constraints apply: ϕ α (q) ≥ 0 (α = 1, 2, . . . , n). This is nothing but an extension of frictionless multibody dynamics along problem P 3 , in the spirit of problem P 4 . Such an extension is presented with details in [START_REF] Charles | The formulation of dynamical contact problems with friction in the case of systems of rigid bodies and general discrete mechanical systems-Painlevé and Kane paradoxes revisited[END_REF].

The analysis introduces a collection C * α,q (α = 1, 2, . . . , n) of friction sets, each of them being associated with one of the unilateral constraints. The corresponding generalized friction law involves the sum of each associated dissipation: S C * α,q (-q+ ). A corresponding evolution problem is precisely stated and it is rigorously proved that any solution is dissipative, that is, fulfils inequality [START_REF] Percivale | Uniqueness in the elastic bounce problem, ii[END_REF].

Hence, the point of view of a generalized friction law applying to generalized reaction forces (instead of the usual point of view of a friction law expressed in terms of the real world reaction forces) yields the following pleasant fallout. It enables to formulate for the first time, a consistent abstract evolution problem associated with frictional unilateral multibody dynamics, in the most general situation. And, in

Figure 1 .

 1 Figure 1. Non-continuous dependence of the solution on the initial datum.

  proper, lower semi-continuous and convex. Its conjugate function (by the Legendre-Fenchel transform) is the support function S C of C : Introducing the closed unit disk B of R 2 (with respect to the canonical Euclidean norm | • |), it is readily checked that: S B (y) = |y|.

	S C (y) = sup	x, y .
	x∈C

an instant t is said to be an atom of the Radon measure r, if r({t}) = {t} r = 0.

addition, energetic consistency of such a formulation is built-in, meaning that any solution is necessarily dissipative.

the more appropriate one in view of deriving an expression of the friction law encompassing the episodes of smooth motion and the impacts, that is, in situations where the reaction force is not pointwise defined but is only a measure (with respect to time). Having these situations in mind, the friction law has to be postulated under the form:

As in the case of the ordinary Coulomb law where C * q is a closed disk, the generalized friction law (7) makes sense whenever the generalized reaction force r is a measure and the generalized right-velocity q+ is a function of time with bounded variation. Equation [START_REF] Percivale | Uniqueness in the elastic bounce problem, i[END_REF] provides the general form of what should be a friction law formulated in terms of generalized reaction force, applying to both smooth episodes of motion and impacts, which is, in addition, energetically consistent. This framework cannot be avoided. Of course, the question remains to know what choice should be made in practice for the friction set C * q . It turns out that a systematic proposal of C * q can always be made on the basis of the paradigm of a pointwise Coulomb cone containing hypothetical real world reaction forces and, then, of the use of formulae (3) to build the cone containing the generalized reaction force r, and therefore to deduce a proposal for C * q (a detailed account of this systematic construction is to be found in [START_REF] Charles | The formulation of dynamical contact problems with friction in the case of systems of rigid bodies and general discrete mechanical systems-Painlevé and Kane paradoxes revisited[END_REF]).

In the above-mentioned systematic construction for the sets C * q , it turns out that unbounded sets C * q can appear. In that case, the function S C * q takes the value +∞, that is, its domain D(S C * q ) is smaller than the entire space T q Q R d . Since the generalized friction law [START_REF] Percivale | Uniqueness in the elastic bounce problem, i[END_REF] entails that the tangential right-velocity q+ t must be in the closure D(S C * q ), it turns out that unbounded friction sets induce forbidden values for the right-velocity q+ t at those instants belonging to Supp r. Another significant difference of unbounded friction sets with bounded friction sets is that they allow atoms of tangential reaction force without any atom 1 in the normal force, which is impossible in the case of a bounded friction set. Typically, this is going to happen when the system q is arriving on the obstacle with zero normal velocity and along a forbidden value of the tangential velocity. In that case, a tangential impact (discontinuity of tangential velocities at an instant where the normal velocity is continuous) will be possible.

The evolution problem

Gathering all the preceding considerations, we are in measure to obtain the general formulation of the frictional dynamics of an arbitrary discrete mechanical system with only one unilateral constraint. As previously, we are given an initial condition (q 0 , v 0 ) ∈ T Q R d × R d compatible with the obstacle: ϕ(q 0 ) ≥ 0 and ϕ(q 0 ) = 0 ⇒ v 0 , dϕ q0 ≥ 0, and v 0 ∈ D(S C * q 0

),

(the last requirement is a restriction only in the case where C * q0 is unbounded and aims at protecting from tangential impact at the initial instant). Problem P 4 . Find q ∈ MMA([0, T ]; R d ) and r ∈ M ([0, T ]; R d ) such that:

q(0) = q 0 , q+ (0) = v 0 , -M(q) • D q Dt = f (t; q, q-) + r, -∀t ∈ [0, T ], ϕ q(t) ≥ 0, -Supp r ⊂ t ∈ [0, T ] ϕ q(t) = 0 , r n ≥ 0,