Ngoc Son Nguyen

Bernard Brogliato

Comparisons of multiple-impact laws for multibody systems: Moreau's law, binary impacts, and the LZB approach

This chapter is dedicated to the comparison of three well-known models that apply to multiple (that is, simultaneous) collisions: Moreau's law, the binary collision law, and the LZB model. First a brief recall of these three models and of their numerical implementation is done. Then an analysis based on numerical simulations, where the LZB outcome is considered as the reference outcome, is made. It is shown that Moreau's law and the binary collision model possess good prediction capabilities in some few "extreme" cases. The comparisons are made for free chains of aligned grains, and for chains impacting a wall. The elasticity coefficient, restitution coefficients, mass ratios and contact equivalent stiffnesses are used as varying parameters.

Introduction

Multiple impacts are very complex phenomena occurring frequently in multibody systems. Roughly speaking, a multiple impact occurs in a multibody system each time the system undergoes several collisions at the same time t imp . In models based on the assumption that the bodies are perfectly rigid at contact, and such that the impacts are instantaneous phenomena, the definition of t imp is clear. When deformations occur, one may consider that an impact is multiple whenever the collisions at the m contact/impact points i, which have non zero durations [t 0,i ,t f ,i ] (with t f ,i = +∞ for some models -think of an overdamped linear spring-dashpot [1, Section 2.1] Ngoc Son Nguyen Institut GeM (UMR CNRS 6183), IUT de Saint-Nazaire, Département Génie-Civil, 58 rue Michel Ange, BP 420, 44606 Saint-Nazaire Cedex, France, e-mail: Ngocson.Nguyen@ univ-nantes.fr Bernard Brogliato INRIA Grenoble-Rhône Alpes, Univ. Grenoble-Alpes, Laboratoire Jean Kuntzman, 655 avenue de l'Europe, 38334 Saint-Ismier, France, e-mail: bernard.brogliato@inria.fr 1 [START_REF] Schwager | Coefficient of restitution for viscoelastic spheres[END_REF]), overlap and consequently may influence each other due to dynamic couplings between the various contact points. One subtlety in the definition of a multiple impact, is that some previously active contacts with zero local relative velocity, may participate into it. This is the case for the two well-known classical systems: chains of aligned balls (like the Newton's cradle) where several balls are in contact before the shock, or the planar rocking block that rotates around one corner. In both cases, one is obliged to take the previously lasting contacts into account even if the multiple impact is triggered at a single contact. In a Lagrange dynamics framework with generalized coordinates q, impacts are associated with unilateral constraints, which are defined from p gap functions f i (q) (signed distances) that define an admissible domain Φ for the generalized position, i.e., q(t) ∈ Φ for all t ≥ 0. Impacts correspond to trajectories hitting the boundary of Φ (denoted bd(Φ)) with a non zero normal velocity, i.e. ∇ T f i (q(t)) q(t -) < 0 if f i (q(t)) = 0. In most cases bd(Φ) consists of co-dimension p ≤ p submanifolds {q ∈ C | f i (q) = 0, for some 1 ≤ i ≤ p}, of the configuration space C Φ. When a co-dimension p boundary submanifold is attained with p ≥ 2 (a kind of singularity of bd(Φ) where two smooth hypersurfaces intersect), one speaks of a p -impact. For instance, the 2-dimensional rocking block with concave base and two corners, undergoes a 2-impact during a classical rocking motion. Consider a chain of n aligned spheres, where one sphere at one edge of the chain hits the other n -1 ones which are at rest and in contact with no pre-constraint: this is an n -1-impact.

Just as for single impacts, several classes of contact/impact models can be used in multiple impacts [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]:

• (i) Algebraic models that relate post and pre-impact velocities as q(t + ) = F ( q(t -))

for some function F , which may be explicitly or implicitly defined. • (ii) First-order dynamics following the Darboux-Keller approach [1, section 4.3.5]: positions are assumed constant, the impact force impulse is used as the new time scale. • (iii) Second-order dynamics that use rheological compliant models with lumped flexibility like spring-and-dashpot linear (Kelvin-Voigt, Maxwell, Zener) or nonlinear models (Kuwabara-Kono, Simon-Hunt-Crossley, etc), discrete-element methods (DEM), or finite-element methods (FEM).

All models have some advantages and drawbacks. It is not our objective in this chapter to classify or to rank models. Rather, we consider three well-known models that belong to classes (i) and (ii), and we compare them in terms of their velocity outcomes, on the benchmark of chains of aligned balls. The results therefore complete [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]Chapter 6] which is restricted to chains of three aligned balls. Our results also indicate when Moreau's and the binary laws may provide realistic outcomes. Since multiple impacts in chains of balls are essentially determined by the nonlinear waves that travel through the chain, we pay attention to characterize, when possible, the waves associated with the domains of applicability of these two impact laws.

Remark 1. In this work we restrict ourselves to frictionless constraints. Remark 2. Multiple impacts are therefore intrinsically different from infinite sequences of single impacts with an accumulation, like in the bouncing ball system. However some approaches for multiple impacts may yield some kind of infinite sequence of impacts, sometimes instantaneously (this may occur for instance in the binary collision model, or with the so-called Han-Gilmore algorithm [1, section 6.1.2] which is not always guaranteed to converge in a finite number of steps, or to converge to a unique solution [3, section 3.4]). This is closely related to another feature of multiple impacts, that is the possible discontinuity of trajectories with respect to the initial data [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF][START_REF] Paoli | Continuous dependence on data for vibro-impact problems[END_REF].

System's dynamics

In this chapter we mainly deal with chains of n aligned balls (or more generally aligned grains not necessarily spherical) with radii R i > 0 whose dynamics is as follows:

           M q(t) = Λ(t) f i (q) = q i+1 -q i -(R i+1 + R i ) ≥ 0, 1 ≤ i ≤ n -1 M = diag(m i ), 1 ≤ i ≤ n, (1)
and Λ(t) ∈ IR n is the vector of generalized contact forces between the balls, q = (q 1 , q 2 , . . . , q n ) T . The gap functions f i (q) are signed distances between adjacent balls and represent the unilateral constraints in the chain. We have Λ = ∇f(q)λ, with λ ∈ IR n-1 is the vector of Lagrange multipliers associated with the unilateral constraints. We obtain the next equalities that will be useful later:

   ∇ T f i+1 (q)M -1 ∇ f i (q) = -m -1 i+1 , ∇ T f i-1 (q)M -1 ∇ f i (q) = -m -1 i ∇ T f i-2 (q)M -1 ∇ f i (q) = 0, ∇ T f i (q)M -1 ∇ f i (q) = m -1 i + m -1 i+1 . (2)
In terms of the kinetic angles θ i j between the submanifold (or hypersurfaces) defined by the equalities f i (q) = 0 and f j (q) = 0, we obtain (see [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF]Equation (6.66)] for the definition of a kinetic angle) when all masses are equal to m > 0:

θ i,i+2 = π 2 , θ i,i+1 = π 6 . (3)
Roughly speaking, and without going into further considerations other than this preliminary geometrical analysis, this means that monodisperse chains of aligned balls may have complex dynamics at impacts because they may not satisfy the conditions that guarantee continuity of trajectories with respect to initial data [START_REF] Paoli | Continuous dependence on data for vibro-impact problems[END_REF]. As shown in [3, Appendix A], the 3-ball chain is equivalent to a particle in the plane hitting in an angle, whose dynamics may be quite complex [START_REF] Towne | One-dimensional collisions and Chebyschev polynomials[END_REF]. As is well-known there is another "natural" set of coordinates for the chain, using conservation of linear momentum. Let z i = f i (q) for each 1 ≤ i ≤ n -1, and z 0 = ∑ n i=1 m i q i . Then z0 = 0 (by adding the n lines of the dynamics which just translate Newton's law of action/reaction). We have z = Nq + L for some easily obtained N ∈ IR n×n and L = (0, R 2 + R 1 , . . . , R n + R n-1 ) T . The n × n mass matrix becomes in the z coordinates N -T MN -1 =     1 0 . . . 0 0 0 M 0     , with M = MT ∈ IR (n-1)×(n-1) positive definite.

Let z = (z 1 , . . . , z n-1 ) T , the dynamics in [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF] then becomes in a reduced form:

   Mz = λ z i ≥ 0, 1 ≤ i ≤ n -1. (4)
If all the balls are in contact at the impact time, then z i (0) = 0. Though the dynamics in (4) looks simpler than (1), this is not necessarily the case because M may not be diagonal.

The multiple-impact models

In this section the three models: Moreau's impact law, the binary collision model and the LZB approach are described, and some of their features are analyzed.

Moreau's impact law

Moreau's impact law belongs to class (i). It is primarily formulated as an extension of Newton's kinematic restitution law, in a Lagrange dynamics framework, and with a global restitution coefficient. Since it can also be expressed in local frames at the contact points, as a linear complementarity problem with unknown the local velocities, it is convenient to implement in event-capturing time-stepping schemes. As such this is the law that is implemented in the software packages SICONOS 1 and LMGC90 2 . It was introduced in [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF].

Let us describe it now. We consider a Lagrangian system with generalized coordinates q ∈ IR n , symmetric positive definite mass matrix M(q) ∈ IR n×n , and a set of unilateral constraints f i (q) ≥ 0, 1 ≤ i ≤ m, defined from the differentiable gap functions f i :

IR n → IR, such that ∇ f i (q) ∆ = ∂ f i
∂ q (q) T = 0 for all q such that f i (q) = 0 (it is assumed that gradients do not vanish on the boundary of the admissible domain). The non negative multipliers associated with the unilateral constraints are denoted λ i , and they are supposed to satisfy complementarity 1 http://siconos.gforge.inria.fr/4.1.0/html/index.html 2 https://git-xen.lmgc.univ-montp2.fr/lmgc90/lmgc90_user/wikis/ home conditions f i (q)λ i = 0. In a compact form one obtains 0 ≤ λ ⊥ f(q) ≥ 0, with λ = (λ 1 , . . . , λ m ) T , f(q) = ( f 1 (q), . . . , f m (q)) T . The right-hand side of the smooth dynamics (outside impacts) is equal to Λ ∆ = ∇f(q)λ with 0 ≤ λ ⊥ f(q) ≥ 0, which under some suitable assumptions and using nonsmooth analysis can be rewritten equivalently as Λ(t) ∈ -N Φ (q(t)), the normal cone being generated by the gradients at the active constraints f i (q) = 0 (the set of active constraints is denoted as I (q) in the sequel).

Remark 3. Readers who are not familiar with convex and nonsmooth analysis, should simply think of normal and tangent cones as a generalisation of normal and tangential subspaces, with normal cones being generated by the gradients of the active constraints on the admissible domain boundary. As we explain next, using the normal and tangential cones is very useful to understand particular features of Moreau's impact law, because they provide a clear geometrical picture of the collision process, a point of view that is lost if these tools are not used.

Moreau goes a step further and replaces the normal cone to the admissible domain N Φ (q), by the normal cone to the tangent cone V (q) = {v ∈ IR n |v T ∇ f i (q) ≥ 0, for all i ∈ I (q)}, computed at the right-limit of the velocity, i.e. the following inclusion is proposed: Λ(t) ∈ -N V (q(t)) ( q(t + )), whose right-hand side we choose to name Moreau's set [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF]. We have to assume that V (q) is non empty, which may be guaranteed by suitable constraint qualification. We also assume that the pre-impact velocity satisfies q(t -) ∈ -V (q(t)). When no constraints are active, i.e. I (q) = / 0, then one sets V (q) = IR n . In this case N V (q) (•) = {0} as expected (contact forces vanish).

In a more general setting, Moreau's set is computed at w(t)

∆ = q(t + )+e q(t -) 1+e
, where e is a global coefficient of restitution (CoR) (global in the sense that it applies to all the contacts), i.e.: Λ(t) ∈ -N V (q(t)) (w(t)) ⊆ -N Φ (q(t)) 3 . One important consequence of using Moreau's set is that since V (q) ⊆ IR n is a convex polyhedral set for velocities (while Φ may be in general non convex and non polyhedral), the calculations of the normal cone are doable as we show next. When an impact occurs at time t, Λ t = ∇f(q(t))λ t is the contact force impulse and the system's dynamics becomes:

M(q(t))( q(t + ) -q(t -)) = ∇f(q(t))λ t ∈ -N V (q(t)) (w(t)). (5)
The objective of the above developments may appear obscure to many readers, however as we show next they pave the way towards a sound and practical impact law. First of all one may use a basic result of convex analysis which states that for a symmetric positive definite matrix M, two vectors x and y, and a closed non empty convex set

K , M(x -y) ∈ -N K (x) ⇔ x = proj M [K ; y]
, where proj M denotes the orthogonal projection in the metric defined by M, i.e.: x = argmin z∈K 1 2 (xz) T M(xz). Using this and after few manipulations we obtain from (5):

M(q(t))( q(t + ) -q(t -)) ∈ -N V (q(t)) (w(t)) q(t + ) = -e q(t -) + (1 + e)proj M(q(t)) [V (q(t)); q(t -)], (6)
where we used that multiplying both sides of (5) by 1 1+e > 0 does not change the right-hand side which is a cone. Other equivalent formulations exist [1, Equations (5.60) (5.61)]. Using now a corollary of the celebrated Moreau's two cones Lemma [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF]Equation (B.18)], it follows that ( 6) is equivalent to:

q(t + ) = q(t -) -(1 + e)proj M(q(t)) [N Φ (q(t)); q(t -)], (7)
where under some constraint qualification (like the so-called Mangasarian-Fromovitz CQ) we can state that N Φ (q) is the polar cone to V (q) (the admissible domain Φ needs not be convex for this). Moreau's law is a global (generalized) law which gives in one compact form the post-impact velocity. The question is then how to compute the projection in a way that is convenient for numerical implementation. Since the projection is done in the metric defined by M(q), we can define the (outwards) normal cone as

N Φ (q) = {w ∈ IR n |w = -∑ i∈I (q) λ i n q,i , λ i ≥ 0}, with n q,i = M(q) -1 ∇ f i (q) √ ∇ T f i (q)M(q) -1 ∇ f i (q)
, the (inwards) normal vector to the submanifold defined by f i (q) = 0 in the kinetic metric. It is however not trivial to calculate the projection onto a cone in the general case of the kinetic metric. We may start directly from the impact dynamics in [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques[END_REF] to get a more tractable expression. Indeed Moreau's set can be expressed as N V (q(t)) (w) = {z ∈ IR n |z = -∑ i∈K (w) λ i ∇g i (w), λ i ≥ 0}, with: g i (w) = w T ∇ f i (q), K (w) = { j ∈ I (q)|g j (w) = 0} ⊆ I (q). Thus, I (q) collects indices of active position constraints, while K (w) collects indices from active velocity constraints inside position active constraints: we see at once that Moreau's set implies a two-stage process: first look at positions, second look at velocities. In a more mathematical language there is a lexicographical inequality imposed at the contact local kinematics. Notice that we can rewrite equivalently N V (q(t)) (w) = {z ∈ IR n |z = -∑ i∈I (q) λ i ∇g i (w), 0 ≤ λ i ⊥ g i (w) ≥ 0}, and we have

∇g i (w) = ∂ g i ∂ w (w) T = ∇ f i (q) = ∂ f i ∂ q (q) T .
Then we obtain:

   M(q(t))( q(t + ) -q(t -)) = ∑ i∈I (q) λ i ∇ f i (q) 0 ≤ λ i ⊥ g i (w) = ∇ f i (q) T w ≥ 0. (8)
In this approach the multiplier λ i has to be interpreted as the contact force impulse at time t, i.e. λ i = λ t,i . Let I (q) = {i 1 , . . . , i l }, and denote f I (q) (q) ∆ = ( f i 1 (q), f i 2 (q), . . . , f i l (q)) T , so that ∇f I (q) (q) = (∇ f i 1 (q), . . . , ∇ f i l (q)) ∈ IR l×n . In the same way we denote λ t,I (q) = (λ t,i 1 , . . . , λ t,i l ) T , and U n,I (q) = (U n,i 1 , . . . ,U n,i l ) T , with U n,i ∆ = ∇f i (q) T q the normal local velocity at contact i. From (8) and using the expression of w(t), we obtain: U n,I (q) (t + ) -U n,I (q) (t -) = D I (q) (q) λ t,I (q) 0 ≤ λ t,I (q) ⊥ U n,I (q) (t + ) + E nn U n,I (q) (t -) ≥ 0 D I (q) (q) = ∇f I (q) (q) T M(q) -1 ∇f I (q) (q), [START_REF] Giouvanidis | Modeling contact in rocking structures with a nonsmooth dynamics approach[END_REF] with E nn = diag(e). This form of Moreau's impact law for normal local velocities is very interesting because it takes the form of a Mixed Linear Complementarity Problem (MLCP), that is numerically tractable. See [1, Lemma 5.2, Corollary 5.1] for existence and uniqueness of solutions to this MLCP. It may be seen as a generalized Newton's impact law, however it is worth noting that it is not a mere application of Newton's law at each active contact. Indeed there is a complementarity condition and inertial couplings through the Delassus' matrix D I (q) (q) ∈ IR l×l . We see from (9) that Moreau's law is kinetically consistent (non negative impulse). If e ∈ [0, 1] it is also energetically consistent [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF]Equation (5.61)] and it can be shown to be kinematically consistent as well (admissible post-impact velocities) using [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques[END_REF]. Indeed we obtain: q(t + ) = proj M(q(t)) [V (q(t)); q(t -)] + e{-q(t -) + proj M(q(t)) [V (q(t)); q(t -)]} (10) The three terms of the right-hand side belong to V (q(t)) and since e ≥ 0 the postimpact velocity belongs also to the tangent cone (a convex cone being closed under addition) and is thus admissible.

Actually, though it is not particularly useful from the calculation point of view, the expression in Equation (6) (or in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]) is valuable to visualize how Moreau's law works from simple geometrical arguments in the plane, as illustrated in Figure 1. This figure demonstrates that Moreau's law outcome is strongly influenced by the (kinetic) angle between the constraints (denoted α on the figure). This is the reason why it can posssess good predictability in the case of multiple impacts where waves play negligible role, but the system's geometry is crucial. For instance, planar rocking blocks follow this intuitive rule: slender blocks have a kinetic angle ≥ π 2 and are likely to rock more easily than flat (or stocky) blocks which have a kinetic angle ≤ π 2 [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF]Remark 6.10]. This is confirmed in [START_REF] Giouvanidis | Modeling contact in rocking structures with a nonsmooth dynamics approach[END_REF][START_REF] Giouvanidis | Nonsmooth dynamics analysis of sticking impacts in rocking structures[END_REF] where tangential effects are added to secure no sliding of the block. On the contrary in chains of balls, the wave propagation is a crucial mechanical effect which is mainly ruled by contact flexibilities. A purely kinematic impact law which does not contain any information on contact stiffnesses will in most cases fail to predict correctly the outcome. It may however in some very particular cases provide good results as shown in Section 5. Remark 4. Other kinematic impact laws have been proposed and studied in the literature [START_REF] Caselli | Collision of three balls on a plane[END_REF][START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibody-multicontact systems, with application to the planar rocking block and chains of balls[END_REF][START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagragian mechanics with impacts and friction[END_REF][START_REF] Gharib | Shock absorption using linear particle chains with multiple impacts[END_REF][START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF][START_REF] Najafabadi | Impacts in multibody systems: modeling and experiments[END_REF][START_REF] Rodriguez | Study of Newton's cradle using a new discrete approach[END_REF][START_REF] Winandy | A maximal monotone impact law for the 3-ball Newton's cradle[END_REF], or using a Poisson coefficient and a two-stage linear complementarity problem [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF]. It would be worth studying them along the q(t -) ∈ -V (q(t))

V (q(t))

q(t) N Φ (q(t)) q(t) q2 (t -) q1 (t -) q1 (t + ) q2 (t + ) = 0 q3 (t -) q3 (t + ) V (q(t)) N Φ (q(t)) Φ Φ α α
Constraints angle α ≤ π 2 : if e = 0 then q(t + ) = 0 Constraints angle α ≥ π 2 : post-impact velocities when e = 0 Fig. 1 Moreau's law and constraints angle (planar case).

same lines as done in Section 5. This is left as a future work. Notice however that as shown in [START_REF] Glocker | An introduction to impacts[END_REF], Poisson-Pfeiffer-Glocker and Moreau's law are equivalent when a unique global CoR is used, though in general Poisson's hypothesis yields multiple impact laws with larger post-impact velocity set than Moreau's one [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]Chapter 3]. Finally, Moreau's law may be in some cases formulated as a quadratic problem under non convex constraints [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]Proposition 3.4] where the cost function represents the energy dispersion. Most of the above results are taken from [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF], an alternate proof of (9) for Moreau's law can be found in [START_REF] Glocker | An introduction to impacts[END_REF]Proposition 5.6]. See also [START_REF] Aeberhard | Theoretical and experimental treatment of perfect multi-contact-collision[END_REF] for a geometric analysis of multiple impacts and a characterization of the domain of admissible post-impact velocities.

Let us come back to [START_REF] Giouvanidis | Modeling contact in rocking structures with a nonsmooth dynamics approach[END_REF]. It implies the following LCP: 0 ≤ λ t,I (q) ⊥ D I (q) (q) λ t,I (q) + (I l + E nn )U n,I (q) (t -) ≥ 0 [START_REF] Caselli | Collision of three balls on a plane[END_REF] If the active constraints are functionally independent (⇒ l ≤ n), then D I (q) (q) is positive definite and this LCP always has a unique solution. Let us calculate it for the chain of balls, using [START_REF] Schwager | Coefficient of restitution for viscoelastic spheres[END_REF], where we assume that during the shock all the balls are in contact, hence l = m = n -1:

D I (q) =                 (m -1 1 + m -1 2 ) -m -1 2 0 0 . . . . . . 0 -m -1 2 (m -1 2 + m -1 3 ) -m -1 3 0 . . . 0 0 -m -1 3 (m -1 3 + m -1 4 ) m -1 4 0 0 . . . 0 . . . . . . . . . -m -1 n-3 0 0 0 . . . 0 -m -1 n-3 (m -1 n-3 + m -1 n-2 ) -m -1 n-2 0 0 . . . . . . 0 -m -1 n-2 (m -1 n-2 + m -1 n-1 ) -m -1 n-1 0 0 . . . . . . 0 -m -1 n-1 (m -1 n-1 + m -1 n )                 (12)
We have the following for an impact occurring at t = 0: Proposition 1. Consider a chain of n aligned balls in [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF]. Let m i = m > 0 for all 1 ≤ i ≤ n. Let also e = 0, and the pre-impact conditions are chosen as q1 (0 -) = 1m/s and qi (0 -) = 0 m/s for 2 ≤ i ≤ n (hence U n,I (q) (0 -) = (-1, 0, . . . , 0) T ). Then the unique solution of ( 11) is

λ 0,I (q) = m n+1        n n -1 n -2 . . . 1       
, which yields U n,I (q) (0 + ) = (0, . . . , 0) T .

Proof: In this case

D I (q) = 1 m            2 -1 0 . . . 0 -1 2 -1 0 . . . 0 0 -1 2 -1 0 0 . . . . . . 0 . . . 0 0 . . . 0 -1 2 -1 0 . . . 0 -1 2            , (13)
which is positive definite as Lemma 1 shows. The result follows by inspection, since there is a unique solution to the LCP.

Remark 5 (Dependent active coordinates). In case the active constraints are dependent, then D I (q) 0 and since it is a symmetric matrix, D I (q) (λ 1 t,I (q)λ 2 t,I (q) ) = 0 for any solutions λ 1 t,I (q) and λ 2 t,I (q) of the impact LCP. Therefore U n,I (q) (t + ) is uniquely defined from the first line in [START_REF] Giouvanidis | Modeling contact in rocking structures with a nonsmooth dynamics approach[END_REF] (see [1, Lemma 5.2, Corollary 5.1] for the same analysis in a slightly more general framework).

Now we have the next result:

Lemma 1. The Delassus' matrix in [START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagragian mechanics with impacts and friction[END_REF] has full rank and is therefore positive definite.

Proof: Consider a matrix as in [START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagragian mechanics with impacts and friction[END_REF] with dimension n × n, and denote it as D n . It is not difficult to show that det(D n ) = 2det(D n-1 )det(D n-2 ), for all n ≥ 3, and letting D 1 = 2. It follows that provided det(D n-1 ) = n and det(D n-2 ) = n -1, then det(D n ) = n + 1. One checks that this is true for n = 3, since det(D 2 ) = 3 and det(D 1 ) = 2. Hence this is true for all n ≥ 3. Due to the fact that the Delassus' matrix is at least positive semi definite, the result follows.

Proposition 1 shows that Moreau's law creates some distance effect with non zero impulse at all contacts, and that all balls are stuck together after the shock (maximal dispersion of the kinetic energy in accordance with [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]Proposition 3.4]). Notice however that λ 0,I (q) > 0 (component-wise) implies from (11) that λ 0,I (q) = -D -1 I (q) (I l + E nn )U n,I (q) (0 -), so that -E nn U n,I (q) (0 -) = U n,I (q) (0 + ). In our case U n,1 (0 -) = -1 m/s so this implies that U n,1 (0 + ) = e m/s: this is true for e = 0 in the above conditions. Calculations for the 3-ball system show that this is also the case when e = 1 [1, p.271].

Proposition 2. Consider the chain of n aligned balls in [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF] with m i = m > 0. Suppose that λ 0,I (q) > 0 (each contact undergoes an impact with positive impulse), with pre-impact relative velocity U n (0 -) = (-1, 0, . . . , 0) T (so that I (q) = {1, . . . , n -1}). Then U n (0 + ) = (e, 0, . . . , 0) T .

Let us now state the following result. We still assume that q1 (0 -) = 1 m/s, and qi (0

-) = 0 m/s, 2 ≤ i ≤ n. Proposition 3. Let e = 1, m i = m > 0, 1 ≤ i ≤ n, U n (0 -) = (-1, 0, . . . , 0) T . As- sume that q1 (0 + ) = 2-n n m/s, qi (0 + ) = 2 n m/s for 2 ≤ i ≤ n (so that U n (0 + ) = (1, 0, . . . , 0) T ).
Then the kinetic energy is conserved, and λ 0,I (q) = D -1 [START_REF] Caselli | Collision of three balls on a plane[END_REF].

I (q) [U n (0 + )- U n (0 -)] > 0 (component-wise) is the solution of the LCP in
Proof: Preservation of the kinetic energy follows from a simple calculation. Notice that U n (0 + ) -U n (0 -) = (2, 0, . . . , 0) T , consequently we need to know only the first column of D -1 I (q) , where D I (q) has the structure shown in the proof of Proposition 1. Let us denote matrices with this structure, and of dimension p, as D p . In fact it can be shown by induction that the first column of D -1 p is equal to 1 det(D p ) (det(D p-1 ), det(D p-2 ), . . . , 2, 1) T , where det(D p ) = p + 1. Therefore since we have

I (q) = {1, . . . , n -1}, the first column of D -1 I (q) is equal to 1 det(D I (q) ) (det(D n-2 ), det(D n-3 ), . . . , 2, 1) T > 0.
Therefore λ 0,I (q) is twice this vector and is positive. We have

D I (q) (q) λ 0,I (q) + (I l +E nn )U n,I (q) (0 -) = U n,I (q) (0 + )-U n,I (q) (0 -)+(1+e)U n,I (q)
(0 -) = 0 which ends the proof since the impact LCP has a unique solution.

It is also checked that the linear momentum of the chain is preserved. Therefore under the stated assumption, Moreau's impact law is unable to separate the balls 2 to n, while ball 1 "rebounds" on the chain and imparts a non zero velocity to the n -1 other balls. It has limited predictability in terms of energy dispersion (see Moreau's line in [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]Figure 2.6] for the 3-ball system). This is what has motivated researchers to extend it while remaning in a rigid-body approach, and this is what motivates us to analyse which are the cases when it does predict correctly the post-impact velocity in section 5. Remark 6. Solving the impact LCP in [START_REF] Caselli | Collision of three balls on a plane[END_REF] allows one to compute the projection in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], i.e. the index set J (q). We could start from the reduced dynamics (4) in which the calculations for the tangent and normal cones are simplified since the constraints define the first orthant. However projections are made in the metric defined by M which is no longer a diagonal matrix. In these coordinates the Delassus' matrix is D I (q) = M and U n,i = żi . Thus there is nothing special to gain using (4) instead of (1).

The binary impact model

Contrary to Moreau's law which handles all impacts at the same time, the binary impact model handles impacts separately. To do so, the multiple impact problem is assumed to be a succession of binary collisions between rigid particles, so collisions are independent of each other. Each binary collision between two balls can be completely solved by using the conservation law of momentum and Newton's kinematic restitution law:

m 1 q- 1 + m 2 q- 2 = m 1 q+ 1 + m 2 q+ 2 , q+ 2 -q+ 1 = -e n ( q- 2 -q- 1 ), (14)
where superscripts (-) and (+) indicate the pre-and post-impact velocities, and e n is the normal restitution coefficient which takes a value from 0 for purely dissipative collision to 1 for purely elastic collision. Note that each binary collision is assumed to be central: the collision occurs only in the normal direction of the contact as illustrated in Figure 2.

q- 1 q- 2 q+ 1 q+ 2 After collision Before collision m 1 m 2 m 2 m 1
Fig. 2 Two particles before and after a binary collision.

The post-impact velocity of each ball is obtained by solving the system of linear equations ( 14):

     q+ 1 = m 1 -m 2 e n m 1 + m 2 q- 1 + (1 + e n )m 2 m 1 + m 2 q- 2 , q+ 2 = (1 + e n )m 1 m 1 + m 2 q- 1 + m 2 -e n m 1 m 1 + m 2 q- 2 . (15)
In the case where the two balls have the same mass and the first ball comes to collide the last one at rest, the post-impact velocities are:

     q+ 1 =
1e n 2 q-1 ,

q+ 2 = (1 + e n ) 2 q- 1 . (16)
If the collision is purely elastic (e n = 1), the first ball stops and the last one moves forward after collision with a velocity equal to the pre-impact velocity of the first ball. This means that the energy and momentum of the first ball are entirely transfered to the last one.

While the outcome of a binary collision is easily obtained, the definition of the succession of binary collisions is not straightforward. One can try to mimic the wave propagation induced by the shock in a granular media to define the sequence of binary collisions. Let us consider a granular monodisperse chain composed of n elastic identical beads as an example. The beads are numbered 1, 2,..., n from the left to the right. When the first ball collides the other balls which are at rest, a solitary wave is initiated and propagates from the left to the right. According to the wave propagation, the succession of binary collisions can be defined as follows: ball 1 collides ball 2, then ball 2 collides ball 3, ..., then ball i collides ball i + 1, ..., and at the end ball n -1 collides ball n. Applying the rule [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF] from the first to the last binary collision, we obtain the impact outcome as follows: balls 1 to n -1 stop and ball n moves forward with a velocity equal to the pre-impact velocity of ball 1. This sequence of binary collisions is also true for a tapered chain where the bead diameter decreases progressively and it has been used by several authors to study the momentum and energy propagation in tapered chains [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF][START_REF] Harbola | Pulse propagation in tapered granular chains: An analytic study[END_REF][START_REF] Machado | Momentum and energy propagation in tapered granular chains[END_REF][START_REF] Rosas | Pulse propagation in granular chains: The binary collision approximation[END_REF]. It is worth mentioning that for elastic monodisperse chains or tapered chains and for the considered particular initial condition, i.e. the first ball collides the other balls at rest, the sequence of binary collisions is uniquely defined. However, this is not true for most of cases. Let us demonstrate this point by considering a monodisperse chain of 10 dissipative beads. We apply the binary collision rule [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF] with the restitution coefficient e n = 0.5 to the sequence of binary collisions defined above. The preimpact velocity of the first bead is equal to 1 and the velocity of each bead after this sequence of binary collisions is shown in Table 1. It can be seen that beads enter into collisions again after the first sequence of binary collisions: there are indeed potential collisions between balls 1 and 2, between balls 2 and 3, and so on. Even for an elastic chain, we can encounter this problem. Let us take an elastic decorated chain (Figure 3) as an example. For this granular chain, three small balls of mass 0.5m are placed periodically between four big balls of mass m. Table 2 shows that there are several potential collisions between balls after the first sequence of binary collisions. A question that arises here is which order of binary collisions should we consider when there are several binary collisions to be handled. We present here two strategies that can be used for a granular chain for which collisions start at the left end and propagate to the right end.

1) Binary collisions are always handled from the left to the right. This means that among the set of possible collisions, the collision at the contact with the least value of the index k is handled first.

2) The order of binary collisions is unimportant so binary collisions can be randomly handled. This strategy has been adopted in [START_REF] Crassous | Impact of a projectile on a granular medium described by a collision model[END_REF][START_REF] Valance | Granular medium impacted by a projectile: experiment and model[END_REF]. Table 2 Bead velocity for a decorated chain after a sequence of binary collisions from the left to the right. It is worth mentioning that the selection of binary collisions with the left-to-right order or the random order presented above is not physically justified. Let us apply these two strategies to a disordered chain of 30 elastic balls. For this kind of granular chain, ball masses are randomly distributed. Sequences of binary collisions are randomly selected with the uniform distribution law. Figure 4 shows a comparison between the impact outcomes obtained with the two considered strategies. One can see that the impact outcome depends strongly on the chosen sequence of binary collisions. In addition, different random sequences of binary collisions lead to different impact outcomes as shown in Figure 5. This is intimately related to the fact that the trajectories are in general discontinuous with respect to initial data as we pointed out in Section 3.1.

Another issue of the binary collision model is that the sequence of binary collisions can tend to infinity before the impact process ends, even for simple cases. For example, Towne and Hadlock [START_REF] Towne | One-dimensional collisions and Chebyschev polynomials[END_REF] have found out analytically that the number of binary collisions for a chain of three balls is infinite if the number z defined in [START_REF] Rodriguez | Study of Newton's cradle using a new discrete approach[END_REF] satisfies z ≥ 1 (see [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF] for more discussions):

√ e n + 1 √ e n 1 1 + m 2 m 1 1 + m 2 m 3 . ( 17
)
The number of binary collisions increases quickly with the number of balls, in particular for dissipative chains (e n < 1) as shown in Table 3 for a disordered chain with e n = 0.9. One can see that the binary collision model is not able to determine the impact outcome with 70 balls for the left-to-right order and with 50 balls for the random order because the number of binary collisions to be handled is too big. In summary, the binary collision model presents three main drawbacks:

1) The impact outcome is possibly not unique, which is related to the discontinuity with respect to the initial data; 2) The impact outcome depends on the chosen order of sequence of binary collisions; 3) The number of binary collisions to be handled is possibly infinite.

The LZB model

This way of treating multiple impacts has been introduced in [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF][START_REF] Liu | Energy dissipation and dispersion effects in a granular media[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems: Part II. Numerical algorithm and simulation results[END_REF] and we briefly summarize it in this section. It has been validated through extensive comparisons between experimental and numerical data in [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF][START_REF] Liu | Impact/contact dynamics in a disc-ball system[END_REF][START_REF] Nguyen | Shock dynamics in granular chains: numerical simulations and comparisons with experimental tests[END_REF][START_REF] Wang | Nonsmooth dynamics of a 3D rigid body on a vibrating plate[END_REF][START_REF] Zhang | Dynamics of planar rocking-blocks with Coulomb friction and unilateral constraints: comparisons between experimental and numerical data[END_REF][START_REF] Zhao | Planar dynamics of a rigid body system with frictional impacts. II. Qualitative analysis and numerical simulations[END_REF] for chains of balls, rocking blocks, bouncing dimers and other setups. This is a model of the class (ii), based on the Darboux-Keller approach [1, section 4.3.5]. As such it is based on the next fundamental assumptions:

1. Forces other than impact forces are negligible during the collision process. 2. Positions are constant during the collision process. 3. Tangential stiffnesses are infinite. 4. The impact consists of a compression phase followed by an expansion phase.

0 compression expansion δ i λ i λ M i M i δ M i δ r i
Fig. 6 The bi-stiffness force/indentation model for the LZB model at contact i.

Then the impact dynamics consist of a first-order dynamics whose state is the velocity, and the time-scale is the impact force impulse. Though the Darboux-Keller shock dynamics have a long history for two-body single impacts, it is only recently that its extension to multiple impacts has been proposed with the use of energetic coefficients of restitution (CoRs) [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems: Part II. Numerical algorithm and simulation results[END_REF]. We summarize the LZB dynamics now, when applied to chains of aligned balls. Let us start from (1): M q(t) = Wλ(t), where W ∆ = ∇f(q) is constant. In this example M and W are constant so the constant position assumption is useless. During the impact we will denote the infinitesimal impulse as dP ∆ = λ, so that the so-called Darboux-Keller dynamics writes Md q = WdP ⇔ M d q dP = W, after a time rescaling has been performed. The next basic assumption is that at each contact i, one has the force/indentation relation λ i = K i (δ i ) η i , where K i is the contact equivalent stiffness and η i is the elasticity coefficient (η i = 1 for linear elasticity, η i = 3 2 for Hertz' elasticity). More precisely, the LZB model may be designed with a mono-stiffness compression/expansion model, or a bi-stiffness compression/expansion model [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF], or even a tri-stiffness model [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]Figure 4.4]. Let us describe the bi-stiffness model as in Figure 6. During the compression phase (from the origin to M i ) one has λ c,i = K i (δ i ) η i , during the expansion (or restitution) phase one has λ e,i = λ M i

δ i -δ r i δ M i -δ r i η i (see [1, section 4.2.1.2]
for a short history about bi-stiffness models). The dashed area corresponds to the dissipated energy during the shock, δ M i is the maximal identation, δ r i is the residual identation. The work done by the contact force during the compression

phase is W c,i = δ M i 0 λ i (δ i )dδ i = 1 1+η i K i (δ M i ) η i +1
, and during the restitution phase W e,i =

δ r i δ M i λ i (δ i )dδ i = -1 1+η i K i (δ M i ) η i (δ M i -δ r i ).
Actually the bi-stiffness model is a piecewise-continuous model which states that

λ i = K i (δ i ) η i ) if δi ≥ 0 (compres- sion), and λ i = K e i (δ i -δ r i ) η i if δi < 0 (expansion), where K e i = K i δ M i δ M i -δ r i η i . Cal- culations show that the energetic CoR at contact i satisfies e 2 i, * = - W e,i W c,i = 1 - δ r i δ M i = K i K e i 1 η i , hence δ r i = δ M i (1 -e 2 i, *
). Perfectly plastic impacts with e i, * = 0 imply that δ r i = δ M i , so that the expansion phase has zero duration and the point (δ M i , 0) is reached instantaneoulsy from the maximum compression point M i .

1 3 4 R M r1 r2 R' M R' M R A 2 R
Fig. 7 The potential energy when the contact point is located at expansion phase. The next step is to calculate the contact force as a function of the potential energy. Starting from λ i = K i (δ i ) η i , and using dλ i dt = λ i dλ i dP i , one finds that:

λ i (P i (t)) = (η i + 1) P i (t) 0 K 1 η i i ∇ f T i q dP i η i η i +1 . (18)
Further calculations not recalled here allow one to show that even in case of precompression (with λ i (0) = 0) one has

λ i (P i (t)) = (1 + η i ) η i η i +1 K 1 η i +1 i (E i (P i (t)))) η i η i +1 , (19)
where E i (P i ) is the potential energy at contact i, i.e. E i (P i ) = E 0,i + P i (t) 0 δi (P i )dP i , where E 0,i is the potential energy due to pre-compression. Taking pre-compression into account is crucial because such multiple impacts usually involve repeated impacts at the same contact point, which correspond to an impact starting again while the zero indentation has not been reached yet, see Figure 8. Repeated impacts render the problem more complex. A crucial result is [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF]Theorem 3.1]. Let us consider Figure 7. Then [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part I. Theoretical framework[END_REF]Theorem 3.1] guarantees that a compression-expansion cycle OR δ r1 (curves 1 and 3) is equivalent, from the energetic point of view, to a cycle ARδ r2 (curves 2 and 4), where the compression would finish ar R (respectively at R). This allows to prove the following. When the contact point i moves from M i to R i along the expansion curve in Figure 7, the recovered energy is δ R,i δ M,i λ i (δ i )dδ i , and we obtain at contact i4 :

δ R i δ M i λ i (δ i )dδ i = δ r2,i δ M i λ i (δ i )dδ i - δ r2,i δ R i λ i (δ i )dδ i = -e 2 i, * δ M i δ 0 λ i (δ i )dδ i - δ r2,i δ R i λ i (δ i )dδ i (20)
where e i, * is the energetic CoR at contact i. According to Stronge [START_REF] Stronge | Impact mechanics[END_REF], the energetic CoR e i, * is defined as e 2 i, * = -W e i /W c i where W c i and W e i are the respective works done by the contact force during the compression and expansion phases. The term premultiplied by -e 2 i, * is equal to the area enclosed by the curve δ r2 Rδ R in Figure 7. Let us assume that the force/identation relationship remains the same for the second compression/expansion phase (i.e., the elasticity properties do not vary). Using this and after manipulations it follows that the potential energy along the repeated impact in Figure 8 is given as follows, where Q denotes a generic point along the force/indentation curve:

E(P(t)) =                        E 0 + P(t)) 0 δ (P(s))dP(s) Q ∈ OM 1 E M 1 + 1 e 2 * P(t) P M 1 δ (P(s))dP(s) Q ∈ M 1 R E R + P(t) P R δ (P(s))dP(s) Q ∈ RM 2 E M 2 + 1 e 2 * P(t) P M 2 δ (P(s))dP(s) Q ∈ M 2 B, (21)
where E M 1 is the residual potential energy at point M 1 , and so on. As a next step, one can use [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF] to derive the distributing law between infinitesimal impulses dP i and dP j at contact points i and j, respectively:

dP i dP j = (1 + η i ) η i 1+η i (1 + η j ) η j 1+η j K 1 1+η i i K 1 1+η j j (E i (P i )) η i 1+η i (E j (P j )) η j 1+η j . ( 22
)
It is noteworthy that if all contacts have the same elasticity coefficient, the distributing law simplifies and shows that the ratio in [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF] depends only on the stiffnesses

ratio K i K j 1 1+η
. It is a well-known fact that the post-impact velocities in chains of aligned balls indeed do not depend on the absolute values of the equivalent contact stiffnesses, but only on their ratio, in case of linear elasticity (see for instance [1, Section 6.1.3]). This result generalizes it. In summary, the potential energy can be calculated along [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF], while the infinitesimal impulse ratio is given by [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF]. Now, contrarily to the case of a single collision where one can make a time-scale change, passing from time t to the impact force impulse dP (since the contact/impact forces are always assumed to be non-negative, and positive for times strictly inside the collision interval), one has dP > 0 and this time rescaling is valid. In case of multiple impacts, one has to choose a so-called primary contact where it is guaranteed that the impulse does not become constant, for otherwise the time rescaling becomes impossible with this impulse. Thus one chooses the primary impulse as the impulse from contact i where the potential energy at this contact E i (P i ) is maximal amongst the various contact points. We obtain the multiple impact Darboux-Keller equations:

1. (contact parameters): K j , η j , e j, * , 1 ≤ j ≤ n -1.

(dynamical equations):

Mdq = WdP, (23)
where the impulse increment dP j at a contact j is related to the impulse increment dP i at another contact i by the distributing law [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF]. The impulse increment dP j can be also related to the time increment dt by the relation:

dP j = λ j dt = (1 + η i ) η i η i +1 K 1 η i +1 i (E i (t)) η i η i +1 dt, (24)
with the contact force λ j computed with equation ( 18). 3. (potential energy (bi-stiffness model)):

E j (P j ) = E Tra, j + 1 Tra P j (t) P Tra (t) ∇ f j (q) T q dP j ( 25
)
where Tra = 1 if δ j > 0 (compression), Tra = e 2 j, * if δ j < 0 (expansion), E Tra, j is the accumulated potential energy at the beginning of the integration, and P Tra (t) depends on the impulse value at the beginning of the subphase (see ( 21)). 4. (impact termination): E j (P j ) = 0 and δ j ≤ 0 at all contacts 1 ≤ j ≤ n -1.

We have δ j = ∇ f j (q) T q = qi+1qi . Notice that (25) could be rewritten in its differential form

dE j = 1 Tra ∇ f j (q) T q dP j ⇔ dE j dP i = 1 Tra δ j (P)Γ ji (E i (P i ), E j (P j )). (26)
with Γ ji = dP j /dP i and initial condition E j (P Tra, j ) = E Tra, j . The multiple impact Darboux-Keller equations is therefore a set of first-order nonlinear and coupled piecewise smooth differential equations, with states q, Γ , E, and state-dependent switching conditions at times of maximum compressions ( δ j = 0, points M 1 , M 2 in Figure 8) or repeated impacts (point R in Figure 8).

Remark 7. 1. The bi-stffness model has several drawbacks: it does not model a bounded maximal contact force, it is a rough representation of plasticity (if plastification is the primary source of dissipated energy), it models dissipation during the expansion phase (while dissipation could occur also during the expansion phase). However it can be improved as described in [3, Section 4.2.4]. 2. The LZB approach can also be formulated with Coulomb friction at contacts [START_REF] Zhao | Planar dynamics of a rigid body system with frictional impacts. II. Qualitative analysis and numerical simulations[END_REF]. 3. The CoRs e i, * can be estimated off-line from pairwise collisions between balls i and i + 1.

4. We employ the word "balls", however the chain may consist of other types of elementary particles than spherical balls, like beads or polyhedral grains. 5. We have written δ j (P) because due to dynamical couplings stemming from M and W in ( 23), the local velocity may depend on several contact impulses. 6. As we shall see in Section 4.3, it is possible to dispense with the distributing law in [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF] which is quite time-consuming during numerical integration (see [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]Chapter 4] for a complete exposition of the event-driven algorithm for the LZB model, in particular the algorithm for the primary impulse selection). The distributing law is nevertheless quite interesting since it highlights in which way the different contacts interact one with each other. 7. We see that the LZB model allows to include the effects of contact flexibilities (which are crucial in chains of balls impacts) while disregrading positions variations. This is done thanks to the distributing law.

Numerical resolution

The numerical algorithms which are used to compute the post-impact velocities, may differ from one impact law to the next. Let us describe now how the above three models of multiple impacts are treated numerically.

Moreau's impact law

As alluded to above the great advantage of Moreau's law is that it is naturally embedded into the discrete-time version of Moreau's sweeping process for Lagrangian systems, using a suitable event-capturing scheme that stems from Moreau's catchingup algorithm. The numerical aspects of the sweeping process applied to mechanical systems, are treated in detail in [START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF]. Let us briefly introduce the catching-up algorithm. We start from the second order sweeping process:

M(q)dv + F(q, v,t)dt ∈ -N V (q) (w), (27)
where v = q dt-almost everywhere, dv is the so-called differential measure associated with the acceleration (which cannot be a function at impact times since the velocity has a discontinuity) so that ( 27) is a measure differential inclusion (MDI).

Outside impacts we have dv = q(t)dt. At an impact time t the MDI ( 27) is equivalent to [START_REF] Towne | One-dimensional collisions and Chebyschev polynomials[END_REF], that is dv = ( q(t + )q(t -))δ t , with δ t the Dirac measure at t. The basic time-stepping method for ( 27) is as follows on [t k ,t k+1 ), with constant time-step

h = t k+1 -t k > 0, k ≥ 0:        M(q k )(v k+1 -v k ) + h F(q k , v k ,t k ) ∈ -N V (q k ) v k+1 + ev k 1 + e q k+1 = q k + hv k+1 , (28)
We can proceed as we did in section 3.1 to transform [START_REF] Towne | One-dimensional collisions and Chebyschev polynomials[END_REF]. We denote

F k ∆ = F(q k , v k ,t k ). v k+1 -v k + h M(q k ) -1 F k 1 + e ∈ -M(q k ) -1 N V (q k ) v k+1 + ev k 1 + e ⇔ v k+1 + ev k 1 + e + h M(q k ) -1 F k 1 + e - ev k 1 + e - v k 1 + e ∈ -M(q k ) -1 N V (q k ) v k+1 + ev k 1 + e ⇔ v k+1 = -ev k + (1 + e) proj M(q k ) [V (q k ); - h M(q k ) -1 F k 1 + e + v k ∆ =v k ] ⇔ v k+1 = -ev k + (1 + e) argmin z∈V (q k ) 1 2 (z -vk ) T M(q k )(z -vk ) (29) Notice that if f i (q) > 0 for all 1 ≤ i ≤ n -1, then v k+1 = v k -h M(q k ) -1 F k .
We infer that the next velocity can be computed by solving a quadratic problem under conic varying constraints. A next step is to compute this projection using complementarity. To this aim we notice first that N

V (q k ) (w k+1 ) = {z ∈ IR n |z = ∑ i∈I (q k ) -λ i ∇ f i (q k ), 0 ≤ λ i ⊥ w T k+1 ∇ f i (q k ) ≥ 0}.
Thus we obtain:

M(q k )(v k+1 -v k ) + h F k = ∇f I (q k ) (q k ) λ I (q k ),k+1 ⇔ v k+1 -v k + hM(q k ) -1 F k = M(q k ) -1 ∇f I (q k ) (q k ) λ I (q k ),k+1 ⇔ ∇f I (q k ) (q k ) T (v k+1 -v k + hM(q k ) -1 F k ) = D I (q k ) (q k ) λ I (q k ),k+1 (30)
where D I (q k ) (q k ) = ∇f I (q k ) (q k ) T M(q k ) -1 ∇f I (q k ) (q k ) is the Delassus'matrix of (position) active constraints at step k. Denoting the local velocities as U n,I (q k ),k we obtain the mixed LCP:

U n,I (q k ),k+1 -U n,I (q k ),k + h∇f I (q k ) (q k ) T M(q k ) -1 F k = D I (q k ) (q k ) λ I (q k ),k+1 0 ≤ λ I (q k ),k+1 ⊥ U n,I (q k ),k+1 + e U n,I (q k ),k ≥ 0 (31)
where we used the expression for w T k+1 ∇ f i (q k ) in the complementarity conditions. The similarity between (31) and ( 9) is obvious. Once the set of active constraints has been computed, one can solve the mixed LCP [START_REF] Liu | Energy dissipation and dispersion effects in a granular media[END_REF] to compute U n,I (q k ),k+1 and λ I (q k ),k+1 . Once λ I (q k ),k+1 is known, one can use the first line in [START_REF] Liu | Energy dissipation and dispersion effects in a granular media[END_REF] to obtain v k+1 and then q k+1 . There exist quite efficient algorithms to solve mixed LCPs, some of which are implemented in the INRIA SICONOS software package.

Remark 8. The rationale behind the above is that the elements inside the normal cone -N V (q k ) v k+1 +ev k 1+e are an approximation of ∇f I (q) (q) λ t,I (q) ([t k ,t k+1 ]), that is the measure of the interval [t k ,t k+1 ] by ∇f I (q) (q) λ t,I (q) . Thus even at an impact time this is a bounded quantity (in fact, the impact magnitude).

In practice the event-capturing method in (28) can be modified to cope with energy conservation, accuracy, etc [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF][START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF]. An important feature is that it is shown to converge [START_REF] Dzonou | A convergence result for a vibro impact problem with a general inertia operator[END_REF], hence for small time-steps the numerical solutions must be close to the analytical ones.

Binary collision model

The binary collision model is solved in an iterative manner until no binary collision is found. For a chain of balls where the impact starts at the left end, the balls are numbered 1, 2, ..., n, and the contacts are numbered 1, 2, ..., s from the left to the right. In this case, we can handle the left-to-right sequence of binary collisions proposed in Section 3.2 by using Algorithm 1. This algorithm can also be used to handle a random sequence of binary collisions by selecting randomly a binary collision in set I instead of getting the minimum value in set I.

LZB impact model

The LZB impact model presented in Section 3.3 can be integrated with respect to the impulse scale. To do so, the contact, at which the potential energy is maximum, is chosen as the primary contact for each integration step. The impulse increment dP j at each contact is related to the one at the primary contact by the distributing law [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF]. Two singularities may be encountered during the integration. The first singularity may occur at the beginning of the impact process where the potential energy is zero at all contacts. The second one may occur during the impact process when a contact, which has left the impact process previously, enters again into the impact process. When a singularity occurs, the distributing law ( 22) must be regularized. The interested reader can refer to [3, section 4.2.8] for the regularization techniques and for the integration algorithm. It is worth mentioning that this integration technique requires a significant computational effort to select the primary contact among all contacts and to handle the singularities at each integration step. In addition, when the primary contact changes from one contact to another, the impulse increment dP j at each contact computed with the distributing law [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF] changes brutally, which might slow down the convergence of the algorithm. Algorithm 1 Handling binary collisions according to the left-to-right order in a granular chain.

Require: qi , m i for all particles i = 1, 2, ..., n Require: e j for all contacts j = 1, 2, ..., s Ensure: q+ i for all particles i = 1, 2, ..., n // Initialize for i = 1 → n do q+ i ← qi end for IsTermination ← f alse N ← 0 number of binary collisions handled //Iterations while IsTermination = f alse do while impact is not yet terminated IsTermination ← true // Find all binary collisions to be handled I ← ∅ set of all binary collisions to be handled

for j = 1 → s do if q+ j+1 -q+ j < 0 then Add j to I q- j ← q+ j q- j+1 ← q+ j+1
IsTermination ← f alse end if end for // Select a binary collision in set I and handle it k ← min(I) get minimum value in I

q+ k ← q- k m k -m k+1 e k m k + m k+1 + q- k+1 (1 + e k )m k+1 m k + m k+1 q+ k+1 ← q- k (1 + e k )m k m k + m k+1 + q- k+1 m k+1 -e k m k m k + m k+1 N ← N + 1 end while
The LZB model can be also integrated with respect to the time scale. To do so, the Darboux-Keller equation ( 23) is first discretized using the Euler explicit method:

qk+1 = qk + M -1 W∆ P k , ( 32
)
where k is an integration step (k = 1, 2, ..., N). The impulse increment dP k j at each contact is obtained by integrating [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF] with the Euler explicit scheme:

∆ P k j = t k+1 t k λ j (t)dt ≈ λ k j ∆t = (1 + η j ) η j 1+η j K 1 1+η j j (E k j ) η j η j +1 ∆t. (33)
A singularity occurs with (33) when a contact enters into the impact process at an integration step k, i.e. E k j = 0. In this case, ∆ P k j can be approximated as:

∆ P k j = t k+1 t k λ j (t)dt ≈ 1 2 λ k+1 j ∆t = 1 2 K j (δ k+1 j ) η j ∆t ≈ 1 2 K j ( δ k j ∆t) η j ∆t. ( 34
)
The potential energy is computed by discretizing ( 25):

E k+1 j = E k j + δ k j + δ k+1 j 2 ∆ P k j , if δ k+1 j ≥ 0, ( 35
)
E k+1 j = E k j + 1 e 2 j, * δ k j + δ k+1 j 2 ∆ P k j , if δ k+1 j < 0. ( 36
)
The impact process can be considered to be terminated at a step k if:

E k j = 0, and δ k j ≤ 0, ∀ j = 1, 2, .., s. (37)
The interested reader can follow Algorithms 2, 3 and 4 to implement the resolution of the LZB model with respect to time into a programming language.

Algorithm 2 Integration up to the end of the impact process.

Require: q0 , M, W, ∆t Require: E 0 j : initial potential energy at all contacts j = 1, 2, ..., s Require: K j , η j , e j, * for all j = 1, 2, ..., s Ensure: q, P at the end of the impact process //Initialize δ0 ← -W T q0 P 0 ← 0 //Integration t ← 0 Time scale IsTermination ← f alse //IsTermination = true: impact is over //IsTermination = f alse: otherwise k ← 0 while IsTermination = f alse do while the multiple impacts not yet terminated Check status of each contact and the termination condition with Algorithm 3 Integrate up to the end of the current step with Algorithm 4 t ← t + ∆t //Advance to the next step k ← k + 1 end while For a comparison between the two above integration algorithms, we consider a monodisperse chain of 1000 elastic beads where the first bead with a velocity of 1 m/s collides the other beads at rest. The CoR e * is then equal to 1.0 for all contacts and the Hertz's contact law (η = 3/2) is used for each contact. The other parameters are: Young's modulus E = 203 GPa, Poisson's coefficient ν = 0.3, ball radius r = 0.01 m and mass density ρ = 7780 kg/m 3 . For this chain, the postimpact velocities of balls must satisfy the energy conservation. The integration with respect to impulse with a step size ∆ P = 10 -6 N.s needs about 2.2 × 10 7 steps and consumes about 380 s CPU time. The resulting post-impact velocities of balls satisfy the energy conservation with a relative error of about 1.5 ×10 -5 . With regard Algorithm 3 Check status of each contact and the termination condition at the beginning of a step k.

Require: δ k j , E k j for all j = 1, 2, .., s Ensure: f lag k j for all j = 1, 2, .., s 1: // f lag k j = 0: contact does not come into collision 2: // f lag k j = 1: contact begins the compression phase 3: // f lag k j = 2: contact is already in the impact process Ensure: IsTermination 4: IsTermination ← true 5:

for j = 1 → s do 6: if E k j = 0 then 7: if δ k j ≤ 0 then 8: f lag k j ← 0 9: δ k j > 0 10: f lag k j ← 1 11:
IsTermination ← f alse 12:

end if 13: else E k j > 0 14:

f lag k j ← 2 15:
IsTermination ← f alse 16:

end if 17: end for to the integration with respect to time using [START_REF] Nguyen | Shock dynamics in granular chains: numerical simulations and comparisons with experimental tests[END_REF], a step size ∆t = 10 -8 s results in about 2.9 × 10 6 steps and about 38 s CPU time. The resulting post-impact velocities of balls satisfy the energy conservation with a relative error of about 2.0 × 10 -7 . The difference between the solutions obtained with the two integration algorithms is about 0.03%. It can be concluded that the integration algorithm with respect to time is about ten time faster than the integration algorithm with respect to impulse for the considered chain. The first one would take more advantages for systems with higher number of particles. Require: M, W, ∆t, qk Require: η j , K j , e j, * , E k j , δ k j , P k j for all j = 1, 2, ..., s Ensure: qk+1 , δ k+1 j , E k+1 j , P k+1 j for all j = 1, 2, ..., s 1: //Compute the impulse increment at each contact ∆ P k j 2: for j = 1 → s do 3:

Comparisons

Free chains of aligned beads

if f lag k j = 0 then Contact does not come into the collision process 4:

δ P k j ← 0 5:

else if f lag k j = 1 then Contact begins the collision process 6:

∆ P k j ← 1 2 K j ( δ k j ∆t) η j ∆t 7:
else if f lag k j = 2 then Contact has been already in the collision process 8:

∆ P k j ← (1 + η j ) η j 1+η j K 1 1+η j j (E k j ) η j 1+η

Varying the elasticity coefficient η

To study the effect of the elasticity coefficient η in the LZB model, a monodisperse chain composed of 100 elastic beads (e * = 1) is considered. The stiffness K i is the same for all contacts, while the elasticity coefficient η is varied. The impact outcomes given by the LZB model for different values of η are shown in Figure 9 and are compared to the impact outcomes given by Moreau's law and by the binary collision law in Figure 10. It is worth mentioning that the impact outcomes given by Moreau's law and by the binary collision model are independent of the elasticity coefficient η. It can be seen that the elasticity coefficient η affects greatly the impact outcome given by the LZB model. For a very small value of η (η = 10 -3 , for example), only the first ball bounces back and the remaining balls move forward with almost the same velocity after impact. This is like the first ball impacting the other balls which are rigidly bonded. It is interesting to note that this particular impact outcome is given by Moreau's law (Figure 10). As η increases, fewer balls move forward after impact as shown in Figure 11. For a high enough value of η (η = 3, for example), only the last ball moves forward after impact with a velocity almost equal to the velocity of the first ball before impact, and the other balls are almost at rest, which is the outcome given by the binary collision model (Figure 10). This is like the first ball impacting the other balls which are separated each other by a gap. Let us take the outcome given by the LZB model as a reference outcome and then quantify the gap between the impact outcome v + obtained with Moreau's law or the binary collision law and the impact outcome v + lzb obtained with the LZB model by the following gap measure:

Φ = v + -v + lzb v + lzb × 100%, ( 38
)
where . is the Frobenius norm of a vector. Figure 12 shows the gap measure Φ defined for outcomes given by Moreau's law and by the binary collision model versus the elasticity coefficient η. It can be seen that, by varying the elasticity coefficient η from a very small value to a big value, the impact outcome given by the LZB model, initially close to the outcome given by Moreau's law, gets away from the latter one but gets closer to the outcome given by the binary collision model. Except for extreme values of η, the impact outcomes obtained with Moreau's law and with the binary collision model are quite far from that given by the LZB model. For spherical homogeneous beads, Hertz's contact model (η = 3/2) is widely adopted in the literature. In this case, the binary collision model gives an approximation of the impact outcome with an error of about 17% compared to the LZB model, while Moreau's law gives an unrealistic outcome.

It is interesting to note in Figure 13 that although the outcomes given by the binary collision model and by Moreau's law are poor for small and big values of η (η < 0.01 and η > 1.0), respectively, they can be considered to be good in term of the number of balls for which these models give a good post-impact velocity compared to the one given by the LZB model. The post-impact velocity v + i of a ball given by Moreau's law or the binary collision model is considered to be good compared to the result v + i,lzb obtained with the LZB model if: where ε is a precision which is chosen to be equal to 0.05 in this study. Let us analyze the wave propagation in the considered granular chain when varying the elasticity coefficient η in the LZB model, and the link between the wave propagation and the impact outcome. Figure 14 shows the potential energy E versus time t at the first 20 contacts (from left to right) for different values of η. It can be seen that the wave propagation is greatly affected by the elasticity coefficient η. Three classes can be observed: (i) strongly localized wave at the first contact for very small values of η (Figures 14.a and 14.b), (ii) attenuated and dispersed wave for intermediate values of η (Figures 14.c, 14.d and 14.e) and (iii) dispersion-free wave for big values of η (Figure 14.f). Herrmann et al. [START_REF] Herrmann | How does the ball chain work?[END_REF] also observed the dispersionfree for the elasticity coefficient η = 3.0. The wave propagation results from the compliance of solid bodies and it is an important dynamical effect which should be taken into account in an impact model. The LZB model takes into account the compliance effect at contacts between particles by using the contact model shown in Figure 6. As a result, it is capable of reproducing the wave propagation induced by a shock and then the impact outcome. On the other hand, Moreau's law and the binary collision model neglect completely this compliance effect and make use of two opposite assumptions: the first one assumes that all collisions occur simultaneously, while the second one assumes that collisions occur in a sequential manner. The first assumption can be justified for the wave propagation category (i) so the impact outcome given by the LZB model coincides with the one given by Moreau's law (Figures 9.a, 9.b and 10). The sequential collisions are observed for the wave propagation category (iii); as a consequence, the impact outcome given by the LZB model coincides with the one given by the binary collision model. The category (ii) corresponds to the wave dispersion for which the shock initiated at the first contact spreads out spatially and the energy induced by impact is shared by many particles. A measure for this dispersion effect in term of post-impact kinetic energies of balls was introduced in [3]:

| v + i -v + i,lzb | v + lzb < ε, (39)
E (J) t (s) t (s) E (J) (a) η = 0.001 (b) η = 0.01 E (J) t (s) E (J) t (s) (c) η = 0.4 (d) η = 1.0 E (J) t (s) E (J) t (s) (
C KE = 1 T + 1 N N ∑ i=1 (T + i -T + ) 2 , (40)
where T + i is the post-impact kinetic energy of ball i (T + i = m i ( q+ i ) 2 /2), and T + is the mean post-impact kinetic energy:

T + = 1 N N ∑ i=1 T + i . (41)
According to (40), the higher the value of C KE is, the lower the dispersion effect is. For a chain of n balls, C KE reaches the maximum value of √ n -1 for the case where the energy after impact is concentrated at one ball, and the other balls are at rest. This chain exhibits zero dispersion effect, also called dispersion-free [START_REF] Herrmann | Simple explanation of a well known collision experiment[END_REF][START_REF] Reinsch | Dispersion-free linear chains[END_REF]. Figure 15 shows the dispersion measure C KE obtained with the LZB model versus the elasticity coefficient η. The maximum value of C KE for the considered chain of 100 balls is √ 100 = 10. It can be seen that the dispersion effect is very weak for very small values of η, and it increases as η increases until η ≈ 0.4 where C KE reaches its minimum value. This means that the dispersion effect is maximum for η ≈ 0.4 for which a strongly dispersed wave propagation can be seen in Figure 14.c. When η increases beyond 0.4, the dispersion effect decreases and almost vanishes for η = 3.0. It is worth mentioning that the dispersion-free obtained for η = 3.0 corresponds to the sequential wave propagation shown in Figure 14.f. It was shown in [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF] that, for a chain of 3 balls, the dispersion measure C KE increases monotonically with η. However, this monotonic dependency of C KE on η no longer exists for a chain with a high number of balls. Figure 15 also shows that Moreau's law and the binary collision model give good impact outcomes for extreme values of η for which the dispersion effect is very weak.

Varying the contact equivalent stiffnesses

It can be seen in the distributing law ( 22) that the impact outcome does not depend on the value of the contact stiffness K i if the latter and the elasticity coefficient η i are the same for all contacts. In this section, we show how the difference in stiffness between contacts affect the impact outcome and for which cases the outcome of the LZB model coincides with the ones given by Moreau's law and the binary collision model. For this study, we set the elasticity coefficient η = 3/2 for all contacts in the monodisperse elastic chain considered in Section 5.1.1 and vary the stiffness K i at each contact according to the following linear law: K i = K i-1 + αK * , with a coefficient α and a reference stiffness K * . If α > 0, the contact stiffness increases progressively from the left to the right of the chain, and the reference stiffness K * is set to the first contact. Otherwise, the contact stiffness decreases progressively, and the reference stiffness K * is set to the last contact. It should be noted that the value of the reference stiffness K * is not of importance.

Figure 16 shows the impact outcome given by the LZB model for different values of α. It can be seen that the impact outcome changes slightly when the contact stiffness is progressively decreased (α < 0) and approaches the one given by the binary collision model. Despite a very strong decrease in contact stiffness (α = -10 4 ), we cannot reach closely the latter one: the two first balls still bounce back after impact. This means that the dispersion-free outcome cannot be reached for the considered granular chain if only contact stiffnesses are varied. In fact, Reinsch [START_REF] Reinsch | Dispersion-free linear chains[END_REF] has developed an analytical analysis for a granular chain with the linear contact model (η = 1) and has shown that the dispersion-free outcome can only be reached if the mass of each ball and the stiffness of each contact are both varied according to some specific laws. On the other hand, the impact outcome changes greatly when the contact stiffness is progressively increased (α > 0) and gets closer to the one given by Moreau's law. The latter one is closely reached for a very strong increase in contact stiffness (α = 10 4 ). Figure 17 shows the gap measure defined in (38) for Moreau's law and the binary collision model versus coefficient α. It is clear that the outcomes given by these two impact laws can be approached by progressively increasing and decreasing the contact stiffness, respectively.

The link between the impact outcome and the wave propagation during impact can be clearly observed in Figure 18. It can be seen that the solitary wave, which travels in a Hertzian monodisperse chain (Figure 18.c), is not significantly disturbed by a progressive decrease in contact stiffness. On the other hand, a progressive increase in contact stiffness affects greatly the wave propagation in the chain: the wave is more dispersed and more attenuated. This makes the impact outcome far away from the one given by the binary collision model. One can also see that for α = 6 and 10 (Figures 18.e and 18.f), secondary collisions occur at each contact, making the wave more scattered. With a very strong increase in contact stiffness (α = 10 4 ), the wave is strongly attenuated (Figure 18.h), leading to the impact outcome given by Moreau's law (Figure 16.f). The wave profiles shown in Figure 18 can explain the non-monotonic dependence of the dispersion measure C KE on α shown in Figure 19. The best dispersion effect (the minimum value of C KE ) is obtained for α = 20. For extreme values of α (a strong decrease or increase in contact stiffness), the dispersion effect is very small, and in these cases Moreau's law and the binary collision model can predict the impact outcome of the chain.

Varying the CoRs (dissipation)

We have studied so far the impact in purely elastic granular chains, i.e. there is no energy dissipation. We will show in the following how these systems behave when ample, for a collision between two beads constituted of chrome steel, which is a very elastic material, the energetic CoR e * defined in Section 4.3 is around 0.95 [START_REF] Nakagawa | Impulse dispersion in a tapered granular chain[END_REF]. Let us vary the energetic CoR e * from 1.0 (purely elastic case) to 0.0 (purely dissipative case) for the monodisperse chain considered in Section 5.1.1 with the elasticity coefficient η = 3/2. Figure 20 shows the impact outcome obtained with the LZB model compared to the one obtained with Moreau's law for different values of the energetic CoR e * . As mentioned in Section 3.2, when using the binary collision model for a dissipative monodisperse chain, there are more than one binary collision to be handled at one time. Two strategies has been proposed to handle these simultaneous collisions. However, the number of binary collisions can be infinite in many cases. Therefore, the binary collision model is not considered for the comparison in this section. We consider that two particles are stuck together if the absolute value of the relative velocity between them is smaller than 0.1% of the pre-impact velocity of the first ball. We define the value of e * under which particles are stuck together after impact. This value of e * is 0.86 for the considered chain of 100 balls, and it increases as the number of balls increases (Table 4). According to Moreau's law, the first ball bounces back and the other balls are stuck together after impact for any value of e * except for e * = 0 for which all the balls are stuck together. Therefore, the outcome given by Moreau's law is very different from the one given by the LZB model except for e * = 0 for which these two models give the same outcome. This result is confirmed in Figure 21 in which the gap measure Φ is plotted against the CoR e * . Figure 22 shows a comparison between the two impact models in term of the kinetic energy ratio KER defined as: KER = T + /T -with T + and T -being the kinetic energies before and after impact, respectively. When a granular chain with multiple contacts is subjected to an impact, the induced energy propagates and disperses in the system (Figure 23), which involves more contacts to participate in the impact process. If the system is dissipative, although each contact dissipates a small amount of energy, the whole system of multiple contacts dissipates a great amount of energy as shown in Figure 22. The wave is damped as it propagates through the system. For the considered chain with 100 balls, the energy induced by the shock is almost dissipated when e * < 0.9. With regard to Moreau's law, it underestimates the energy dissipation. Indeed, this impact model neglects completely the wave propagation in a system with multiple contacts so the impact is only localized at the first contact. It should be noted that Moreau's law describes the impact in the considered granular chain as a single impact between the first ball and another solid composed of the other balls. The only case where this impact law gives the same outcome as the one given by the LZB model is the purely dissipative case (e * = 0), for which the wave is strongly damped and the energy induced by the shock is almost localized at the first contact (Figure 23.f). Let us consider a decorated chain to investigate how the distribution of particle masses affects the impact outcome. The considered chain is composed of 101 balls whose masses are distributed as follows: the masses of balls with an odd number (1, 2, 3,...) are equal to m and the masses of balls with an even number (2, 4, 6,...) are equal to εm. The mass ratio ε is varied from 0.01 to 100. The Hertz's contact model (η = 3/2) and the energetic CoR e * = 1.0 are used for simulations performed with the LZB model. Figure 24 shows the impact outcome obtained with the LZB model for different values of the mass ratio ε. One sees that placing small balls between big balls makes the energy more distributed in the chain after impact (Figures 24.b and 24.d) except for very small or very big values of ε. When ε is very small, if we look only at the velocity of the big balls in Figure 24.a, the decorated chain behaves similarly to a monodisperse chain composed of the big balls (Figure 24.c). This means that separating big balls by very small balls does not significantly change the impact outcome of the big balls. In this case, curves representing the evolution of the potential energy at the two contacts on each small ball almost overlap and we find again the solitary wave which was observed for a monodisperse chain (Figure 14.d). When ε is very big (the first ball is very small compared to the second ball), the first ball bounces back with most of energy after impact. Concerning Moreau's law, its impact outcome for the considered chain is similar to the one of a single impact between the first ball and the remainder of the chain, independently of the mass distribution. Because the mass of the first ball is very small compared to the remainder of the chain, the first ball bounces back after impact with a velocity almost equal to its velocity before impact. This means that Moreau's law is not capable of predicting the effect of the mass distribution on the outcome of the decorated chain. The only case where this law gives the same impact outcome as the one given by the LZB model is for a very big value of ε (Figure 24.f). This is due to the fact that the big mass of the second ball compared to the first ball prevents the wave from propagating in the chain so the collision process is almost localized at the first contact as shown in Figure 25 Remark 9. The binary collision model is not used for this kind of chain because it leads to undefined impact outcomes after a huge number of binary collisions for several values of the mass ratio ε.

Figure 26 shows the dispersion measure C KE obtained with the LZB model and Moreau's law versus the mass ratio ε. It can be seen that the dispersion of postimpact kinetic energies of balls obtained with Moreau's law is weak and does not change significantly with the mass distribution. On the contrary, the LZB model predicts a strong effect of the mass distribution on the energy dispersion of the decorated chain. For this kind of chain, the energy induced by the shock is the best dispersed in the chain for ε = 0.64 and 1.55. As stated in [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF], the energy dispersion and the force transmission in a granular chain are related together. The first value (ε = 0.64) is quite close to the characteristic value ε = 0.59 shown in [START_REF] Jayaprakash | New family of solitary waves in granular dimer chains with no precompression[END_REF] for which the force transmission in a decorated chain is minimum.

Conclusions

For the tested systems of chains of aligned balls: Moreau's law has good predictive capabilities for small CoR (big dissipation), very small elasticity coefficient, big stiffness increase through the chain, or high mass ratio in decorated chains. In terms of waves, Moreau's law has good prediction capabilities when the wave is localized at the first contact. The binary collision model has good predictive capabilities for large elasticity coefficient, or large stiffness decrease through the chain. However it is very hard to draw conclusions with the binary collision law due to intrinsic issues like the impossibility to choose a unique order of collisions (different sequences usually yield different outcomes), and the lack of a criterion that guarantees its convergence (an infinite number of impacts is possible in some cases). For these reasons this approach should be disregarded most of the time.

Chains impacting a rigid wall

We have considered so far free granular chains where the first ball impacts the other stationary balls. In this section, we consider a monodisperse chain of 100 balls where all balls move with the same velocity and impact a rigid wall as illustrated in Figure 27. It is noteworthy that contrarily to the free chains, in this case the linear momentum of the 100 balls is not conserved. This kind of impact has been experimentally studied in [START_REF] Falcon | Collision of a 1-D column of beads with a wall[END_REF] and a good agreement between the numerical results obtained with the LZB model and the experimental results has been shown in [START_REF] Liu | Frictionless multiple impacts in multibody systems: Part II. Numerical algorithm and simulation results[END_REF]. It was observed that when the chain impacts the wall, the collision process starts at the bottom and then propagates to the top of the chain. The top ball leaves the chain first, is then followed by the next one and so on. The considered chain is composed of 100 elastic balls and the elasticity coefficient is varied. The balls are numbered from 1 at the top to 100 at the bottom. According to Moreau's law, all the balls are still stuck together and the chain moves upward after the impact with the same velocity as the one before impact. The same impact outcome is obtained with the binary collision model with a sequence of binary collisions from the bottom to the top of the chain to mimic the wave propagation. The post-impact velocities of balls obtained with the LZB model for different values of the elasticity coefficient are shown in Figure 28. It is shown that when Hertzian contact model is used (η = 3/2), balls are detached from each other after impact except few balls in the middle, and the top ball doubles almost its velocity. However, when the linear contact model is used (η = 1), about 70 balls in the middle are almost stuck together after impact, and this number of balls is about 80 for η = 0.1. It is expected that for a very small value of η, all the balls are stuck together after impact, which corresponds to the impact outcome given by Moreau's law. However, we were not able to simulate this impact problem for a very small value of η. In fact, when a chain of balls collides a wall, contacts undergo many repeated collisions as shown in Figure 29.b. As a consequence, integrating such an impact process is much more difficult than integrating the impact in a free monodisperse chain where each contact undergoes only one collision (Figure 29.a). One would expect that a value of η higher than 1.5 makes the top ball bounce back with a higher velocity. However, this is not the case for a monodisperse chain impacting a wall as shown in Figure 28.d where the post-impact velocity of the top ball for η = 2.0 is lower than for η = 1.5.

For a granular chain colliding a wall, a small dissipation at each contact can lead to a large damping effect as each contact undergoes many repeated collisions as mentioned above. As shown in Figure 30.a, a marked change in the impact outcome is observed for e * = 0.99, compared to the elastic case (Figure 28.c), and about 50% of energy is dissipated in this case. When e * = 0.96, the whole energy is dissipated and the chain is stuck to the wall after impact as shown in Figure 30.b. It is interesting to note that when the considered monodisperse chain impacts a wall, the whole energy is dissipated at higher value of e * than when it is free: as shown in Section 5.1.3, most of energy is dissipated for e * < 0.9 for the free chain.

Conclusions

In this chapter we have led comparisons between three classical multiple-impact laws: the binary collision model, Moreau's impact law and the LZB approach. The comparisons have been made on chains of aligned beads (free or impacting a wall) in terms of the post-impact velocities, kinetic energy dispersion, when the coefficients of restitution, the elasticity coefficients, or the contact stiffnesses are varied. The results given by the LZB model are considered as the reference. Waves propagation is known to be a crucial effect in such systems. We found that Moreau's law and the binary collision model, can predict with accuracy the impact outcome, only in few "extreme" cases (like very low or very high elasticity coefficient, mass ratio). Moreau's law applies well when the wave is localized at the first contact. Its advantage is that it is easy to implement, even in case of a great number of bodies and contacts. The binary collision law suffers from severe drawbacks like the possible infinity of impacts, different outcomes for different sequences of impacts, which make it very delicate to use for reliable computations in most cases. Future studies should focus on two-dimensional granular systems, with Coulomb's friction at contacts.

Fig. 3

 3 Fig. 3 Illustration of a decorated chain.

  -0.4444 0.2963 -0.3951 0.2634 -0.3512 0.7023

Fig. 4

 4 Fig. 4 Ball post-impact velocity versus ball number for a disordered chain obtained with the leftto-right and random sequences of binary collisions.

Fig. 5

 5 Fig. 5 Ball post-impact velocity versus ball number for a disordered chain obtained with four random sequences of binary collisions.

Fig. 8

 8 Fig. 8 Repeated impact: the contact point with two compression phases.

In sections 5 4 . 4

 544 .1.1, 5.1.2 and 5.1.3, chains with beads of equal masses are studied. Decorated chains are analysed in section 5.1.Algorithm Integration up to the end of each step k.

Fig. 9

 9 Fig. 9 Post-impact velocities obtained with LZB model versus ball number for different values of η.

Fig. 10 Fig. 11

 1011 Fig. 10 Post-impact velocities obtained with Moreau's model and binary collision model versus ball number.

Fig. 12

 12 Fig. 12 Gap measure Φ of Moreau's law and of the binary collision model versus η used in the LZB model.

Fig. 13

 13 Fig. 13 Number of balls for which the binary collision model and Moreau model give a good post-impact velocity.

Fig. 14

 14 Fig. 14 Potential energy E at the first 20 contacts versus time t for different values of η.

Fig. 15

 15 Fig. 15 Dispersion measure C KE versus elasticity coefficient η.

  α = 10 2 (f) α = 10 4

Fig. 16 Fig. 17

 1617 Fig. 16 Post-impact velocities obtained with LZB model versus ball number for different values of α.

α = 10 4 Fig. 18

 418 Fig. 17Gap measure Φ of Moreau's law and of the binary collision model versus α.

Fig. 19

 19 Fig. 19 Dispersion measure C KE versus coefficient α.

Fig. 20 Fig. 21

 2021 Fig. 20 Post-impact velocities obtained with LZB model versus ball number for different values of the energetic CoR e * .

Fig. 22 Fig. 23 Fig. 24

 222324 Fig. 22 Kinetic energy ratio KER versus coefficient of restitution e * .

Fig. 25

 25 Fig. 25 Potential energy E obtained with the LZB model at the first 20 contacts versus time t for (a) ε = 0.01 and (b) ε = 100.

Fig. 26

 26 Fig. 26 Dispersion measure C KE versus the mass ratio ε.

Fig. 27

 27 Fig. 27 Illustration of a granular chain impacting a wall.

Fig. 28 Fig. 29

 2829 Fig. 28 Post-impact velocities for an elastic monodisperse chain impacting a wall obtained with LZB model versus ball number for different values of η.

Fig. 30

 30 Fig. 30 Post-impact velocities for a monodisperse chain impacting a wall obtained with LZB model versus ball number for (a) e * = 0.99 and (b) e * = 0.96.

Table 1

 1 Bead velocity for a monodisperse chain after a sequence of binary collisions from the left to the right.

	1	2	3	4	5	6	7	8	9	10
	0.25 0.19 0.14 0.11 0.08 0.06 0.04 0.03 0.026 0.08

Table 3

 3 Number of binary collisions N c obtained with the left-to-right and random orders versus the number of balls n in a disordered chain with e n = 0.9.

	n	10 20 30 40	50	60	70
	N c -left-to-right order 53 189 991 18476 4731360 38936068 -
	N c -random order	51 153 397 1316	-	-	-

Table 4

 4 Value of e * under which particles are stuck together after impact for different values of the number n of balls.

	n 2 3 10 20 30 40 50 100 500
	e * 0 0.1 0.5 0.66 0.75 0.79 0.8 0.86 0.88

  .b.

		0.01							0.02
	E (J)	0.006 0.008						E (J)	0.01 0.015
		0.004
		0.002							0.005
		0	0	0.0001	0.0002	0.0003	0.0004	0.0005	0	0	0.0002	0.0004	0.0006	0.0008	0.001
						t (s)							t (s)

These developments make sense under some well-posedness conditions of the dynamics, which are assumed to hold here. In particular positions q(•) are absolutely continuous, velocities q(•) are right continuous of local bounded variations -hence possess right and left limits everywhere-, and accelerations are measures, as well as λ . See[START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF] Theorem 5.3] and[START_REF] Dzonou | A convergence result for a vibro impact problem with a general inertia operator[END_REF].

In Figures7 and 8the subscript i is not indicated. Thus R i is R, and so on.

Comparisons of multiple-impact laws for multibody systems