Jean Jacques Moreau

An introduction to Unilateral Dynamics

The paper is devoted to mechanical systems with a finite number of degrees of freedom. After showing how inequality requirements in evolution problems can be handled through differential inclusions, one introduces dynamics by an elementary example of unilateral mechanical constraint. Then a general setting is constructed for multibody multicontact systems. The description of unilateral interaction at each possible contact point is formalized, with account of possible friction. This generates the numerical time-stepping policy called Contact Dynamics. The treatment of collisions or other frictional catastrophes in this framework leads to measure-differential inclusions, an essential tool in nonsmooth dynamics. The energy balance of nonsmooth evolutions is discussed. Two illustrations of the proposed numerical methods are presented. The former concerns the mechanisms of collapse of a bridge arch under local forcing. In the latter, the construction of a conical pile of grains is simulated, in order to investigate stresses in the bulk and the distribution of pressure on ground.

1 Purpose

Computation in multibody dynamics

).

The techniques used in multibody dynamics may be classified into the three following categories.

Event-driven methods

The methods so qualified (abbr. ED) are practical mainly when the con cerned time-interval equals the union of not too many subintervals, a priori unknown, over which the status of the various contacts remain unchanged, i.e. no collision which would create new contacts occurs, no contact either gets loose, nor any critical situation needing a change in the analytical expression of the Coulomb law of dry friction is met. On each of these subintervals, the same numerical techniques as in the investigation of machines with classical bilateral, possibly frictional, constraints may be used [START_REF] Haug | Computer Aided Kinematics and Dynamics[END_REF]. As integration pro ceeds, there only is to watch the evolution of some indicators. In particular, the contact forces will be calculated. If, after a certain instant, some of the values computed for these forces are found to have directions incompatible with the unilaterality of the non-interpenetrability constraints (here we ne glect adhesive, i.e. gluing effects), the programme decides that the motion has to be calculated otherwise. But one should keep in mind that the contacts which get loose after the critical instant are not necessarily those for which an unfeasible contact force has just been evaluated [22]. A popular approach to such discussions consists of reducing them to complementarity problems, similar to what is commonly met in constrained optimization.

More embarrassing is the calculation of the new velocity in case the criti cal instant is that of a collision. Contact forces take very large values during a very short time-interval and only poor phenomenological information is usu ally available about such physical circumstances. Even in the simplest case, that of the collision of two otherwise free members of the system, the tradi tional coefficient of restitution is known to depend not only on the materials these bodies are made of, but also on their shapes and relative orientations at the collision locus [START_REF] Stoianovici | A critical study of the concepts of rigid body collision theory[END_REF]; only the case of spherical beads appears relatively comfortable [START_REF] Foerster | Measurements of the collision properties of small spheres[END_REF]. Furthermore, if some of the colliding bodies belong to clus ters of already contacting ones, percussional reactions should be expected at all existing contact points. Though rigid body collisions currently are an active domain of research [START_REF] Brogliato | Nonsmooth Mechanics[END_REF][13] [START_REF] Cholet | Exper imental study of collisions of angular particles[END_REF][27] [START_REF] Stronge | Contact problems for elasto-plastic impact in multi-body systems[END_REF], computation has to rely on prag matic rules whose validity has to be checked in each domain of application.

Smoothing methods

'Nonsmoothness' is the salient feature of the problems in view. In fact, after the set of the possible positions of the investigated system has been para metrized through an element q of R n , the geometric restriction that the non-interpenetrability constraints impose on q are expressed by a set of in equalities. Hence, instead of running in a smooth submanifold as in tradi tional analytical dynamics, the point q is confined in a region of R n whose boundary is made of a lot of pieces of hypersurfaces ( millions or billions of them in current applications to granular materials): this is nonsmoothness in space. Furthermore, collisions are expected to induce velocity jumps: this is nonsmoothness in time. To end, the contact forces or 'reactions' associated with the non-interpenetrability constraints are governed by highly irregular laws. These forces vanish as soon as the corresponding contacts break while, if contact holds, the commonly stipulated mechanical conditions do not ex press them as fonctions of q. If, in addition, dry friction is taken into account (most usually in the form of Coulomb's law) it introduces some irregular re lationships between contact forces and the sliding velocities. All this may be called nonsmoothness in law.

In such a state of the affairs, a natural move is to replace, approximately, the nonsmooth governing relationships by some regularized ones. First the non-interpenetrability constraints will be replaced by some stiff repulsion laws which take effect as soon as two members of the system corne close to each other. This automatically handles the possible collisions, as far as one con siders them as 'elastic', while the dissipativity of collisions may be accounted for by adding some damping actions or also by using different repulsion laws in the episodes of approach and of separation. Similarly, frictional contact may be somewhat regularized through the introduction of local elastic mi crodeformation and of viscosity-like effects. The dynamics of the approximate system is then governed by differential equations with sufficient regularity to be handled through standard numerical techniques. The drawback is that the need of precision requires the use of very stiff appproximate laws. Hence the time-stepping procedures applied have to resort to very small step-length and possibly also have to enforce numerical stability by introducing artificial damping or artificially increasing inertia. When treating dynamical applica tions, the effect of such an artificial alteration of the mechanical data may blur the picture. Significant simulations of loose ( collisional) flows of granular materials have been obtained in that way, but when dense collections of bod ies are concerned (pieces of masonry or compact granula tes) the method is mainly applied to quasi-static evolutions in which only a succession of equi librium states is looked for. Dynamical computation is then used only as a way of attaining each of these equilibrium states, a strategy referred to as 'Dynamic Relaxation'.

The pioneering work of P. Cundall [START_REF] Cundall | A computer mode! for simulating progressive large scale movements of blocky rock systems[END_REF] was precisely based on the regu larization strategy, today implemented in the majority of commercial pieces of software intended to handle non-interpenetrability. Because such computa tion techniques are close to those applied in molecular simulations, they are commonly referred to as 'Molecular Dynamics' methods (abbr. MD), specially in the domain of granular mechanics [START_REF] Walton | Numerical simulation of inelastic, frictional particle particle interactions[END_REF].

Contact Dynamics

This is the technique (abbr. CD) advocated by the author, to which the present paper is meant to provide an introduction. It originated from [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dy namics[END_REF], where the (unilateral) contact, possibly frictional and/or collisional, between rigid bodies received a formulation in terms of elementary convex analysis which proves suitable for computation. Mathematically, the resulting evolu tion problems are governed in smooth cases by differential inclusions. General information on the latter concept may be found in [START_REF] Aubin | Differential inclusions[END_REF] [START_REF] Deimling | Mémoire sur la théorie des liaisons finies unilatérales[END_REF], but the need of treating also non-smooth evolutions calls for its extension to that of a mea sure differential inclusion. The Sweeping Process introduced earlier by the author, with motivation in the quasi-static evolution of elastoplastic systems [START_REF] Moreau | On unilateral constraints, friction and plasticity[END_REF] [START_REF] Moreau | Evolution problem associated with a moving convex set in a Hilbert space[END_REF] seems to have provided the first occurence of measure differential in clusions in literature. This process and some related evolution problems are still today the objet of mathematical research; see references in [START_REF] Kunze | An introduction to Moreau's sweeping process[END_REF]. ln Sect.2 of the present paper a description of the process, of purely kinematical na ture, has been estimated able to make a tutorial introduction to the handling of unilaterality in mechanics.

A directing idea of the Contact Dynamics approach is that the main ob ject of computation is the velocity function t f---+ u E R n . Time-stepping algorithms essentially have to determine the evolution of this fonction. by applying the principles of dynamics and the specified force laws. The po sition fonction t f---+ q is only to be updated at each step through adequate integration. Possibly, some members of u are 'pseudo-parameters' such as the components of the spin vector of a solid, instead of the time-derivatives of po sition parameters. From the geometrical standpoint, one should observe that, in analytical dynamics, the position q ranges in some differential manifold, while for each q the possible velocities are elements of the tangent space to the manifold at this point. The latter is a vector space relatively to which such concepts as linearity, convexity, etc. make sense, while no algebraic concept of this sort is generally available in the position manifold.

Contact Dynamics procedures rest on drawing the balance of momentum of the investigated mechanical system over each time-step. No estimation of the acceleration is needed and the burden of calculating mathematically the curvatures of the involved surfaces is avoided. These curvatures are sim ply accounted for by the fact that, from one step to the other, the normal directions to the detected contacts vary.

The method results in time-stepping schemes which, at least in what concerns the velocity fonction, are of the imp licit type. For this reason, a single computation step is needed in particular to decide whether equilibri um in a given position is a possible motion.

In principle, the Contact Dynamics algorithms are ready at each step to treat collisions on the same footing as persistent contacts but, of course, the need of physical information about such phenomena is not overcome.

Contents of the paper

The contents of Sect. 2 is purely kinematical, in other words it pertains to differential and integral calcul us in the variables of time and space. There is shown how inequality requirements imposed on a moving point in R n may be enforced through conditions involving its velocity vector. The example of the Sweeping Process, which is governed by a differential inclusion of first order, provides an introduction to the mathematical aspects of Unilateral Dynamics.

In Sect. 3 is developed the elementary example of the dynamics of a punc tual particle unilaterally confined by a frictionless material boundary with prescribed motion. The traditional aim of rational mechanics, that of elim inating the unknown reactions of frictionless constraints, is attained in the form of a differential inclusion. In this simple mechanical example, the decisive step leading to the Contact Dynamics method is introduced. It consists in re placing the familiar formulation of frictionless contact by an apparently more complicated relationship involving velocities. The advantage of this transfor mation is illustrated by its application to a time-stepping scheme for the numerical approximation of the motion.

The framework for treating multibody multicontact systems is constructed in Sect. 4. This starts, as usual, with the parametrization of the system posi tions by means of an element q of R n , after what the geometrical conditions of non-interpenetration of the system members and their possible confinement by external obstacles are assumed expressed by a finite set of inequalities. A contact corresponds to one of these inequalities being satisfied as equality. At every instant such that the derivative dq/dt = u E R n exists, every possible contact involves a vector U of local relative velocity of the contacting abjects which, generically, is expressed as an affine fonction of u. As far as the contact actions consist of a simple force R, calculating the element r of R n made of its generalized components is only the matter of matrix transposition. If the traditional assumption of differentiability of t f----t u is made, the equations of the system dynamics may be written clown through Lagrange's formalism or any other method of stereodynamics.

Exploiting these equations requires some information about the possible contact forces. Such an information receives a codified form in Sect. 5 under the name of a contact law, generically a relationship between R and U, usually depending on the actual configuration attained by the system. Here again, the special case of frictionless contact is taken as an introduction. It allows one to stress the importance of the concept of a 'prospective' contact law. Roughly speaking such a law, rather than connecting the actual values of R and U, pertains to their limits on the right of the concerned instant. The Coulomb law of dry friction receives an adaptation to such a formalism.

Then cornes in Sect. 6 the application of the preceding to the construc tion of a time-stepping numerical scheme for the numerical approximation of solutions, first developed in the traditional framework of smooth motions, i.e. the unknown fonction t f----t u is assumed locally absolutely continuous. This scheme is of the implicit type with regard to u. Due to the form given to contact laws and to their discretized forms, this time-stepping procedure automatically secures the preservation of the non-interpenetration conditions and manages the possibility of contact breaking.

In contrast, the event of the sudden introduction of new contacts, i.e. collisions, leads in Sect. 7 to giving up the framework of smooth dynamics, since u is expected to exhibit jumps. Henceforth, the fonction u is assumed to be of (locally) bounded variation on the considered time-interval. With such a fonction, there classically is associated an R n -valued measure du on the time-interval, which may be called the differential measure of u. Smooth dynamics is retrieved as the special case where du possesses a density fonction with regard to the Lebesgue measure dt.

The measure du presents an atom at every instant of discontinuity of u, but this fonctional framework may also accomodate other sorts of frictional catastrophes, as referred to in 7.1. Instead of a classical differential equation, dynamics is now governed by a measure-differential equation. Contact actions are not necessarily represented as time-dependent forces but as impulsion measures. Their integrals on a time-interval constitute impulsions, a familiar concept in traditional dynamics.

Since the time-stepping numerical procedure of Sect. 6 consists of drawing the balance of impulsions on each time-step, it remains formally applicable to non-smooth evolutions. There does not seem objectionable to apply the precedingly identified contact laws so long as the measure du is diffuse, i.e. u continuous. The case of a collision occuring at some instant t e is more em barrassing since each contact-impulsion throughout the system is liable to possess an atom at point t e . The weight of such an atom is nothing but the vector traditionally called a percussion, say P. In 7.5 an efficient computation trick is introduced, by which P is connected, through Coulomb's law in the prospective form, with some weighted mean of the pre-and post-collision val ues of the local velocity U. The resulting calculation is definitely richer than the traditional laws of restitution in that all contacts present in the system at instant t e are collectively taken into account, but its phenomenological validity has to be tested in each application context.

N onsmooth dynamics does not allow one to draw energy balances as pre cise as those traditionally available for smooth evolutions. The reason is that some rules of the differential and integral calculus have to be replaced by a calculus of differential measures, shortly presented in Sect. 8. The thermody namic correctness of collision models requires dissipativity, a property which is not always confirmed by energy balances, as drawn there.

As final illustrations, two examples of the numerical application of the foregoing are presented.

In Sect. 9 a two-dimensional model of a stone bridge is used to show how the collapse mechanism under some localized forcing depends on the value admitted for the interblock friction coefficient.

Section 10 is devoted to the numerical simulation of the construction of a conical pile of grains. Eager controversy took place in recent years about the experimental finding of a local minimum of ground pressure at the vertical of the pile apex. The numerical exploration of stress in a numerical simulation conducted by the CD method appears to remove any paradoxical features from this subject.

The differential handling of inequality conditions

Expressing viability in terms of velocity

Let us understand by a moving point a mapping q of a ( time) interval I into R n . We want q(t) to comply for every ti n I with some inequality

f(t,q(t)) s 0, ( 1
)
where f : I x R n --+ R denotes a C 1 real fonction. In other words, q(t) is required to belong for every t to tf>(t

) := {x E R n I f(t,x) S O}.
It is assumed that, fort in I and x in R n , the gradient Vf(t,x)

:= ( â f / âx 1 , ... , â f / âx n ) is a nonzero n-vector.
Let some t be such that the right-side derivative q'+ ( t), called the right side velocity of the moving point, exists. Through the chain rule, the real fonction T r-+ f ( T, q( T)) is found to possess at T = t a right-side derivative equal to ff ( t, q( t)) + q'+ (t). V f (t, q( t) ). This derivative should be S O if [START_REF] Abadie | Dynamic simulation of rigid bodies : modelling of frictional contact[END_REF], assumed to hold everywhere in I, is satisfied at t as equality. In contrast, if inequality holds strictly at t, no sign condition cornes to restrain right-side derivatives.

For t in I and x in R n , put I'(t x) := { {v E R n I ff � t,x) +v.Vf(t,x) S O} if f(t,x);::: 0 ' R n otherw1se. so that the above observation means q'+(t) E I'(t, q(t)). What follows may be viewed as providing a converse. Suppose that the interval I, nonnecessarily compact, contains its origin t0 and that q is locally absolutely continuous on I. Equivalently, the (two side) derivative dq/dt exists almost everywhere in I and equals a fonction u : I --+ R n , the velocity function, which is locally integrable with regard to the Lebesgue measure on I; notation : u E .CfojI; R n ), meaning that u is Lebesgue-integrable over every compact subinterval of I. And, for every t in I, one has

q(t) = q(to) + l t u(s) ds. to (2)
The following is established in [START_REF] Moreau | Sorne basics ofunilateral dynamics[END_REF]:

Viability Lemma. Assume that q is locally absolutely continuous on I and that the inclusion

dq dt E I'(t, q(t)) ( 3
)
holds for almost every t in I. If inequality ( 1) is verified at the initial instant t 0 , it is verified for every subsequent t.

The term viability is used in system theory and control, in particular when economic systems are concerned, to express that the trajectory of a process should remain in a specified set (one may refer to [START_REF] Aubin | Viability Theory[END_REF], a book actually developed in a too sophisticated topological context for such an elementary statement as the above Lemma to find place in it).

Selectors

A condition of the form (3) is called a differential inclusion [START_REF] Aubin | Differential inclusions[END_REF] [START_REF] Deimling | Mémoire sur la théorie des liaisons finies unilatérales[END_REF]. Since the right-hand sicle is set-valued, uniqueness is a priori not expected to hold for the solutions (if any) to initial value problems. By a selector of the 'multifonc tion' (t, x) H I'(t, x), one means a single-valued fonction, say (t, x) H ,(t, x), such that ,(t, x) E I'(t, x) for every t and x. Then

dq dt = ,(t, q(t)) (4)
is a differential equation whose (locally absolutely continuous) solutions, if any, consequent to some initial condition verifying (1), meet the assumptions of the Viability Lemma, making q(t) belong to <I>(t) for every subsequent t.

A basic example is provided by taking as I the 'lazy selector' of I', i.e. by defining ,(t, x) as the element of I'(t, x) with minimal Euclidean norm. When f ( t, x) < 0, i.e. when x belongs to the interior of <I>( t), the set I' ( t, x) consists of the whole of R n , so that 1 ( t, x) equals 0, the zero of R n . Otherwise, I'(t, x) equals a half-space which contains O if ff (t, x) s; 0, in which case 1 = 0 again. If JI(t,x) > 0, one finds 1 = -(ff/llv'fll 2 )v'f, a vector oriented in the direction of decreasing f, normal at x to the hypersurface f(t, .) = Const drawn through this point.

For such a choice of 1 , a solution to ( 4) consequent to some initial position q(t 0 ) in <I>(t 0 ) may be described as follows. The point q(t) belongs for every t to the moving region <I>(t). As long as it lies in the interior of <I>(t), q stays at rest. It is only when the boundary of <I>(t), i.e. the hypersurface with equation J(t, .) = 0, moves inward and reaches q that the point takes on a velocity in inward normal direction, so as to go on belonging to <I>(t). The magnitude of the velocity vector equals the 'normal speed' of the hypersurface.

We have proposed to call Sweeping Process the above kinematical associ ation of point motions to the given motion of a set (in R n or in a real Hilbert space). See [42][59] for references and recent developments of the subject.

Characterizing the process by a differential inclusion

If, at time t, a point x lies on the hypersurface f(t, .) = 0, the vector v' J(t, x) ( we have assumed it nonzero) normal to this hypersurface is directed outward of the region <I>( t). The half-line emanating from the origin of R n , generated by v'f(t,x), is said to constitute the (outward) normal cane to <I>(t) at point x; notation: N P ( t)(x). The definition of a normal cone for less regular sets than <I> is a vast subject; some other cases will be met in the sequel. For x in the interior of <I>(t), it proves consistent to view N P ( t )(x) as reduced to the zero of R n , while the cone shall be defined as empty if x (/. <I>( t). By discussing the various cases occuring in the calculation of the lazy selector, one sees that if 1( t, x) equals this selector, every solution q to the differential equation ( 4) verifies, for almost every t, the differential inclusion dq -dt E N P ( t)(q(t)).

(

) 5
Unexpectedly the converse is true, i.e.( 5) in spite of its multivalued right hand side actually is equivalent to the differential equation ( 4), as far as locally absolutely continuous solution are concerned.

In fact let q : I -+ R n , be such a solution to [START_REF] Aubin | Differential inclusions[END_REF]. For almost every t, the two-side derivative q' = dq/dt exists, so that the right-hand member is non empty and therefore q(t) E <I>(t); the same is true for every t, by continuity.

For t such that q( t) happens to lie in the interior of <I>( t), ( 5) implies q' = 0, which makes that (4) is also satisfied. Otherwise, suppose that q(t) belongs to the boundary, i.e. the fonction T r--+ f ( T, q( T)) vanishes at T = t. Then the right-derivative Jf(t,q(t)) + q'+(t).v'f(t,q(t)), if it exists, is :S O while, symmetrically, the left-derivative is 2:: O. Therefore q' ( t), when it exists, sat isfies Jf (t, q(t)) + q'(t).v' f(t, q(t)) = 0, i.e. it belongs to the boundary of the half-space I'(t, q(t)). Furthermore, (5) entails that q'(t) is directed along the inward normal to the half-space. All this elementarily characterizes q' ( t) as the proximal point to Oin I'(t,q(t)), namely 1(t,q(t)).

It was under the formulation (5) that the Sweeping Process was primi tively introduced [START_REF] Moreau | Problème d'évolution associé à un convexe mobile d'un espace hilbertien[END_REF][50], with <I>(t) denoting a nonempty closed convex sub set of a real Hilbert space H. The motivation then was in the quasi-static evolution of elastoplastic systems [START_REF] Moreau | On unilateral constraints, friction and plasticity[END_REF] [START_REF] Moreau | Application of convex analysis to the treatment of elasto plastic systems[END_REF]. The convexity assumption allows one to establish the existence of solutions under rather mild conditions con cerning the evolution of <I>(t), even discontinuous. Another consequence of this convexity is that the multifonction x r--+ N P ( t)(x) is monotone in the following sense (see e.g. [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contrac tion dans les espaces de Hilbert[END_REF]) : whichever are x1,x2 in H, Y1 in NP(t)(x1), Y2 in NP(t)(x2), one has (x1 -x2).(y1 -y2) 2:: 0, with the dot denoting the scalar product of H. By elementary calculation, this inequality entails that, if t r--+ q 1 (t) and t r--+ q2(t) are two solutions to (5), the Hilbert distance llq 1 -q2II is a non-increasing fonction of t. From this non expansion property, it follows that at most one solution to (5) can agree with some initial position q(t 0 ).

Another source of interest of the formulation (5) is to render evident that the successive positions of the point q are connected with those of the given region <I> in a rate-independent way. In fact, because the right-hand member is a cone, the differential inclusion is found invariant under any non-decreasing differentiable change of variable.

lmplicit versus explicit time-stepping

Coming to the numerical approximation of solutions through time-stepping schemes, let us denote by [t i , t f ], with length h, a time-step ('i' as in initial, 'f' as in final). From an estimate q i of q(t i ), obtained as the result of the antecedent time-step, computation has to deliver an estima te qr of q( tr).

The formulation ( 4) naturally leads to take ui = 1( t i , qi) as an estima te of the velocity throughout the time-step, generating the prediction qr = qi +hui. This is a computation scheme of the explicit type.

If ( 5) is discretized by viewing ( qr -q i ) / h as a representative of the velocity, a strategy of the explicit type would not allow one to express qr, since the right-hand member is multivalued. In contrast, the implicit strategy consists in invoking the value that this right-hand member would take at the unknown point, so one has to solve (6) (the positive factor h has been dropped since N <l> (t f ) is a cone). This qualifies qr as an orthogonal projection of q i onto <I>( tr). In the case where <I>( tr) is convex, the projection is unique and this characterizes qr as the nearest point to q i in <I>( tr). In particular qr = q i when q i happens to belong to <I>( tr). We have proposed to call this procedure the catching-up algorithm [START_REF] Moreau | Evolution problem associated with a moving convex set in a Hilbert space[END_REF].

Complementarity

From the description made of the Sweeping Process in 2.2 it is clear that velocity may be discontinuous. As for the explicit prediction qr = q i + hui, it only requires of Ui to be the derivative of the fonction q on the right of t i . The following observation makes an introduction to analogous, but more complicated, situations we are to meet in Dynamics.

Let q, associated with u through (2), verify [START_REF] Aubin | Differential inclusions[END_REF] almost everywhere in /. Let t 1 E / and assume that the fonction u possesses a limit on the right of t 1 , say ut; in view of [START_REF] Acary | Contribution à la modélisation mécanique et numérique des édifices maçonnés[END_REF] this limit also provides the right derivative cj+(t 1 ). As a stronger assertion than ( 4), which pertained to bilateral derivative and was only declared to hold almost everywhere, let us prove that ut= 1(t1 ,q(t i )).

Trivially, both members of this equality are

O if fi •-f(t 1 ,q(t i )) < O.
Otherwise, i.e. when fi = 0, it was seen that [START_REF] Baraff | Issues in computing contact forces for non-penetrating rigid bodies[END_REF] Inclusion ( 5) means the existence of a fonction t H ,\(t) :s; 0 such that u(t) = >-(t)Vf(t,q(t)). Since Vf is continuous and nonzero, the assumed existence of ut secures that of the right-limit >-t and ut= >t V J (t 1 ,q(ti)).

(

) 8
If i: < 0, instant t 1 is followed by an interval throughout which f < O.

This has been observed to imply u = 0, so that À vanishes on this interval and consequently also its right-limit .Xt. Summing up, one has f.

+ < 0 \ + < 0 1• 1 + \ + 1 = 0, 1 -, "1 -, A (9 )
a system of complementarity conditions. This is a popular formalism in many research domains where inequality requirements are faced. Solving (8) ( 9), with Ïi defined as in [START_REF] Baraff | Issues in computing contact forces for non-penetrating rigid bodies[END_REF], constitutes a linear complementarity problem.

Through arguments from Convex Analysis, such problems are shown to be equivalent to finding the critical points of some quadratic fonctions over poly hedral convex sets. In the present setting, where inequality CJ.i � 0 simply ex presses that ut belongs to the half-space I'( t i , q( t 1 )), one readily checks that the above system of conditions characterizes ut as the minimizing point of the fonction x H llxll 2 /2 over I'(ti, q(t1)), namely 1 (t1, q(t i )) as announced.

A hydromechanical illustration

Assume that n = 2 and that t, x 1 , x 2 are Cartesian coordinates in physical space, with the t-axis vertical and directed downward. Picture the region f (t, x 1 , x 2 ) � 0 as an underground cavity and the curve x 1 = q 1 (t), x 2 = q 2 (t) as a stationary waterstream dripping down into it. The differential equation ( 4), on account of the diverse circumstances met in the definition of the lazy selector, expresses that: i) any part of this stream which happens to be detached from the cavity wall is rectilinear and vertical; ii) when water runs on the wall, it follows a line of steepest descent (this agrees with hydrodynamics under the simplifying assumption that inertia effects are negligible with regard to gravity and to liquid/wall friction); iii) the dependence of I on the sign of f f makes that the stream can run only on a part of the wall exposed upward: when it reaches the rim of a possible overhang, water gets loose and falls vertically down as described in i).

In this example, under the complication typically added by unilaterality, the comparison of ( 4) and ( 5) merely reflects the classical equivalence between the two standard properties of the lines of steepest descent in a surface: at each point on such a line i) the slope is maximal; ii) the direction is orthogonal to the level curve of the surface.

Frictionless confinement of a particle

Primary formulation

Notations are the same as in Sect.2, with n = 3. The element q := (q 1 , q 2 , q 3 ) now consists of the orthonormal coordinates of a material point Q, with mass m, moving under the action of a given force field (t, x) H X(t, x) and constrained in the region '1> ( t) by the impenetrability of its boundary, assumed to be realized as a material surface with imposed motion. This is the occasion of recalling that, in Mechanics, the description of a constraint never reduces � as it does, for instance, in Optimization or in Abstract System Theory � to imposing a geometrical restriction on positions. Sorne information is needed about the mechanical process through which this restriction is enforced. For instance, using some servomechanism in order to secure (1) could result in quite different motions than those obtained under the present assumption of confinement by the contact action of a given material boundary.

On a time interval I throughout which the motion is smooth enough for the velocity function u in (2) to be locally absolutely continuous ( this precludes collisions, events to which we shall corne back later), the motion of the particle Q obeys, almost everywhere in I,

du m dt = X(t, q(t)) + r(t), ( 10
)
where the force r := (r 1 , r 2 , r 3 ) denotes the unknown reaction possibly ex erted by the confining boundary. In this Section, the confinement process is assumed to comply at every instant with the following model:

• this is a contact process f(t,q)<O =} r=O, (11)
• the possible contact is frictionless f(t,q)=O =} :3>..ER:r=>..VJ(t,q), (

• without adhesion

>..«;O. ( 13
)
If compared with the definition given in 2.3 for the normal cone N P (t) (x) at a point x, the above system of three conditions is found equivalent to

-r E NP(t)(q). ( 14
)
Therefore, the traditional aim of Analytical Dynamics, namely the elimina tion of the unknown reactions of the so-called perfect constraints is attained by rapproaching conditions (10) and ( 14)

du X(t,q(t)) -m dt E NP(t )(q(t)), ( 15
)
an integro-differential inclusion, as the unknown fonctions q and u are es sentially connected by [START_REF] Acary | Contribution à la modélisation mécanique et numérique des édifices maçonnés[END_REF]. Formally, this inclusion implies q(t) E <P(t) for almost every t in I, since otherwise the right-hand member would be empty. By continuity q( t) E <P( t) holds for every t.

The 'Contact Dynamics' approach

Ftom the theoretical standpoint as well as in the elaboration of approximation schemes, the differential inclusion (15) proves difficult to handle (see however [START_REF] Paoli | Mouvements à un nombre fini de degrés de liberté avec contraintes unilatérales: cas avec perte d'énergie[END_REF], under the assumption of convexity for the function f). In what concerns the position function q, the problem at hand is of differential order 2, as expected in a dynamical context, while in the purely kinematical setting of Sect. 2, q was the unknown of a differential problem of order 1. However, some of the arguments used then will be transposed in what follows.

Let I' be defined as in 2.1. Provided the initial position q( t 0 ) lie in <P( t 0 ), the solutions of (15) are the same as those of du

X(t, q(t)) -m d t E Nr ( t ,q ( t )) (u(t)) . ( 16
)
In fact ( 16), assumed to hold for almost every t in I, implies u(t) E I'(t, q(t)). Since q(t 0 ) E <I>(t 0 ), the Viability Lemma secures q(t) E <I>(t) for every tin I. The definition of I' then makes that the normal cone N r ( t ,q ( t ) )(u)

is contained in N P ( t ) (q(t)) whatever is u, so that ( 15) is a fortiori satisfied.

Conversely, let q be a solution to [START_REF] Cholet | Chocs de solides rigides[END_REF], hence q(t) E tf>(t) for every ti n I.

When f(t,q(t)) < 0, the set N P ( t ) (q(t)) reduces to the singleton {O} and the same is true for the right-hand member of ( 16). In contrast, for t such that f(t,q(t)) = 0, the set I'(t,q(t)) is a half-space and the argument used in 2.3 proves that the element u(t) = u + (t) = u-(t) belongs to its boundary plane.

Consequently, N r ( t , q ( t )) ( u) consists of the cone generated by v7 f (t, q(t)) hence equals NP(t ) (q(t)).

First example of a CD numerical scheme

As before, let us denote by [t;, trl, with length h, a time-step. Ftom the ap proximate values q;, u; obtained for q and u at t;, one has to calculate qr, ur, pertaining to tr. The given force field X is assumed to depend smoothly on its arguments, so that one chooses to approximate it throughout the time-step by the value it takes at t m := t; + h/2 and q m := q; + hu;/2. It is also at the point (t m , q m ) that f is calculated, in order to decide whether boundary contact is in effect or not and to determine the set I' accordingly. Depending on the sign of f(t m , q m ), the latter equals the whole of R 3 or a half-space with v7 f(t m , q m ) as normal vector. Inclusion ( 16) is thus discretized in the form i.e. in view of the right-hand member being a cone,

h U; + -X(t m , q m ) -Uf E Nr ( t q i (ur). m m, m (17)
This classically characterizes ur as the proximal point to the known element Ui + hX(t m , q m )/m in I'(t m , q m )-One finishes the calculation with This algorithm automatically handles the possible breaking of contact this happens if Ui + hX(t m , q m )/m falls into the interior of I'(t m , q m ) .

Remark 1. -Provided q(to) E <f>(t 0 ), inclusion ( 16) has been precedingly ob served to secure f(t,q(t)):::; 0 at every consequent t. But, in the above time stepping procedure, it is only ur which, at each step, is constructed as an element of I'(t m , q m ), so that the Viability Lemma is just involved through time-discretization. One thus may fear that some violation of the inequal ity would build up from step to step. Actually, if the step-length is not too large, numerical experiments show some self-corrective effect which, on the contrary, tends to reduce violations. This effect seems related to <f>(t) possess ing a nonempty interior. In contrast, in industrial softwares devoted to the dynamics of machines, some bilateral constraints are commonly introduced, leaving a set of feasible positions with empty interior. Their treatment in terms of velocities [START_REF] Haug | Computer Aided Kinematics and Dynamics[END_REF] then requires corrective procedures to prevent cumu lative errors.

Remark 2. -The calculation of ur from ui is based on mechanical elements evaluated at the mid-position q m ; in turn, ur is used to calculate from q m the final position of the current step and, from there, the mid-position of the subsequent time-step. This interleaving makes the above time-stepping procedure resemble the policy called 'leapfrog' in Molecular Dynamics simu lations. In order to figure out how it improves precision, compared with an ordinary Euler explicit time-stepping scheme, one may apply it to calculating the parabolic unconstrained motion of a particle in a uniform gravity field. If h is constant, the values found for q at the successive steps coïncide with the exact solution, while the Euler scheme generates cumulative errors. Of course, when an algorithm of the above sort is implemented with constant h, one may calculate each q m from the antecedent one by a single incrementa tion. The proper output of the computation, namely qr, may not be needed at each step.

Remark 3. -In contrast with the approach of the motion of a point in a surface through the traditional methods of dynamics, the above time-stepping procedure does not require calculating the curvature of the boundary. This curvature is implicitely accounted for by the fact that the direction of v' f evolves from one step to the other. The question of existence of solutions to [START_REF] Cholet | Chocs de solides rigides[END_REF] or ( 16) is not addressed in this paper; one naturally expects that f has to be twice dfferentiable.

Parametrization

Let the possible configurations of a body collection be parametrized (at least locally) through generalized coordinates, say q := (q 1 , ... , q n ). For the sake of reducing the number n, this parametrization may be constructed with ac count of possible permanent, frictionless linkages imposed on the members of the system. After that, the constraints of non-interpenetrability are addi tionally considered. The geometric restriction consequently imparted on the system positions is assumed expressed by a finite set of inequalities fc,(t, q ) S:. 0, a E {l, ... , 1î:}, [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF] where fi, ... , f "' are given fonctions. Through the presence of t as an argu ment off a, provision is made for the inequality to describe the confinement of a member of the system by some external obstacle or boundary with pre scribed motion. U nder this parametrization, a motion of the system consists of a mapping t E / H q(t) E R n and, as in the foregoing, this mapping is assumed locally absolutely continuous, i.e. there exists a locally integrable velocity fonction u : I -+ R n from which q may be retrieved in the form (2).

As an example, one may consider a pair of members of the system whose positions in a chosen reference frame are well located as soon as the value of the element q of R n is known (together with the time t in case of a time dependent parametrization). Then, one may take as fa the expression, as a fonction of (t, q), of some measure of the overlap of the two bodies. This overlap should be understood as a directed quantity, so that it becomes neg ative in case the bodies lie apart from each other. The convention applied in [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF] of characterizing the permitted configurations through the Sc O inequal ity, cornes from Convex Optimization theory, where such a sign convention offers technical advantages. No convexity hypothesis is made here concerning the fonctions f a : since such an assumption would not be preserved under a change of parametrization, it cannot in general have any mechanical mean ing. If one prefers to deal with the � 0 symbol, there only is to consider, instead of the overlap the opposite quantity, usually called the gap between the considered bodies.

The above formalism is not limited to collections of strictly rigid bodies, since q may also include parameters accounting for a finite-freedom approxi mation of deformability. Such additional parameters possibly arise from some modal representation of the deformation dynamics or from the Fini te Element discretization of deformable parts [START_REF] Jean | The Non Smooth Contact Dynamics method[END_REF].

Contact kinematics

Suppose that inequality f a S:. 0 expresses the local non-interpenetration of some pair of members of the system, say B and B ', so that equality f a = 0 corresponds to these bodies touching each other at some point of space denoted by Ma. This we shall assume here to be an isolated contact point, but other contacts, associated with different values of a, may also be in effect between the same bodies at the same instant. For every imagined motion t H q( t) bringing the system through the considered contacting position for some t, with a definite value of u = dq/dt E R n , the velocity vectors Va and V�, relative to the chosen reference frame, of the respective particles of B and B' passing at point M a let themselves be expressed as affine fonctions of u. The same is thus true for the relative velocity Ua = V a -V � of body B with respect to body B' at this point, say [START_REF] Cundall | A computer mode! for simulating progressive large scale movements of blocky rock systems[END_REF] where G a : R n -+ E 3 (the space of the vectors of physical space) denotes a linear mapping, depending on t and q. No attention is paid at this stage to the imagined motion preserving contact or not. The term W a E E 3 , a known fonction of t and q, vanishes in the familiar case of a time-independent parametrization.

Similar formula holds if inequality f a ::; 0 expresses the confinement of a member B of the system by some external material boundary with prescribed motion. Assume that equality f a = 0 corresponds to contact taking place at some point, here again denoted by M a . The local velocity, at this point, of body B with respect to the boundary has the same form as U a in [START_REF] Cundall | A computer mode! for simulating progressive large scale movements of blocky rock systems[END_REF], where W a now reflects the known velocity of the boundary (for a time-independent parametrization, W a equals the negative of this velocity vector).

At the contact point M a , we assume that a common tangent plane to the respective surfaces of the concerned bodies has been defined. This does not require of both surfaces to be smooth; for instance, contact may take place between a smooth body and some corner point or sharp asperity of the other. Let n a denote the unit vector normal to this plane, directed toward B. In computation, as well as in existential studies, it proves usefol that the definition of the above elements would be conventionally extended to a neighbourhood of the concerned value of (t, q) in R x R n . This allows one to express as a fonction of (t,q) the normal gap, say g a (t,q), between B and B', counted as negative in the case of overlap. Classically, the derivative of the fonction t H g a (t,q(t)) is found equal to U a .n °' , the normal component of the relative velocity of the contacting bodies at point M a .

Sometimes, in Computation literature, the second time-derivative of the gap is improperly referred to as the 'normal relative acceleration'. Actually, since the material particles involved in the definition of U a are not the same from an instant to the other, this second derivative has in general nothing to do with the relative acceleration vector. As an example, one may consider a body of circular or spherical shape: a variety of rotations may be imparted to such a body without altering its overall location, so yielding the same gap while the normal acceleration is changed.

Remark. -The representation of non-interpenetration through of a finite set of inequalities is operative in most practical situations. However, non interpenetration cannot be described in that way in the neighbourhood of a confi gu ration where two sharp asperities or corners corne into contact by their points. This is evidenced by the fact that, in this case, the set of the values of the local right-velocity U;; which are compatible with non-interpenetration is not a convex cone anymore. In contexts where the probability of such an event cannot be treated as negligible, numerical techniques have to resort to adequate procedures (possibly involving U a ) for the identification of a mechanically plausible contact plane.

Contact forces

Assume that the contact actions that body B experiences at point Ma from body B' are described as a simple force n a (there would be no conceptual difficulty in adding to this description some local torque, accounting for a resistance to rolling). Then B' experiences from B the force -R n . The stan dard machinery of Analytical Dynamics needs a representation of this pair of forces, in regard to the chosen parametrization, through its covariant compo nents ( or 'generalized components'), namely the element r n of R n expressed as [START_REF] Curnier | A theory of friction[END_REF] with a; : E 3 -+ R n denoting the transpose of G a .

The convention of implicit summation will never be applied to Greek in dices.

If inequality f n :s; 0, expresses the confinement of a member B of the system by some external obstacle with prescribed motion, (19) still holds with U n denoting the local velocity of B relative to this obstacle. Then it is found that r n in (20) consist of the covariant components of the force R n alone, acting on B at the contact point. Its counterpart, exerted by B upon the obstacle, is not in this case a force experienced by the system. Incidentally observe that the term W n does not appear in [START_REF] Curnier | A theory of friction[END_REF].

The equation of Dynamics

As before, the context here is that of standard dynamics, involving the second derivative of the fonction t H q(t). Therefore the velocity fonction t H u(t) is required to be locally absolutely continuous. Using Lagrange's technique or any other tools from classical solid dynamics, one obtains a differential equation, to be read as an equality of elements of R n A(t,q) � � = F(t,q,u) + I)n, [START_REF] Deimling | Mémoire sur la théorie des liaisons finies unilatérales[END_REF] a where A denotes the n x n inertia matrix. The expression F comprises cer tain standard terms (commonly referred to as 'centrifogal' and 'gyroscopic') and also the covariant components of some applied forces, supposed given as fonctions of time, the position of the system and its velocity. The elements , a , a E {1, 2, ... , ,.,, }, are made of the covariant components of the respective contact forces, as expressed in [START_REF] Curnier | A theory of friction[END_REF].

The same formalism remains more generally valid with q related to some velocity fonction u by other kinematical relations than (2). For instance, when dealing with 3-dimensional rigid bodies, it is usual to attach to each of them a frame of principal axes of inertia emanating from its center of mass. Then one may choose to enter, among the constituents of the R n -valued fonction u, the three components relative to these axes of the spin vector of the rigid body, instead of the time-derivatives of some directional parameters. This of fers the considerable advantage of generating a contribution in the matrix A which is diagonal and constant with regard to t and q. Retrieving from these spin components the evolution of some directional parameters of the con cerned body is only the matter of integrating adequate kinematical formulas. Correlatively, if forces are applied to the rigid body, the total moments of these forces about the same axes should be entered as covariant components into the corresponding lines of the right-hand side of [START_REF] Deimling | Mémoire sur la théorie des liaisons finies unilatérales[END_REF].

Since contact forces vanish when contact is not in effect. the summation in the right-hand side of ( 21) may be restricted to the values of a belonging to J(t,q) := {a E {l, ... ,,.,,} 1 f a (t,q) 2: O}.

(

) 22
The geometric conditions (18) of non-interpenetrability, joined to the dif ferential equation of dynamics [START_REF] Deimling | Mémoire sur la théorie des liaisons finies unilatérales[END_REF], clearly are not enough for determining the motion consequent to initial data. Sorne phenomenological information should be added, concerning the contact forces. Since the contact phenom enon takes place in physical space, this information is expected to involve the vectors R "' , U a , as well as the values of t and q specifying the actual configuration of the system. Hence, for every a labelling a possible contact, a relationship of the form law a (t,q, U a , R "' ) = true,

called a contact law, should be available.

Contact laws

Frictionless contact

Under the present notations we are to meet the same circumstances as in Sect. 3. If the contact labelled a concerns two bodies denoted by B and B ', with common normal unit n "' directed toward B, the assumptions of no friction an no-adhesion mean :3p "' ;:::: 0 : 'R, "' = P a n "' . lt has been agreed in the foregoing to extend the definition of n "' , at least in a neighbourhood of the concerned values of t and q, to cases where g a , the normal gap, takes nonzero values and to state R "' = O if g "' > O. Define K (t ) •= {{VE E 3 1 V .n "' ;:::: O} if g a (t,q) :S 0 "' ' q • E 3 otherwise. This is the set of the values of the local right-velocity of B relatively to B' ( the latter may be a member of the system or an external obstacle with prescribed motion) which are compatible with non-interpenetration. In the first line, K a equals a half-space, hence the normal cone Nx: ,, , evaluated at the origin 0 of E 3 , equals the half-line generated in this vector space by -n "' . Otherwise K a = E 3 , so that the cone Nx:,, (0) reduces to the set {0}.

Therefore, at time t, the no-friction and no-gluing assumptions (including the case of no-contact) are equivalent to assert [START_REF] Génat | New results on Painlevé paradoxes[END_REF] The move made in Sect.3 of replacing inclusion (15) by inclusion ( 16) admits as a counterpart here the replacement of inclusion ( 24) by [START_REF] Foerster | Measurements of the collision properties of small spheres[END_REF] In fact, in this context where u is continuous, the same argument as in 2.3 shows that U a belongs to the boundary plane of the half-space K a , hence Nx; ,, (U a ) = Nx; ,, (O), while in case of no-contact Nx: ,, (U a ) = {O} whatever is U a . In short, (25) contains all the stipulations implied when a contact is declared frictionless. But, in addition, (25) entails U a E K a (t, q), since otherwise Nx: ,, ( U a ) would be empty. If g a (t, q) > 0 this actually imparts no restriction on U a while if g a (t,q) :S 0, i.e. a E J(t,q), this implies U a .n "' ;:::: O. It has been precedingly recalled that U a .n "' equals the derivative of the fonction t r-+ g a (t, q(t)). This allows one to invoke the Viability Lemma (2.1), with f = -g a , in order to prove that the assumption of (25) being verified for almost every t in I entails : non-interpenetration holds for every t > t0 , provided it holds at to.

The latter statement applies more generally to any contact law which, among other phenomenological stipulations, secures the following

• in all cases U a E K a , • if U a E interior K a , then R "' = O.
In other words, one has the implications (26) [START_REF] Frémond | Rigid body collisions[END_REF] We propose to say that a package of information, concerning the possi ble contact labelled o:, if it possesses these two properties, is a contact law of prospective type (or, in the terminology precedingly used by the author [START_REF] Moreau | Sorne basics ofunilateral dynamics[END_REF][59], a complete contact law). The underlying idea is that such a law does not properly govern the values of Ua and n a at the actual instant, but their limits on the right of this instant, assumed to exist. In fact, if U :x .n a > 0, the concerned instant is followed by a contactless time-interval. Since n a vanishes over this interval, the same is true for its right-limit.

Multicontact frictionless Dynamics

We are now to see how the observations made in Sect.3 let themselves be transposed into the present setting. Even the simple case then considered, of a particle confined by a frictionless material boundary, may exhibit multicon tact features if the boundary consists of several parts with smooth equations, each of them enforcing an inequality of the form [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF]. If two of these smooth surfaces meet to form an edge, the particle, when lying on this edge, expe riences contact forces from bath parts, the resultant of which may take any value in the convex cane generated by the respective inward normals. There fore, the writing in ( 14) remains valid provided N .i; ( t ) (q) is defined as the convex cane generated by the outward normals to the surfaces which form the edge.

When coming to general multicontact systems, one has to connect nor mality in the linear Euclidean space E 3 of the vectors of physical space with normality in the space R n of the abstract components. The following rela tionship is found to hold [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dy namics[END_REF] between the element V f a of R n and the normal unit vector n a at point M a to the contacting bodies, directed toward B [START_REF] Haug | Computer Aided Kinematics and Dynamics[END_REF] The proof of this rests on a unilateral version of the algebraic theorem of Lagrange multipliers, known in Convex Analysis as Farkas' lemma [START_REF] Rockafellar | Convex Analysis[END_REF].

In all the sequel, we shall assume that the mapping Ga is surjective of R n onto E 3 ; equivalently, its transpose c; is injective of E 3 into R n . Then À a in (28) is nonzero. Sorne special positions of a multibody system may give rise to 'wedging' effects which contradict this assumption.

In view of ( 19) and [START_REF] Curnier | A theory of friction[END_REF] this allows one to replace the laws of frictionless contact, either [START_REF] Génat | New results on Painlevé paradoxes[END_REF] or [START_REF] Foerster | Measurements of the collision properties of small spheres[END_REF], by equivalent relationships involving only the abstract components u and r a instead of Ua and n a . Under the definition (22) of J, put W(t,q) := {v E R n /Va: E J(t,q) Ôft a + v.Vfa::; O} ,

a polyhedral closed convex set. One then finds [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF] that a value r of the sum La r <> is compatible with the contact law (25) holding for allo:, if and only if -r E N w ( t ,q) ( u). Consequently, the elimination of frictionless reactions from the dynamical equation ( 21) is achieved in the writing d u F(t,q,u) -A (t,q) d t E Nw(t,q) (u). [START_REF] Heyman | The stone skeleton, structural engineering of masonry archi tecture[END_REF] This allows one to derive a time-stepping scheme for the numerical approx imation of solutions quite similar to the one presented in 3.3. The occurrence of the matrix A in (30) at the place occupied in ( 16) by the scalar factor m does not constitute an essential complication. It only means that, instead of the standard Euclidean metric of R 3 , one is using in R n the Euclidean metric defined by the positive definite matrix A.

Handling inequality conditions in terms of acceleration

The differential inclusion (30) formally resembles (5) which has been found to characterize the Sweeping Process. In fact, in the special case where F = 0 and A = 1, (30) makes the fonction t H u appear as a solution to the Sweeping Process by the moving (closed, convex) set t H W(t, q(t)). Of course, the latter is not given, since it depends on q which itself is connected to u by ( 2), but one may infer from this analogy that, similarly to the equivalence of ( 5) to (4), the inclusion (30) could be replaced by a differential equation whith right-hand member defined by a minimization property. This is the object of the forthcoming.

A time-stepping scheme of the implicit type, rests on predicting the veloc ity whithout resorting to any expression of the acceleration. In contrast, what follows is aimed at determining the right-acceleration u + , which is needed when a scheme of the explicit type is being planned, as well as an Event Driven calculation [START_REF] Abadie | Dynamic simulation of rigid bodies : modelling of frictional contact[END_REF] [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF].

Let ti denote an instant preceded by some time-interval throughout which the motion, with absolutely continuous u, satisfies [START_REF] Heyman | The stone skeleton, structural engineering of masonry archi tecture[END_REF]. Assume that u re mains continuous at instant t 1 , i.e. no collision occurs. Non-interpenetration entails that, for every a in J(t 1 , q(t 1 )), the right-derivative of the fonction t H f a(t, q( t)) at t i is :S: 0, while the left-derivative is ?: O. Hence this fonction has zero derivative of order one at t 1 .

Let us assume in this Su bsection that the fonctions fa are C 2 . The investigation of the motion by explicit time-stepping or by an ED policy rests on the assumption that t i is followed by a nonzero interval through out which (30) is verified again. It just may happen that some contacts break at ti, inducing a change of J(t,q). If u possesses a derivative on the right of t1, the fonction t H fa(t,q(t)) possesses a second derivative on the right of t1, expressed through the chain rule with a certain aa, in the form aa(t1, QI, u1) + ü"t .'V fa(t1, q1). Since the fonction is zero at t1, as well as its first derivative, non-interpenetration requires of this second derivative to be :S O. If it is strictly negative, f a ( t, q( t)) becomes strictly negative on a sub sequent interval, making the corresponding quantities r a = µ a '\7 f a vanish on this interval, as well as their limits for t .J, t 1 , assumed to exist. One thus obtains a set of complementarity conditions By joining them to the equation of dynamics OE one reduces the determination of ut and of the multipliers µ a to a Linear Complementarity Problem in standard form.

Because the matrix A is positive definite, such a problem is classically equivalent to minimizing a convex quadratic fonction in a closed convex poly hedral subset of R n . lt has been shown in [45] [START_REF] Moreau | Quadratic programming in mechanics: dynamics of one sided constraints[END_REF] that this extrema! charac terization of the acceleration may be viewed as extending Gauss' Principle of the Least Deviation to mechanical systems subject to unilateral frictionless constraints. Sorne dual minimization property is also found to characterize the contact forces.

Coulomb friction

The presence of dry friction, governed by the law of Coulomb, at the pos sible contact with label a, is expressed by a relationship of the form [START_REF] Saxcé | New inequation and functional for contact with friction[END_REF]. Traditionally, the law of Coulomb is only invoked for persistent contact, but when devising numerical schemes, making it meaningful for U a .n a � 0 and securing that a contact law of the prospective type is so stated, is just the matter of writing the code adequately.

Apart from the numerical success, the consistency of the concept of pros pective type is illustrated by the following feature.

Dropping the label a for brevity, one defines the (non adhesive) Coulomb friction at some contact point by giving the Coulomb cane C, a closed con vex conical region of E 3 to which the contact force R exerted by B' upon B should belong in any circumstance. ln the standard case, C is rotationally symmetric about the normal vector n and contains it, but more general situ ations, accounting for anisotropie friction, are possible. The law consists in a relationship between the force R and the local velocity U of B relative to B' which resembles a plasticity law in that the values of R lying in the interior of C are compatible with U = 0 only. But the 'flow rule' which characterizes the values of U compatible with a value of R lying on the boundary of C does not involve the normality of U to this boundary. In short, this is not an 'associated' flow rule ( except in the case of zero friction coefficient).

The concept of bipotential has been introduced by De Saxcé and Feng [START_REF] Saxcé | New inequation and functional for contact with friction[END_REF] as a tool for handling non-associated force/velocity (or stress/strain-rate) laws, from the theoretical standpoint, as well as in numerical techniques. Here is the aspect this concept takes in the case of Coulomb friction.

As usual, one decomposes the vectors R and U into their normal and tangential components. Let T denote the subspace of the vectors in E 3 which are orthogonal to n.

R = R T + RN n, R T E T, RN E R, U =U T +UNn, UT ET, UNE R.
Classically, if the Coulomb law is restricted to situations where RN is known, say RN = 1, the law turns out to reduce to a relationship between R T and U T of the associated type. Let D1 := {RT ET I R T + n E C} (the 'unit section' of the cone C) and define in T the real fonction (the 'dissipation fonction' of the said restricted law) TET H cp1(T) := sup{S.T I SE -Di}.

In the traditional case of isotropie friction with coefficient , , one simply has cp1(T) = , II T II-

Similarly to 5.1, define /C(t ) • = { {V E E 3 1 V.n � O} in case of contact or overlap ' q • E 3 in case of no-contact.
The Coulomb cone depends ont and q; put C = {O} in case of no-contact.

Using arguments from Convex Analysis, one establishes that the relation between the elements U and R of E 3 consisting of the system of conditions

U E /C, RE C, -U.R = cp1(UT) RN (31)
is a contact law of the prospective type which, in the standard situation, reduces to the law of Coulomb. Furthermore, one may prove \f V in /C, \f S in C : V. S + 'Pl ( V T ) SN � 0 so that [START_REF] Isac | Complementarity Problems[END_REF] expresses that the real fonction (V, S) H V. S + cp 1 (V T )SN, separately convex with regard to V and S, attains at point (U, R) its minimal value relative to the product set /C x C and that this minimal value is zero.

6 Time-stepping

Numerical dynamics of multicontact systems

As before, denote by [ti, tr], with length h, a time-step. From the approximate values qi, ui obtained for q and u at ti as the result of antecedent computation, one has to calcula te qr, ur, corresponding to tr.

The inertia matrix A(t, q) is assumed to depend smoothly ont and q, so that one chooses to approximate it throughout the time-step by the value A m it takes at the 'midpoint' t m := ti + h/2, q m := qi + hui/2. Similarly, the term F(t, q, u) in ( 21) is assumed to vary slowly enough for allowing one to replace it by the value Fm it takes at t = t m , q = q m , while the value Ui of u is used.

Thereby it is understood that the force fields or pull-back actions which, among other terms, are compounded into the fonction F, do not depend on q in too steep a way. If such is not the case, a trick to overcome this difficulty is presented in [START_REF] Jean | Frictional contact in collections of rigid or deformable bodies : numerical simulation of geomaterials[END_REF]; it results in replacing A by an artificial inertia matrix in the construction of which the possible stiff elasticity coefficients are entered.

It is also at the point ( t m , q m ) that the fonctions f o. are calculated, in order to decide which contacts are to be treated as effective. The set J( t, q(t)) is thus estimated to equal J m := J( t m , q m ) throughout the time-step.

The dynamical equation ( 21) is then discretized in the form

A m (ur -ui) = hF m + L P o. , ( 32)
aEJrn where the element p °' of R n is made of the covariant components of the impulsion at contact a, i.e. the integral P" over [ ti, tr] of the contact force n o. . Throughout the time-step, one ascribes to the linear mapping G� of ( 20) its value computed at ( t m , q m ), hence [START_REF] Jean | Dynamics in the presence of unilateral contact and dry friction[END_REF] One does the same with the linear mapping G o. and with the rheonomic term W o. of [START_REF] Cundall | A computer mode! for simulating progressive large scale movements of blocky rock systems[END_REF]. Since a discretization scheme of the implicit type is being planned, the final velocity ur is invoked, so as to calculate [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF] In the exact problem, the contact law with label a should connect U a to the contact force R °' at every instant. Our approximation scheme of the implicit type consists in connecting, through the same relationship, the final value U a .f of U a to the average of this force, namely p a / h. In usual situations such as the no-friction case described by [START_REF] Foerster | Measurements of the collision properties of small spheres[END_REF] or also the case of Coulomb friction, the contact law happens to be positively homogeneous with degree zero relatively to the contact force, so the discretized law reads law a ( t m , q m ,U a ,f, P o. )= true.

( 35)

Nonlinear Gauss-Seidel iterations

Solving the system of conditions ( 32) to ( 35) is the hard part of the compu tation. From ( 32) and ( 33), one obtains ur = u r +A;;; 1 � G�P f3 , ( 36)

(]EJm
where Ur := U i + hA;;; 1 Fm may be called the relaxed velocity. Here is an iteration technique à la Gauss-Seidel which amounts to treating a succession of single-contact problems.

Let an estimated solution ( u ;: sti , P! t J, (3 running through J m , be ob tained with [START_REF] Jean | Frictional contact in collections of rigid or deformable bodies : numerical simulation of geomaterials[END_REF] satisfied. One expects to obtain a corrected estimate, say ( uf °rr , Pfo rr ), by selecting a label a and altering only P °' , i.e. by making Pf0rr = P! ti for (3 i-a. The new estimate is astrained to verify (36); equiva lently, since the old estimate satisfies the same, [START_REF] Jean | The Non Smooth Contact Dynamics method[END_REF] and to satisfy the discretized contact law [START_REF] Jean | lmplementation of unilateral contact and dry friction in computer codes dealing with large deformations problems[END_REF]. By applying G n to both mem bers of ( 37), one gives to [START_REF] Jean | lmplementation of unilateral contact and dry friction in computer codes dealing with large deformations problems[END_REF] the form where H n := G n A;;; 1 G� is a symmetric positive definite 3 x 3 matrix. Solving [START_REF] Johansson | A linear complementarity algorithm for rigid body impact with friction[END_REF] in the unknown P::" orr is easy in some significant cases [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dy namics[END_REF], such as two-dimensional Coulomb contact or the case where Hn is axissym metric about n °' . Generally, some iterative procedures may be applied, in which the normal and tangential components of the contact force are alter natively treated as known. Anyway, when a solver has been devised, the above procedure of correcting successive estimates is iterated, with a ranging cycli cally in Jm• The decision of stopping iterations may be made on observing the magnitude of the vector P::O rr -P::'s, ti ; this actually is equivalent to check ing the precision at which each pair ( u ;: stï , P::'s, ti ) satisfies the corresponding contact law.

u corr = U esti + A -le* (P °' -P °' •) f f m a corr esti
Observe that, provided such a numerical convergence check is made, the linear operator H n in (38) may be replaced by any other mapping of E 3 into itself with zero limit at the origin, with the possible advantage of mak ing resolution easier. This replacement is also used in tricks for accelerating convergence.

The mathematical convergence of algorithms of this sort, in the case of Coulomb contact, as well as the very existence of solutions to the problems addressed, has only been established in special situations [44][40]; uniqueness cannot be expected in general.

If all the contact laws invoked are of the prospective type, one sees that, for each a, the final relative velocity Unr belongs to K n (t m , q m )-Similarly to 3.3, Remark 1, this ensures non-interpenetration with a welcome self-corrective effect if the step-length is not too large. Contact laws of the prospective type also manage automatically the possible breaking of contacts, without resorting to any analysis of complementarity conditions.

To start iterations, one needs an estimate ( uj' sti , P! t J verifying [START_REF] Jean | Frictional contact in collections of rigid or deformable bodies : numerical simulation of geomaterials[END_REF]. One may take uj'8 t i = Ur and ail P! t i equal to zero. In cases where the set Jm of the active contacts does not change much from one step to the other, it could be much more efficient to take as first guess for the P! ti the values calculated at the antecedent time-step for the contacts already present and to adjust ur ti by means of [START_REF] Jean | Frictional contact in collections of rigid or deformable bodies : numerical simulation of geomaterials[END_REF]. Anyway, iterations do not require the handling of large matrices since, in the representation of each G o. for a rigid body, only a 3 x 6 submatrix is nonzero (2 x 3 in the case of a two-dimensional multibody system).

Also observe that, due to the computation being of the implicit type relative to velocity, a single computation step is needed to check whether a given position of the system is that of a possible equilibrium. One just have to launch iterations with ui = 0 ; if Uf is found zero, the corresponding set of computed values for the contact forces is compatible with equilibrium. ln this sort of application, the inertia matrix A only provides the backbone of the computation and its value may be chosen so as to ease convergence.

N onsmooth dynamics and collisions

Discontinuous velocity functions

The preceding Sections were developed in the context of usual 'smooth' dy namics, in which the velocity fonction u is assumed locally absolutely con tinuous on the time interval J. The properly 'unilateral' feature was only the possibility offered to contacts of breaking at any instant.

In contrast, the sudden occurrence at some instant t e of a new contact, i.e. a collision is expected to generate a velocity jump. Even without such an event, the presence of dry friction at a contact point has long been known, in some cases, to forbid the existence of a smooth solution beyond some catastrophic instant, on the left of which the contact force, as well as the derivative u', may become unbounded [START_REF] Génat | New results on Painlevé paradoxes[END_REF]. This is a dynamic analogue to the locking situations familiarly met in the statics of frictional systems. Around year 1900, such a lack of smooth solution for an apparently well set problem seemed inadmissible to P. Painlevé and induced him to question the very concept of a contact force. Hence the inadequate denomination of 'Painlevé's paradox' which has been given to the observation. Today, one is accustomed to see models reaching the limits of their validity domain. Such is the case for smooth dynamics at a catastrophic instant of any sort, so that one is induced to enlarge the fonction space where u is looked for [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dy namics[END_REF] [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of Painlevé's problem[END_REF].

In order to develop Nonsmooth Dynamics, the widely accepted extension of the fonctional framework is that of R n -valued fonctions with bounded variation. More precisely, since the time interval J has not been supposed compact, we consider the linear space of the fonctions with locally bounded variation, i.e. they have bounded variation over every compact subinterval ; notation lbv(J, R n ). Any u in this space is locally Lebesgue-integrable, so that (2) remains meaningfol. The classical property of a fonction with bounded variation to possess a derivative almost everywhere has no relevance in the present situation, for such a pointwise defined derivative of u (generating values for the accelerations of the system elements in physical space) does not allow one to retrieve the fonction through integration. In fact, the possible discontinuity points of u, which make at most a countable, hence Lebesgue negligible, subset of J may be responsible in an essential way for the evolution of this fonction. For instance, a step-fonction possesses almost everywhere a derivative equal to zero without having to be a constant.

The role of the fonction u' is now played by an R n -valued measure on I, called the Stieltjes measure or differential measure of u, that we shall denote by du. In the smooth case, where u is locally absolutely continuous, one has du = u' dt, with dt denoting the Lebesgue measure (this is nothing but the differential measure of the real fonction t H t).

Classically, a fonction u in lbv(J, R n ) possesses a left-and a right-limit at every point t of I, say u-(t) and u + (t) (with the convention u-(t 0 ) = u(t 0 ) and symmetrically at the possible other end of J) and one has for every compact subinterval [a, T] of J In particular, the integral of du over the singleton {a} equals the possible jump of u at point a.

A measure-differential equation

Since one cannot rely on the second derivative of the fonction t H q(t) anymore, the accelerations of the elements of the system in physical space are also missing, as well as forces. The latter have to be replaced by E 3 - valued measures, the integrals of which over every compact subinterval of J constitute impulsions. In smooth dynamics, the impulsion measures admit, relatively to dt, some density functions whose values, at a given instant, are properly forces, so that impulsions, in accordance with the traditional definition, equal their time-integrals.

In particular, instead of a contact force R a , there will more generally be invoked a contact impulsion, an E 3 -valued measure on the time interval I. For the homogeneity of notations, let us denote by dS a this measure, without having to pay attention to any fonction s a admitting it as differential measure. The covariant components of dS ,, , relative to the parametrization in use, are R n -valued measure on I, say ds 0 , which, similarly to [START_REF] Curnier | A theory of friction[END_REF], are defined as [START_REF] Johansson | Study of frictional impact using a nonsmooth equation solver[END_REF] This expression makes sense if the linear mapping G� depends continuously ont and q, an assumption usually verified in multibody dynamics.

The differential equation of dynamics ( 21) has then to be replaced by

A(t,q)du = F(t,q,u)dt+ Lds °' , (40)
a an equality of R n -valued measures on I, called a measure-differential equa tion.

Concerning the connexion of this writing with a general formulation of classical dynamics, refer to [START_REF] Moreau | An expression of classical dynamics[END_REF]. It may be asked why ( 40) is given the re stricted form of an equality of measures, instead of an equality of distributions of arbitrary order on the interval J. Our answer, in the present context, is that nonsmoothness cornes from the contact actions which, due to unilater ality are subject to inequality requirements ( even in the more general case of adhesive contact [26][68]). This is a basic fact that distributions verifying inequalities are necessarily measures.

Such an equality of measures may equivalently be exploited as an equality of fonctions, thanks to the following observation. Given a finite collection of R-valued or R n -valued measures, such are dt, du and ds °' above, there exists (non uniquely) a nonnegative real measure on I, say dµ, relatively to which these measures respectively admit density fonctions t� E Cf oc (J, dµ; R), u� E Cf0 c (I, dµ; R n ), (s °' )� E Cfo c (J, dµ; R n ). Then ( 40) is equivalent to the equal ity of R n -valued fonctions A(t, q) u� = F(t, q, u) t� + L(s °' )�, [START_REF] Knight | Vibration-induced size separation in granular media : the convection connection[END_REF] a holding dµ-almost everywhere in J. Replacing the base measure dµ by another admissible one amounts to multiply densities by positive real fonctions [START_REF] Moreau | Bounded variation in time[END_REF].

There remains to precise how contact laws may be entered into this for malism. As observed before, in the case of Coulomb friction as well as in the no-friction case, contact laws in the form (23) are conic, i.e. positively homogeneous of degree zero with regard to the argument R/'. If the density of contact impulsion (S °' )� is introduced at the place of R °' in the contact law, the relationship so stated between the contact impulsion measure and the local velocity fonction U 0 is thus indifferent to the change of the base measure. Such a policy does not seem questionable as long as the measure dµ and, consequently, the measure dS °' are diffuse on the considered time interval. In this case u is continuous, so that the local velocity U 0 is well defined for every t.

Collisions

In contrast, assume that, at some instant t e , a velocity jump occurs as the result of a collision.

The measure du then presents an atom at the point t e of I and the same should be expected from the other measures involved in the calculation. Nec essarily the base measure dµ in use presents an atom at t e and dµ may be ad justed in such a way that this atom has unit mass, i.e. the said atom is a Dirac measure located at t e . With dµ so specified, the value P 0 = (S 0 )�(t e ) of the density of contact impulsion for the contact labelled a, is nothing else than the vector called percussion in the traditional formalization of shocks. Conse quently, the measure equation ( 40) readily applies to the standard question of determining the post-collision velocity u + (t e ), knowing the pre-collision value u -(t e ), provided some information is available about contact impulsions.

An information of this sort is implicitely included in the various pragmatic approaches proposed in the form of collision laws, since the very times of Descartes or Newton. These laws usually involve empirical parameters meant to account for the 'bounciness' of the collision ( coefficients ofrestitution) and for the possible role of friction during the process. It is today recognized that such parameters can only be identified in narrowly delimited situations [START_REF] Stoianovici | A critical study of the concepts of rigid body collision theory[END_REF] such as the collision of two, otherwise free, spherical beads. The coefficients possibly determined for a pair of free bodies are not valid anymore if any of them is subject to additional constraints [START_REF] Ivanov | On the problem of constrained collision[END_REF]. It should be kept in mind that if one of the bodies involved in a collision is part of a cluster of contacting objects, nonzero contact percussions should be expected at all contacts in the cluster.

Viewing a collision as a strictly instantaneous process is a schematic model, pertaining to the dynamics of strictly indeformable bodies. In con trast, if some slight deformability of the colliding bodies is evoked, a nonzero time-interval, say [t e , t e + 8], should be ascribed to the interaction. Large values of the contact forces are expected, the time-integral of which over this interval makes the contact percussion vector, but at this scale smooth dynamics remains applicable.

Studies of the latter sort are usually developed through a multiple scaling approach : a micro-time is introduced, in order to investigate the variation of the velocity fonction, while the variation of the position q is neglected. The direction of the sliding velocity U 0 at some of the involved contacts may be found to exhibit large variations as a fonction of the micro-time, preventing the identification of any representative value which could, phenomenologi cally, be connected with the total impulsion p a . The situation is better if U a remains zero : because the Coulomb cone, closed and convex, is constant with regard to the micro-time (see however the discussion in 8.1), the condition of n a belonging to it commutes with the integration invoked in defining P 0 • Deeper investigation of collision processes has also been conducted, by taking into account the whole deformation of the interacting bodies, either analytically [START_REF] Maw | The rebound of elastic bodies in oblique impact[END_REF] or numerically thanks to the finite element representation of deformability [START_REF] Walton | Numerical simulation of inelastic, frictional particle particle interactions[END_REF]. This also cornes to confirm that collision laws can only roughly approximate the complexity of the process.

Frictionless inelastic collision

In the time-stepping scheme of Sect. 6, contact forces were involved only through their total impulsions on each time-interval, so that it remains ap plicable to problems formulated in the nonsmooth framework. This numerical procedure handles collisions in an automatic way. It only happens that, in the course of computation, a time-step exhibits some contacts not detected at the preceding step. Since the velocity ur has to comply with some non interpenetration conditions from which ui is exempt, the contact impulsions p a involved in [START_REF] Ivanov | On the problem of constrained collision[END_REF] are expected to take values of larger magnitude than in the case of constant J m , But calculation remains the same and delivers kinematically admissible post-collision velocities. In short, because the con tact laws [START_REF] Jean | lmplementation of unilateral contact and dry friction in computer codes dealing with large deformations problems[END_REF] have been assumed positively homogeneous with degree zero in their last arguments, the algorithm treats collisions on the same footing as permanent contacts.

There remains to precise which collision model the algorithm turns out to approximate. Assume for simplicity that a single contact has corne to increment l m , i.e. there occured a single impact at some unknown instant t e of [ti, tr]. The endpoint values Uï and Uf are viewed as approximating u-(t e ) and u + (t e ) respectively.

We first have a look at the calculation to be made in the special frictionless setting of 5.2. The polyhedral set W defined in [START_REF] Heyman | The stone skeleton[END_REF], namely the set of the values of u + compatible with all non-interpenetration conditions, including the new one, is approximated by As in 6.2, one constructs the relaxed velocity Ur := Uï + hA;; -,1 F m . If the positive definite matrix A m is used to define in R n a Euclidean metric, the core of the computation simply consists in constructing ur as the nearest point ta Ur in W m• By imagining the step-length h arbitrarily close to zero, one sees that the calculation approximates the following collision law :

The post-collision velocity u + (t e ) equals, in the sense of the kinetic me tric A(t c ,q(t c )), the nearest point ta u-(t c ) in W(t c ,q(t e )).

Let us refer now to the contact locus, in the notations of 5.1. The above calculation amounts to admit that, for each a, the contact impulsion p a is connected, through the law offrictionless contact [START_REF] Foerster | Measurements of the collision properties of small spheres[END_REF], with the post-collision velocity u:;. Due to this law being of prospective type, the implication p a =/= 0 =? u:;. n a = 0 holds, meaning that, if the contact labelled a takes an effective part in the nonsmooth process, it exhibits the feature traditionally formulated by saying that the Newton restitution coefficient is zero. But, depending on the circumstances met in constructing the projection, one may also have P °' = 0, with U;t. n °' only astrained to be nonnegative.

The present model (it was introduced in [5 1 ] under the name of Standard Inelastic Shock) is somewhat more realistic than the traditional Newton state ment, in that all the contacts present at instant te are treated collectively. This produces plausible results for instance in the case, popular in the Earthquake Engineering literature, of a rectangular block rocking on a horizontal plane: if the traditional formulation, with zero restitution coefficients, was applied to all contacts, no rocking could be found. In contrast, the above collision law allows one to discuss the outcome with regard to the aspect ratio of the block. A slender block presents a succession of oscillations, progressively damped by the inelasticity of collisions, while, below some critical aspect ratio, a single episode of motion leads to permanent contact.

A three-parameter collision law

In the same line, there remains to show how non-zero friction and non-zero restitution can be entered into a computationally efficient collision model. At an instant where the system experiences a velocity jump we propose, for every contact a, to relate through the Coulomb law of friction (stated in prospective form) with friction coefficient 'Y a , the contact percussion P °' with some average value ui of the local velocity. This artificial value is defined as a weighted mean of the (known) pre-collision velocity U;;_ and of the (unknown) post-collision velocities U;t. Introducing as before the normal and tangential components of the concerned vectors, put

The empirical parameters Po: and T a: will be called the normal coefficient of restitution and the tangential coefficient of restitution at the contact labelled a, denominations justified by what follows.

Since the invoked contact law is of the prospective type, implications [START_REF] Frémond | Adhérence des solides[END_REF] and [START_REF] Frémond | Rigid body collisions[END_REF] show that P °' can be nonzero only if ui N = 0, i.e. u-;+; N = -p 0 U;: N , which formally is Newton's restitution law. But the present formulation is richer than stipulating a normal restitution law separately for each contact ; it also allows p o: = 0, in which case only the inequality U� N � 0 happens to be asserted. It is the global calculation, involving all the contacts together through the equation of dynamics, which decides between these two alterna tives.

Similarly, the global calculation, if friction is large enough, may end in the zero sliding case of Coulomb's law at contact a. Then u-:; T = -Ta u;; T , which is a law of tangential restitution.

The construction of the average velocity may readily be incorporated into the time-stepping computation scheme of Sect.6, so that the algorithm is ready to face possible collisions at every instant. At the price of a few extra lines of code at each time-step before launching the Gauss-Seidel iterations, no computation cost at all is added to the iterations themselves Of course the same collision model may be used in conjonction with other solvers [START_REF] Johansson | A linear complementarity algorithm for rigid body impact with friction[END_REF] [START_REF] Johansson | Study of frictional impact using a nonsmooth equation solver[END_REF].

The time-stepping scheme also yields plausible results when some of the frictional catastrophes referred to in 7.1 are met [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dy namics[END_REF].

It should however be admitted that relating the percussions p a to the respective average velocities U:1;, is nothing but a pragmatic trick, generating a collision law of the sort which precisely was criticized in 7.3. Before relying on it in a specific context, one has to check its practical value against any available experimental measurement.

In the simplest case of all, that of the collision of two otherwise free spherical beads, the resulting law turns out to be mathematically equivalent to another three-parameter collision model anteriorly proposed on the basis of quite different arguments [START_REF] Walton | Numerical simulation of inelastic, frictional particle particle interactions[END_REF] and which, experimentally, has been found in fairly good agreement with reality [START_REF] Foerster | Measurements of the collision properties of small spheres[END_REF].

The various circumstances met when a ball bounces against a fixed plane are also convincingly reproduced, in particular the exotic behaviour of a 'su perball' [55] [START_REF] Moreau | Numerical experiments in granular dynamics: vibration induced size segregation[END_REF].

The rocking of a slender block supported by an oscillating table has pro vided another occasion of testing the model. Sorne sensitive features of the rocking regimes have been correctly predicted by CD computation [START_REF] Raous | Experimental analysis of the rocking of a rigid block[END_REF].

It is in the domain of the numerical simulation of granular materials that the computational simplicity of the procedure proves invaluable. In fact, one may then have several ten thousands of contacts to investigate at each time step.

Sorne satisfactory comparisons of computation with experiments are pre sented in [55] [START_REF] Moreau | Numerical experiments in granular dynamics: vibration induced size segregation[END_REF], pertaining in particular to the circulatory currents observed in shaken granular samples and their connection with the segregation of grain sizes [START_REF] Knight | Vibration-induced size separation in granular media : the convection connection[END_REF].

The study of the flow of grains over a rough incline has also benefited from joining experimental investigation with CD numerical simulations [START_REF] Chevoir | Dense granular flows down an inclined plane[END_REF].

The efficiency of the method in the dynamics of granular materials mainly stems from that, once a time-discretization has been chosen, all the collisions which have been detected as occurring on a given time-step are treated to gether. This of course entails some trade-offs which have to be assessed. It does not seem to matter that the ordering of collisions, which mechanically should be successive, is only internal to the algorithm. In fact the problems in view are physically undeterministic, since a slight change in the initial con ditions is liable to generate after a short time a completely different sequence of positions and collisions. A more critical issue is that each contact point is only counted once in the analysis. An accumulation of bounces, such as those of a ping-pong ball coming to rest, is thus viewed as a single collision as soon as the successive bounces are all comprised in the time-step. This could result in underestimating the total energy loss arising from inelastic collisions. This source of error may be checked by repeating computation with step-length reduced to half: if no appreciable change is found in the energy versus time curve, one may conclude that no harm was done.

It has been stressed in the foregoing that, in compact assemblies, the bod ies involved in a collision may be part of clusters of already contacting objects, between which percussions are likely to be induced. The described procedure of global dynamical calculation does take this into account. Physically how ever, the transmission of impulses in clusters should involve elasticities in a way analogous to sound propagation, an effect which cannot be analyzed through the present model of strictly rigid objects.

Anticipating on the question of energy balance, to be studied in the next Section, one may finally rise another objection. Like other popular collision laws, the trick of average velocities offers no security against the possible violation of the laws of Thermodynamics. In fact, for colliding objects of not too usual shapes, the calculated outcome may entail some energy creation. This is a physically inacceptable conclusion, except of course for artificially boosted collisions as in electric billiard games. A systematic thermodynamic formalism has been proposed by M. Frémond [START_REF] Frémond | Rigid body collisions[END_REF], securing that the empirical collision laws one may insert in it are dissipative [START_REF] Cholet | Chocs de solides rigides[END_REF] [START_REF] Cholet | Exper imental study of collisions of angular particles[END_REF].

Energy balance and calculus for differential measures

Departing from smooth dynamics

Drawing the energy balance of a dynamical evolution is essentially relevant to the case where the external obstacles or boundaries are fixed in the (Galilean) reference frame in use. ln the smooth dynamical context of Sect.4 it was al luded to the possibility of constructing the parametrization (q) with account of linkages or permanent hinges through which members of the system would be connected with some external supports. Also these external supports have to be assumed fixed, in order to obtain a 'scleronomic', i.e. time-independent parametrization. This makes the inertia matrix A(t, q) constant with regard to its first argument and gives to the kinetic energy of the system an ex pression which, for every q, is a positive definite quadratic form in u, namely Ek = A i ju i u i /2. ln this framework of smooth dynamics, one classically es tablishes that the time-derivative of Ek equals the total power of the forces experienced by the system, including the contact forces between members and the possible action of external supports. If contacts are frictionless as well as the (fixed) possible external hinges the corresponding terms in the power vanish. In the familiar case where, additionally, the other forces in presence derive from a potential energy, the integration of derivatives yields the well known property of energy conservation.

In contrast, the elementary example of a particle hitting 'inelastically' a fixed frictionless plane shows that, in nonsmooth situations, the absence of friction does not secure energy conservation anymore. Mathematically, the reason is that the formula for the time-derivative of a product of (locally ab solutely continuous) fonctions, used in establishing the energy balance in the smooth case, has to be replaced by Jess precise relationships to be presented in 8.2 below. The physical background needs some comments.

The discussion sketched in 7.3, in which a nonzero time-interval [t e , t c +B] was, at least qualitatively, assigned to the process, makes one understand that the small deformations that the contacting bodies undergo during the collision may generate values of the local velocity vector U a disagreeing with condition U a .n a = 0, while the vector n a itself may exhibit transient varia tions. The instant power of contact forces can then be nonzero, yielding for their total work over the time-interval a nonzero, commonly negative, value. The same remark applies to the external linkages asserted to connect the system with fixed supports : even frictionless, they may transmit work to the outside world. This induces one to question, when nonsmooth circum stances are met, the legitimity of involving these linkages in the construction of the parametrization : the fact that such linkages are perfect in the sense of the traditional smooth dynamics does not secure that the percussions they impart in a collisional situation should have zero generalized components.

Nonsmooth differential calculus

Sorne rules of a calculus for l.b.v. vector fonctions of a real variable t are established in [START_REF] Moreau | Bounded variation in time[END_REF]. ln particular, one finds expressions for the differential measures of fonctions constructed from elements of l.b.v. through multilinear operations. For instance one has ( 45) and various similar formulas where left-and right-velocities are differently combined. If the scalar product in R n is denoted by a dot, this yields [START_REF] Moreau | Quadratic programming in mechanics: dynamics of one sided constraints[END_REF] Furthermore, thanks to the quadratic form u ,-+ u.u = llull 2 being nonnega tive, one establishes the following inequality, in the sense of the ordering of real measures on I [START_REF] Moreau | Problème d'évolution associé à un convexe mobile d'un espace hilbertien[END_REF] This holds true, more generally, in any Euclidean or real Hilbert linear space.

By integrating the concerned measures on the singleton {t e }, one may derive from these rules of extended differential calculus a formula for the jump of kinetic energy at the instant of a possible collision 2 applied to the Euclidean metric associated with the inertia matrix A = A(q(t c )). If ( 40) is used here to express A(u + -u-) (recall that the measure dt has zero integral on a singleton), this yields

1 1 E;t -t:; = 2 (é + u-) . A(u + -u-) = 2 (u + + u-). I>a• ( 4 9 )
a where p a = G� p a is the n-vector made of the generalized component of the percussion p a occurring at the contact labelled a. The right-hand member equivalently writes down as

� L (u + + u-) .c: p a = L �(G a u + + Gau-).P a a a " l ( _ + ) a = � 2 U a + U a . P . a (50)
Therefrom the jump of [k is expressed as a sum of terms corresponding to the respective contacts effective in physical space. Mathematically, it may be convenient to view each of these terms as defining the energy or work of the corresponding percussion. But one should keep in mind that nothing in this analysis can justify the conception that such terms could, each for its part, express a separate energy balance for some physical phenomenon occurring at the corresponding contact locus.

As an illustration, one may corne back to the frictionless collisions pre sented in 7.4. Since each percussion p a was then assumed connected with the corresponding U-;+: through the law of frictionless contact [START_REF] Foerster | Measurements of the collision properties of small spheres[END_REF], one has u-;+:. p a = O. From the right-side inequality in ( 4 7), one concludes that such collisions essentially involve some energy loss. In contrast, one might stipu late that the same law of frictionless contact should hold between p a and the arithmetic mean (U;; + U-;+; ) /2. Then (50) would entail energy preservation.

Note that if some of the contacts involved take place between a member of the system and some external obstacle, the latter must be fixed.

When invoked in the space R n of the parameters, the norm and scalar product should be understood in the sense of the kinetic metric. More gen erally than in the above investigation of an instant phenomenon, the energy balance of the possibly nonsmooth motion of the system over a time inter val may be stated as an equality of measures, in which the evolution of the matrix Ais taken into account [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dy namics[END_REF] [START_REF] Moreau | Bounded variation in time[END_REF].

Behaviour of a masonry structure

The dynamics of block assemblies has been submitted to CD computation, in order to simulate their response to seismic actions [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF]. In this domain of ap plication, it may be necessary to take also into account a slight deformability of the blocks and the presence of mortar joints [START_REF] Acary | Contribution à la modélisation mécanique et numérique des édifices maçonnés[END_REF]. N umerical simulation is applied here to a two-dimensional mode! of stone bridge in order to display its collapse under some localized forcing. The pur-pose is pedagogical, with a view to make clear that the ruin is not the result of the crushing strength of stone beeing exceded, but has to be discussed in terms of geometry and of the unilaterality of non-interpenetrability con straints [29][30]. If mortar joints are present between stones, their strength is neglected and the contact described as obeying the Coulomb law of friction. Coulomb himself, when investigating this issue, observed that, in practical instances, the friction coefficient was large enough for no sliding to occur between archstones so that the evolution toward collapse consisted of the opening of some joints through the relative rotation of the corresponding archstones about hinge points.

The structure, initially in equilibrium under its own weight, is submitted to the action of an external object with imposed motion, forcing the central archstone clown. Computation of the subsequent motion is performed twice, with the respective values of 0.6 and 0.3 for the friction coefficient at all contacts.

10 Stress distribution in a conical grain pile 10.1 Constructing a pile N umerical simulation is a popular mode of investigation in granular mechan ics. In particular, dry cohesionless granular materials are naturally modelled as collections of bodies which internet by contacts affected with Coulomb fric tion. Provided the model is sufficiently validated by comparing some measure ments, feasible in laboratory on physical experiments, with the correspond ing numerical results, computer simulation may bring precious informations about intimate mechanisms and quantities otherwise inaccessible.

The following observation has, in recent years, started a lot of speculation and controversy [11][70]. If a conical pile is created by pouring grains from a point source onto a rough rigid horizontal ground, some experimentalists have been surprised to find that the distribution of pressure at various ground contact points was not proportional to the height of the material above and even that some local minimum of ground pressure could be present at the vertical of the apex.

The distribution of stress in a granular material not only depends on the forces actually exerted on it, but strongly also on the way the granular mass has been prepared. Numerical simulation thus has to reproduce the preparation process.

Figure 3 displays a view of the following numerical experiment [START_REF] Moreau | Application des algorithmes "Contact Dynamics" aux milieux granulaires[END_REF], per formed with the Contact Dynamics algorithm of Sect. 6. A pile of about 14000 spherical grains is constructed by depositing grains, one by one, at the contact of already existing ones. Grains diameters are distributed at random, uniformly from 0.25 cm to 0.5 cm. The fixed horizontal ground roughness is simulated by a random pavement of grains with the same distribution of A large part of the deposited grains run down the pile slope before stop ping, sometimes triggering avalanches. The velocity of running grains may reach 24 cm/s. For sufficient precision in calculating these motions, including the avalanches, the step length is fixed at 2 x 10-4 s. The duration of the simu lated phenomenon equals 72 s, so that 360 000 steps have been needed. In the final stages the number of contact points is about 33 000. The Gauss-Seidel procedure remains practical on a microcomputer (but a matter of weeks) for a system of such a size because the largest part of the pile stays in quasi equilibrium : the contact impulsions calculated in the antecedent time-step for the contacts already effective are used as first guess in iterations.

Figure 4 shows, in four stages of the pile building, a slice 2 cm thick containing Oz. The grains deposited between t = 15.5 and t = 18.7 and which, at the considered stage, have their centers in this slice are represented in black. One observes that a noticable proportion of the deposited grains do not run out on the slopes, but accumulate in the central region so as to cause some plastic deformation of the existing granular mass. This deformation may be assessed by visualizing the distribution of these black grains at the successive stages. Let us call a fossile layersuch a collection of grains observed at a certain date and deposited during a specified anterior period.

Stress distribution

The programme allows one to choose a line segment in a plane drawn through Oz. By rotation about Oz, this segment generates a curved strip in the shape of a truncated cone or of a cylinder, across which force transmission will be analyzed. To this end, the contacts occurring in grain pairs with centers That a granular material, observed at large scale, admits a Cauchy stress field like any continuous medium is a familiar fact, for instance in Soi! Me chanics. When some numerical simulation has revealed the values of inter granular forces, the assessement of the average stress over a delimited sample is also a classical matter but the question of the representativity of the sam ple, in regard to its size, needs investigation. In fact, photoelastic experiments as well as numerical simulations reveal that the transmission of forces in a granular mass in equilibrium exhibits a certain inhomogeneity. Loads hap pen to be carried by privileged grains forming force chains. A sample has naturally to exceed the size of these chains in order to yield a representative value of stress. For an investigation of this subject through Contact Dynamics simulations, see [START_REF] Moreau | Numerical investigation of shear zones in granular mate rials[END_REF] [START_REF] Moreau | Évolutions en présence de liaisons unilatérales : notions de base[END_REF].

The purpose of Fig. 5 is to check that, in the present numerical experi ment, the numbers of grains is large enough for the above reciprocity property to hold reasonably true, in spi te of sampling aleas ( a two-dimensional exam ple, with irregular polygonal grains is presented in [START_REF] Moreau | Évolutions en présence de liaisons unilatérales : notions de base[END_REF]). In each point of the pile, the average density of meridian force Tn cor responding to a cylindrical strip with axis Oz has the same direction as the meridian section of a fossile layer containing this point.

This distribution of meridian forces expresses an arching effect, present even in the central region of the pile.

In the peripheral bank, which has been constructed by the deposition of avalanching grains, the fossile layers are nearly parallel to the free surface. Through the reciprocity property, this is equivalent to saying that, on a con ical eut parallel to the free surface, the average density of transmitted force is vertical. This feature looks natural since the bank is made of l ay ers suc cessively deposited by avalanches. When such a layer stops, the supporting forces it experiences from the material beneath should equilibrate its weight, and this verticality of transmitted forces is expected to persist after other l ay ers possessing the same property have been superimposed.

Ground pressure

Ground contacts: 2367 The central core of the pile clearly has part of its weight supported by the arching effect of the peripheral bank and this explains that in some experi ments, a local minimum of ground pressure has been found at the vertical of the pile apex. Figure 7 displays the results obtained in the present numerical simulation. The supporting ground surface has been divided into annular re gions with axis Oz and uniform width. For ail the contact points between the pile and the rough gound surface wich project in one of the annular regions, the total vertical component of the contact forces is divided by the area of the region, yielding the corresponding average pressure. The graph shows the respective values of this pressure (in deciPascal) plotted against the distance to Oz. Regions of small radii are naturally more subject to sampling aleas than the larger ones, but a central local minimum is in evidence. The pressure at center is smaller than the hydrostatic pressure which would be found in a liquid column of the same height and the same average volume mass, namely 3320 dPa corresponding to the height of 5.5 cm and a volume mass of 0.615 (in computation, the volume mass of the material the grains are made of has been taken equal to 1).

Anyway, the statics of a pile of dry grains has nothing to do with hy drostatics. This is evidenced by the fact that the free surface of the pile at rest is not horizontal. The merit of the reciprocity property disclosed in the foregoing is to establish a mathematical connexion between the arching effect and the inclination of the free surface.

1 -

 1 u a = �u-+ - 1 -u +aN 1 + Po: aN 1 + Po: o: N u a = --2:::_ u-+ -

( 48 ) 2 -

 482 This result, sometimes called Kelvin's theorem, actually is of purely algebraic nature and one may recognize in it the elementary identity (a+ b).(ab) = llall llbll

Fi g . 1 .

 1 Friction coefficient 0.6 Fi g . 2. Friction coefficient 0.3

Fi g . 3 .

 3 View of the pile sizes. Friction coefficient : 0.4 everywhere. Restitution coefficients : O. Grav ity : 981 cm/s 2 • The vertical of a grain center at the time of its deposition is chosen at random in the neighbourhood of the Oz axis at a maximal distance of 0.8 cm. Deposition frequency : 200 per second.

  Fi g . 4. Fossilelayers separated by the conical or cylindrical surfaces are reviewed. A contact force R is retained in the list if the corresponding contact point M has a projection onto the separating surface which belongs exactly to the strip. The meridian half-plane II containing Mis used, in order to decompose R into a component R z in the Oz direction and a component R r in the axifugal direction (the component of 'R orthogonal to II is not investigated). The sum of all R z in the list, as well as the sum of all R r , are divided by the area of the strip. This yields the components in II of a vector Tn which may be viewed as the average density of meridian force transmitted across the strip.

Fi g . 5 .

 5 Reciprocal cutsIn a classical continuous material, possessing a Cauchy stress tensor field, the similar construction would deliver an average density of meridian force related to the strip normal unit vin the form T II = O"[]V, where O"[] denotes a symmetric two-dimensional tensor. Symmetry implies a reciprocity property for a pair of (conical or cylindrical) strips drawn through a given point, with respective normal units v and v' and force densities T II and Th : one readily finds v. Th = v'. TI I, For a common value equal to zero, there cornes out that the strip with normal v is parallel to Th if and only if the strip with normal v' is parallel to T II .

Fig. 6 .

 6 Fig. 6. Reference to fossile layers

Figure 6

 6 Figure6shows the same 2 cm thick slice as before, with two fossile layers outlined. This suggests the following approximate law :In each point of the pile, the average density of meridian force Tn cor responding to a cylindrical strip with axis Oz has the same direction as the meridian section of a fossile layer containing this point.This distribution of meridian forces expresses an arching effect, present even in the central region of the pile.In the peripheral bank, which has been constructed by the deposition of avalanching grains, the fossile layers are nearly parallel to the free surface. Through the reciprocity property, this is equivalent to saying that, on a con ical eut parallel to the free surface, the average density of transmitted force is vertical. This feature looks natural since the bank is made of l ay ers suc cessively deposited by avalanches. When such a layer stops, the supporting forces it experiences from the material beneath should equilibrate its weight, and this verticality of transmitted forces is expected to persist after other l ay ers possessing the same property have been superimposed.

Fi g . 7 .

 7 Average ground pressure in annular regions