A time-stepping scheme for inelastic collisions Numerical handling of the nonoverlapping constraint

We propose here a numerical scheme to compute the motion of rigid bodies with a non-elastic impact law. The method is based on a global computation of the reaction forces between bodies. Those forces, whose direction is known since we neglect friction effects, are identiﬁed at the discrete level with a scalar which plays the role of a Kuhn-Tucker multiplier associated to a ﬁrst-order approximation of the non-overlapping constraint, expressed in terms of velocities. Since our original motivation is the handling of the non-overlapping constraint in ﬂuid-particle direct simulations, we paid a special attention to stability and robustness. The scheme is proved to be stable and robust. As regards its asymptotic behaviour, a convergence result is established in the case of a single contact. Some numerical tests are presented to illustrate the properties of the algorithm. Firstly, we investigate its asymptotic behaviour in a situation of non-uniqueness, for a single particle. The two other sets of results show the good behaviour of the scheme for large time steps.


Introduction
The present work was initially motivated by the necessity to handle particle collisions in the direct simulation of fluid-rigid particle mixtures.By direct simulation we mean that each particle is treated individually, and the Navier-Stokes (or Stokes) equations are solved in the moving domain occupied by the fluid.Many methods have been proposed recently to compute such flows.Some are based on an embedding of the solid phase in a global domain which is covered by a cartesian mesh (see Glowinski [9,Chapter 8] or Bertoluzza et al. [2]); the other class of methods consists in using a conforming mesh of the fluid domain (see Hu [12], Johnson et al. [13], Maury [16]).All these approaches face the problem of body overlapping, which makes it difficult to guarantee the robustness of the computation.Different strategies have been used.In Hu [12], the mesh is refined in the neighborhood of the interparticle gap, so that lubrication forces can be approximated with accuracy and consequently can be expected to prevent overlapping.Indeed, as the interparticle force between two smooth bodies separated by a viscous fluid behaves like −ε/ε (see Kim [14]), where ε is the distance, there cannot be any contact in finite time for reasonably regular external forces.Another commonly used strategy consists in adding short range repulsive forces between particles, which tend to prevent particles to overlap (see Glowinski et al. [10]).Those methods have proved to behave quite satisfactorily in many situations, but they necessitate a fine tuning of some numerical parameters, and the actual minimal distance between bodies cannot be controlled a priori.In [16], we introduced a heuristic method to control this minimal distance, by running at each time step a minimization procedure on a global functional of the particle positions, such that the minimal distance corresponding to the resulting configuration is greater than a prescribed "safety distance" ε > 0. This method, although robust in practice, even in the case of multiple contacts, is not consistent from the energy point of view, and its long-term effects on the computed flow are difficult to estimate.The approach we proposed in [15] is based on a first-order approximation of the lubrication forces exerted by the fluid in the interparticle gap.It is more respectful of the underlying physics, but the singular character of those forces (proportional to the reciprocal of the distance) introduces new constraints on the time step, and we must add that this approach, despite its very respect of local phenomena, does not seem to improve the overall modeling of fluid-particle mixtures significantly.It made us look for more basic tools, yet respectful of some basic physical properties like momentum balance and decreasing character of energy.
In the present work, we propose a new algorithm to handle the nonoverlapping constraint, which fits into the general class of Contact Dynamics methods, introduced by Moreau [19].The action of the surrounding fluid shall be simply described by a given external force field, so that from now on we consider "dry" systems of rigid bodies.Given the nature of the interaction between particles moving in a viscous fluid, which is dissipative, we chose to consider here purely inelastic collisions.The evolution problem on which our approach is based fits into the general framework of differential measure inclusions (see Schatzman [20], or Moreau [19]).From the theoretical point of view, the critical point is uniqueness, which necessitates strong regularity assumptions on the data: analyticity is required.Uniqueness is lost as soon as the external force is no longer analytic (see in Schatzman [20] or Ballard [1] counter-examples to uniqueness with a C ∞ force).This question of uniqueness in regard with the numerical algorithm will be addressed at the end of the paper.We refer to Stewart [23] for a general presentation of mathematical issues raised by those models, and a wide presentation of numerical methods.The most general result concerning convergence analysis of a numerical scheme for granular flows can be found in [22], where Coulomb friction is taken into account.
The situation we consider here is less general than the aforementioned one, but the scheme we propose presents the following features: 1. Collisions are not treated as events 1 .All contacts which are likely to occur during a time step are indeed handled globally, without any prior prediction of actual violations of the non-overlapping constraint, and reactions are computed as a field of Kuhn-Tucker multipliers associated to first-order approximations of those constraints, expressed in terms of velocities.2. As a consequence, this method can be proved to be unconditionally stable.For any choice of the time step, the kinetic energy of the approximate solution dissipates (when external forces are reduced to 0).Since there is no parameter to tune up (except for the stopping criterium associated to the Uzawa algorithm), this algorithm can be applied to very different situations with no modification.3. Provided the quadratic minimization step is solved exactly, the scheme produces feasible configurations only (in the case of spherical bodies).4. Non-smooth forcing terms can be taken into account, as soon as their mean values on subintervals can be computed.This makes it possible to handle multiple solutions numerically (see section 6.1).
The paper is organized as follows.In section 2 we present the rigid-body problem, and we introduce the natural functional spaces for positions, velocities and reaction forces.We then present the time-stepping scheme (section 3).In the next section we prove some general properties of the scheme, all of which are valid for any time step.Section 5 is devoted to the asymptotic behaviour of the scheme: as the time step h goes to 0, a subsequence of the computed solutions is shown to converge in some sense to a solution of the original problem, in the case of a single contact.Some numerical tests are finally presented, to illustrate the behaviour of the algorithm in different contexts: non-uniqueness, sticky particle model, and many-body motion (section 6).

Evolution model
We consider the mechanical system of N rigid spheres (or discs, or line segments) in R d (with d = 3, 2, or 1) of radii (r i ) 1≤i≤N and masses (m i ) 1≤i≤N .We introduce the configuration space, and the associated feasible set is the signed distance between spheres i and j (see Figure 2.1).Note that D ij is negative as soon as spheres i and j overlap.
As we do not consider rotations here, our configuration space is endowed with a natural Euclidean structure, so that the tangent space at any q ∈ Q can be identified with Q itself.We shall nevertheless denote by T Q = R dN this tangent space in order to distinguish velocities from positions.We define G ij = ∇D ij as the gradient of the distance between the two spheres i and j : and C q as the set of feasible directions at q ∈ Q 0 , and the outward normal cone to Q 0 at q as its polar cone: Let u = ( q1 , q2 , . . ., qN ) ∈ T Q denote the generalized velocity vector.We finally denote by f the force field acting on the spheres, and by M the mass matrix.The problem we are interested in may be phrased formally: where u − (resp.u + ) is the left (resp.right) limit of the velocity vector at time t, and P C q is the Euclidean projection onto the closed convex cone C q .This model will be given a proper mathematical framework in the next section.As regards its mechanical sense, equation (2), which can be read formally expresses the fact that overlapping is prevented by repulsive forces acting on each sphere along the normal vector at the contact point.When there is no contact, N q reduces down to {0}, so that (2) reads as the ordinary differential equation M q = f.Equation ( 3) provides the collision model.This is a particular case of a more general model (which includes possibly an elastic behaviour) where e is a restitution coefficient, set to 0 in the case of a non-elastic impact law.Indeed, as N q and C q are mutually polar, it holds I − P N q = P C q , where I is the identity.(see Moreau [18]).Note that, in the non-elastic model we consider, some energy is lost during each collision.This kinetic energy is actually transformed into heat within the bodies.The reader may refer to Frémond [8], where this increase in temperature is explicitly integrated into a general thermomechanical model for collisions between rigid bodies.
Remark 1 For 1D problems (like the one which is presented in Section 6.2), it is in fact more natural to define Q 0 as a connected component of the set of all configurations q with no overlapping.Indeed, the evolution problem does not allow q to pass from a connected component to another (the bodies cannot leap across each other).With this convention, the feasible set is closed and convex, and in this particular situation, N q identifies with the subdifferential of the indicatrix of Q 0 : and hence q −→ N q is a maximal monotone operator.In general Q 0 is not convex, but it is not far from being convex: it can be proved easily to be uniformly η-prox-regular in the sense defined in Colombo [7], as soon as 1/η < 4 min(r i ).

Functional framework
We consider the time interval I = (0, T ) and we introduce the following functional spaces: W 1,1 = set of all R dN -valued functions which are absolutely continuous over the interval I ; BV = set of all dN vector-valued functions having bounded variation over I : BV is the set of all functions t → u(t) ∈ T Q , such that each component u of u verifies sup where S = (t 0 , t 1 , . . ., t N S ) runs over the set of increasing subdivisions of the time interval I ; M 1 = set of N(N −1)/2 vector2 -valued bounded measures on I : it contains all μ = (μ ij ) 1≤i<j ≤N such that μ ij is a continuous linear functional over the set C 0 (I ) of continuous functions over I , vanishing at 0 and T .The set of component-wise positive measures will be denoted by As regards the forcing term, we shall consider functions of t and q, with Carathéodory type regularity conditions (see Coddington et al. [6]) : let the function f be defined in Q 0 × I , measurable in t for fixed q and continuous in q for fixed t, such that and f is uniformly Lipschitz with respect to the space variable in the feasible set: In order to avoid a special treatment of the initial time, we shall pick an initial position q 0 in the interior of Q 0 , so that N q = {0}.It follows that the velocity is continuous at 0 + , and u(0) is defined without ambiguity.
We may now state the problem: Given a time interval I = (0, T ), a force f (which meets regularity assumptions ( 6)-( 7)), and initial conditions (q 0 , u where P N q is the euclidean projection onto the closed convex cone N q .Equation (10) is to be understood in the sense of distributions (it identifies two 0-th order distributions).
Remark 2 Note that f(•, t) is required to be in W 1,∞ (Q 0 ) and not necessarily in W 1,∞ (Q).This makes it possible to apply this approach to model interaction forces like where V : (0, +∞) → R + is smooth, decreasing and convex, possibly singular at 0 : since the distance between any q ∈ Q 0 and 0 is kept apart from 0 (it is always greater than twice the smallest radius), such a function is indeed lipschitzian.For instance an electrostatic potential, V (d) = K/d, meets the requirements.
Remark 3 As q has the W 1,1 regularity in time, we can consider that q is continuous.Similarly, any u ∈ BV (or u h in the next section), shall be identified with the right-continuous element of the class.With this convention, all pointwise identities or inequalities will be meant in the classical sense (i.e.everywhere).
Remark 4 Obstacles like walls or fixed bodies may be taken into account.For the sake of simplicity, we shall not introduce here extra notations for that purpose, but the set of Kuhn-Tucker multipliers μ could be extended to include reaction forces against obstacles.

Approximation spaces
Let h = T /N h be the time step, and let P 0 h denote the set of all those functions which are piecewise constant over I according to the uniform subdivision (0, h, 2h, . .

. , (N h − 1)h, T ). Similarly, P 1
h denotes the set of continuous, piecewise affine functions according to the same subdivision.We introduce the following approximation spaces for trajectories, velocities, and reaction forces, respectively: Similarly, q n h denotes q h (nh), for any q h ∈ X h .As for external forces, we shall use the following mean value approximation f(q, s) ds.

Time-stepping scheme
The sequence of approximated fields (q h , u h , μ h ) ∈ X h ×V h ×R h is built according to the following scheme, to which we shall refer as scheme (S): 2. Compute u n+1 h as the solution to the constrained minimization problem min with The approximate reaction field is the dual component of a solution to the associated saddle-point problem 3. Update the positions Remark 5 The velocity field u n+1 h is uniquely defined (see (14): it minimizes a strictly convex fonctional over a closed convex set), whereas μ n+1 h is not.So in a strict sense scheme (S) is an algorithm as far as q and u are concerned only.
The following remark explains why the previous scheme can be expected to approximate the original problem.

Remark 6
Step 2, which is the essential (and time-consuming) part of the scheme, can be written where ∂I K h (q n h ) is the subdifferential of the indicatrix function of K h (q n h ).Now observe that, for any q ∈ Q 0 , N q is the sum of all half lines −R + G ij for indices verifying D ij (q) = 0, and ∂I K h (q) (u) can be written as the same sum for indices i and j such that As the latter quantity is a Taylor expansion of D ij (q + hu), the set ∂I K h (q n h ) (u n+1 h ) can be seen as a prediction of N q n+1 h .Therefore the scheme we propose can be considered, at least formally, as a semi-implicit time discretization of inclusion (2), Note that the collision law (3) does not appear explicitly in the scheme.It is actually implicitly contained in (20), which selects the velocity corresponding to the collision law we considered here.For this reason, this scheme is intrinsically dedicated to inelastic impact models, and we must say that there is no obvious way to adapt it to other impact laws.

Saddle-point problem
In the present approach we solve the saddle-point problem ( 16) by a Uzawa algorithm.Note that this choice is not mandatory, since the way the minimization procedure is performed is completely independant from the time-stepping scheme itself.Other algorithms could be implemented to perform this task, and the theoretical results presented in the next sections would remain unchanged.We introduce, for this section only, the vectors and matrices which are involved in the Uzawa procedure: where r = N(N − 1)/2 is the number of constraints.The problem can be put in the classical saddle-point form The Lagrangian of the problem is and any saddle-point of L, i.e. any couple (u, μ) verifying is such that u minimizes J over K.
The Uzawa algorithm (see Ciarlet [5]) consists in approximating a saddle-point (u, μ) of L by sequences (u k ) and (μ k ).A step decomposes into two substeps: 1. Solve the primal problem 2. Update the Kuhn-Tucker multipliers where + is the orthogonal projection onto R r + : Proposition 1 We suppose that where α is the smallest eigenvalue of the matrix M. Then the sequence (u k ) converges to the solution u to the constrained minimization problem.The sequence Proof The convergence of the primal sequence (u k ) is a classical result (see Ciarlet [5]).As for the Kuhn-Tucker multipliers, the sequence (μ k ) can be shown to meet the assumptions of Opial's Lemma (see Haraux [11]), and thus it converges (weakly) to some μ in the dual solution set.As the problem is finite-dimensional, the convergence is strong.
We conclude this section by two remarks on the actual implementation of this algorithm.
Remark 7 It is not necessary to assemble the whole matrix C. At a given time step, positions and velocities are known, so that most constraints (which correspond to particles far away from each other) can be eliminated a priori, in the spirit of what is presented in Sigurgeirsson [21].For monodisperse situations (all radii are identical), the total number of potentially active constraints (i.e.number of rows of C) is κN/2, with κ = 6 in two dimensions (each disc is close to 6 other discs, at most), and κ = 12 in three dimensions.
Remark 8 A sufficient condition on ρ for the algorithm to converge can be obtained explicitly in standard situations.We give here a rough lower bound for ρ max in the monodisperse case.In this situation, the maximal number of likely contacts per body is smaller than the aforementioned parameter κ.The ellipticity constant α is simply min(m i ) 1≤i≤N .Using expression (1) of G ij , we can write C T C ∈ M dN (R) by blocks: where C ij is a d × d matrix defined by with J i denoting the set of indices of the spheres which may get into contact with sphere i.By Gerschgorin circle theorem, the spectral radius of C T C (i.e.||C T C|| 2 ) is smaller than the maximum among the 1-norms of its rows.Since e ij = 1, it follows that ||C T C|| 2 ≤ 2κh 2 √ d, where κ is the upper bound of (J i ) (cardinal number of J i ).Finally, where κ is 12 for 3D problems, as mentioned in the previous remark.It is to be noted that this lower bound is independent of N , the number of bodies.

Properties of the scheme
This section is devoted to some general properties of the algorithm in the general case (multiple contacts are allowed).Some of those properties will be used in the next section to prove convergence in the case of a single contact.
Proposition 2 (Feasibility) If the minimization step is solved exactly, then the scheme (S) produces feasible configurations only: Proof As Step 2 of Scheme (S) is supposed to be solved exactly, then for any h, n, i = j , We established that q n h ∈ Q 0 for any h, n.Now for any t = nh+θh, with θ ∈ [0, 1], one has Note that this property is no longer true if the bodies are not spherical, because the functions q −→ D ij (q) are not convex in general.
Proposition 3 (Stability) Let (q h , u h ) h>0 be a sequence built according to Scheme (S), for a given f.Assuming again that the minimization procedure (Step 2) is performed exactly, then u h verifies the following estimate: where |•| M is the euclidean norm associated to the symmetric positive definite matrix M: F is the function which dominates f (see condition ( 6)), and α is the smallest eigenvalue of M (α = min(m i )).
Proof Let us first establish that, for every n, i = j , It is a direct consequence of the fact that a Kuhn-Tucker multiplier may be different from zero only if the corresponding constraint is active.Indeed, either μ n+1 ij = 0 (the constraint is non-active), or, if μ n+1 ij > 0, then necessarily (the constraint is active) We perform the scalar product of the relation (18) with the velocity field u n+1 h : We finally get, by summing up over steps 0, 1, . . ., n, which gives inequality (24).
Remark 9 In the case f ≡ 0, the previous proposition expresses the fact that kinetic energy decreases at the discrete level, for any time step.
Proposition 4 (Momentum balance) If the forcing term does not depend on q, then the total momentum balance is exactly verified at the right end of each subinterval by the approximated solution.
Proof Let us denote by P(u) the total momentum associated to the velocity field u = (u 1 , . . ., u N ): and by f i the force acting on sphere i.As the approximated reaction forces verify the law of action and reaction, it follows from (18) that where f = f i is the resultant force, so that which is exactly the momentum balance.
Remark 10 Note that the previous property remains valid even if the saddle-point problem is not solved exactly.
Remark 11 As a direct consequence of the previous property, the computed trajectory q h of the center of mass of the system turns out to be an interpolation of its exact trajectory (which is uniquely defined, even if particle trajectories are not), which can be expressed with obvious notations We shall complete this section by a proposition concerning the Kuhn-Tucker multipliers.As u n+1 h is uniquely defined, so is the global reaction term As for the vector μ n+1 h itself, uniqueness does not hold in general, as shown by the following example: considering a two-dimensional collection of N identical discs in a cristal-like configuration such that most discs are in contact with 6 others, the number of active constraints is asymptotically 3N, whereas the number of degrees of freedom is 2N .This non-uniqueness of the vector of Kuhn-Tucker multipliers is known to lead to numerical instabilities in some situations.We show next that the risk of instabilities is actually limited, because the solution set is bounded, as will be deduced from the following lemma.

Lemma 1 We consider
+ , and we define The set Proof Let us first establish the uniqueness for the homogeneous problem.For Let i 0 denote the index of an extremal vertex of the convex hull conv{q i , 1 ≤ i ≤ N}.By Hahn-Banach's theorem, the compact {q i 0 } and the closed convex set conv{q i , 1 ≤ i ≤ N, i = i 0 } can be separated in a strict sense by a plane in R d .We denote by x an element of this plane, and by v a normal vector to it.One has so that (q i 0 − q j ) • v > 0 for j = i 0 .Now the balance of contact forces exerted upon sphere i 0 in the direction v reads which is positive unless λ ji 0 = 0 for all j = i 0 .Therefore all multipliers associated to a contact with sphere i 0 are equal to 0, and this approach can be iterated for the reduced family (q j , j = i 0 ).By downward induction on the number of active spheres, we establish that 0 is actually reduced to {0}.
As for the non-homogeneous problem, F is obviously convex.The asymptotic cone of F is defined as (see e.g.Bourbaki [3]) where λ is any element of F (for example μ).It can be checked that (the proof, which is elementary, can be found in Maury [17]), in a finite-dimensional space, the asymptotic cone of a convex set is reduced to {0} if and only if the set is bounded.If F is not bounded, then its asymptotic cone C F contains a half line R + ξ , with ξ ∈ R r + , which yields μ + R + ξ ⊂ F , and consequently ξ is a solution to the homogeneous problem, so that ξ = 0. Therefore, F is bounded.

Proposition 5 At each time step of Scheme (S), the solution set for
Proof This is a direct consequence of lemma 1, with q = q n h and

Convergence theorem (case of a single contact)
We establish in this section a convergence result for the proposed time-stepping scheme, in the case of a single constraint (two spheres) in R 3 .As uniqueness does not hold, neither does continuity with respect to initial conditions.Thus we cannot expect any error estimate for the scheme (S).Moreover, the necessity to extract a subsequence to establish convergence is not purely technical: different subsequences may actually converge to distinct solutions to the same problem (see Section 6.1).
In order to alleviate notations, we consider the system of two spheres with unit mass and same radius r, but this assumption is not needed in the proof.The result remains valid for the case of two different spheres, or a sphere and an obstacle.The constrained minimization problem is again supposed to be solved exactly at each time step.

Two-body system
We consider here the mechanical sytem of two identical, hard spheres, described by and the associated feasible set where D(q) = |q 2 −q 1 |−2r is the distance between the two bodies.As previously, the tangent space at any q ∈ Q is denoted by T Q = R 6 .The convex cone N q is now where G = ∇D is the gradient of D. In order to avoid special care of the initial instant, we shall pick an initial configuration q 0 such that D(q 0 ) > 0.
We may now state the single-contact problem: given a time interval I = (0, T ), a force field (q, t) → f(q, t) ∈ R 6 verifying the regularity conditions ( 6)-( 7) and the initial conditions where ( 28) is understood in the sense of distributions, and P N q is the Euclidean projection onto the closed convex cone N q .
Let us now rewrite the overall time stepping-scheme in the present case of a single contact.The discretization spaces for pathlines (X h ), velocities (V h ), and reactions (R h ), are For a given h > 0, we shall consider built according to Scheme (S): 1. Initialization 2. Time step where u n+1 h is the solution to the constrained minimization problem min and μ n+1 h ∈ R + is the associated Kuhn-Tucker multiplier.

Convergence
Theorem 1 Let (q h , u h , μ h ) h be a sequence of solutions obtained by Scheme (S) (equations ( 31) to ( 35)), with h → 0. Then there exists a subsequence of time steps (still denoted by h), and and (q, u, μ) is a solution to problem ( 26)-( 30).
Proof The proof will be decomposed into 9 steps, which we shall briefly describe below.The critical ones are steps 3 and 9.

The family
3. The fields u h have uniform bounded variation, i.e.,

The familly (u h
) is relatively compact in L 1 (I ).One can extract a subsequence (still denoted q h ) such that q h −→ q in W 1,1 .The limit velocity u = q = lim qh is in BV, and the limit motion q is feasible.5.The sequence (μ h ) h is bounded in L 1 , so that one can extract a subsequence which converges weak-to a vector-valued bounded measure μ ∈ M 1 .6.The pair (u, μ) verifies (28).7. The complementarity slackness condition (29) holds true. 8.The initial condition (26) is verified.9.The jump equation (30) is verified.
Step 1 (Feasibility) This is exactly Proposition 2, which was established for any number of contacts: Step 2 (u h is uniformly bounded) By Proposition 3, Step 3 (u h has uniformly bounded variation) We shall first establish a lemma, which can be seen as a mono-dimensional version of the property.

Lemma 2 We consider sequences (u h ) h , (ρ h ) h , and
We suppose furthermore that there exists a constant C ≥ 0 such that Then the familly (u h ) has uniform bounded variation: Proof A straightforward calculation leads to with λ k h ∈ [0, 1] for any h, k.Therefore function u h can be written as a sum of the homogeneous solution v h defined by

and a function w h determined by
The first one v h is monotonous and keeps a constant sign, so that var(v h ) ≤ |u 0 |.

The second one verifies var(w
which ends the proof of Lemma 2.
Remark now that, as u n+1 h is defined as the Euclidean projection of where α ≥ 0 is some real parameter.If the constraint is non-active (i.e.u n h + hf n+1 h (q n h ) ∈ K h (q n h )), then α = 0. Otherwise, as D(q n h ) ≥ 0, where λ n h is a non-negative real number.Let us show that λ n h ≤ 1.Again, as D(q n h ) is non-negative, 0 lies in K h (q n h ), therefore and therefore λ n h ≤ 1.We established that

Let us now introduce the Euclidean decomposition
where v n h is the projection of u n h onto G⊥ .We verify now that U n h meets the assumptions of Lemma 2 : • G(q n+1 h ) − G(q n h ) can be written as hg n h , and the corresponding family of piecewise constant functions (g h ) verifies (36).As ρ n h ∈ [0, 1], all assumptions of Lemma 2 are met, which implies that var(U h ) is bounded uniformly in h.
We denote by P n h the projection onto G(q n h ) ⊥ , so that v n h = P n h u n h .We conclude this step by writing Step 4 (q h converges to a feasible trajectory) As (u h ), which is bounded in L ∞ , is now proved to have uniform bounded variations over I , it is relatively compact in L 1 , so that one can extract a subsequence, still denoted by (u h ), which converges in L 1 to u.As u h is the Sobolev derivative of q h (in W 1,1 ), then q h converges to a q in W 1,1 , and q = u.The limit velocity u is in BV.Indeed, for any ϕ ∈ (C 1 c (I )) 6 , Finally, as the convergence q h −→ q is uniform (W 1,1 → L ∞ ), and D(q h (t)) ≥ 0 for every h, t, with D continuous with respect to q, then D(q(t)) ≥ 0 for every t ∈ I .

Step 5 (μ h converges weak-to μ)
For any h, so that (μ h ) is bounded in L 1 (|G| is the modulus of G(q h ), which is a positive constant).Therefore, by the De La Vallée Poussin compactness criterion, one can extract a subsequence which converges weakly-to μ in M 1 .
Step 6 (The balance of momentum is verified by the limits u and μ) Let us first remark that the derivative of u h can be defined3 in the sense of distributions as where δ n is the Dirac mass at t n = nh.For any ϕ ∈ (D(I )) 6 , it comes when h goes to 0. The integral I uh • ϕ may also be expressed where ϕ h ∈ V h is defined by taking the left value of ϕ on each subinterval [t n , t n+1 ), and similarly q h stands for the function of V h which is equal to q n h on [t n , t n+1 ).As f is uniformly Lipschitz with respect to q, and as q h converges to q in L ∞ , f h (q) converges to t → f(t, q(t)) in L 1 .Similarly ϕ h converges to ϕ in L ∞ , so that For the second term, we remark that (μ h ) can be seen as a bounded sequence of linear functionals on L ∞ .According to this remark, one may write (G • q designs t −→ G(q(t)), which is the uniform limit of G

Finally we have
From ( 37) and (38), it follows for all ϕ ∈ D(I ): this is equation (28) in the sense of distributions.
Step 7 (The limit reaction field μ is non-active when there is no contact) For any t 0 ∈ I such that D(q(t 0 )) > 0, there exists η such that so that, by uniform convergence of q h to q, there exists h 1 such that Since the velocities u h are uniformly bounded in L ∞ , there exists h 2 ≤ h 1 such that, for any h ≤ h 2 , D(q h ) + hG(q h ) • u h ≥ a/4 > 0, so that the numerical constraint D + hG • u ≥ 0 is not activated in (t 0 − η, t 0 + η).
Step 8 (The initial condition is verified) As q 0 is chosen apart from the boundary of Q 0 and thanks again to the uniform boundedness of u h in L ∞ , there exists η such that μ h vanishes in (0, η) for all values of h.It holds then Step 9 (The collision law is verified) For t 0 such that D(q(t 0 )) > 0, we already established (see Step 7) that μ ≡ 0 in a neighbourhood of t 0 , so that u is absolutely continuous in this neighbourhood (as it solves (28)), hence continuous at t 0 : The interesting case is of course D(q(t 0 )) = 0. We first establish from (28) that the discontinuity occurs along direction G : Lemma 3 At any contact time t 0 there exists λ ∈ R such that where G is the normalized gradient G/ |G| introduced in step 3.
Proof For η > 0, and H ∈ T Q such that H • G = 0, we consider the test function ψ η H where ψ η is continuous, piecewise affine, zero outside [t 0 − η, t 0 + η], taking values 0, 1, 0 at t 0 − η, t 0 , t 0 + η, respectively 4 .Now writing (28) against ψ η H yields The left-hand side can be expressed The force integral I ψ η H • f goes to zero with η.The last term can be bounded where the sup goes to zero because G • q(t) tends to G(q(t 0 )) as η tends to 0, with H • G(q(t 0 )) = 0.So finally the jump (u + − u − ) • H is zero for any H orthogonal to G, which ends the proof of Lemma 3.
Unless explicitly mentioned, fields are now taken at t 0 , where t 0 is a time at which contact occurs, and G designs G(q(t 0 )).We first check that u − • G ≤ 0. If otherwise, then D(q) = 0 implies that is negative for some ε > 0, meaning that q(t 0 − ε) / ∈ Q 0 , which is impossible.So which implies similarly that q enters the forbidden domain Q C 0 right after t 0 .Therefore Let us now notice that the condition (30) characterizes the velocity after the collision as the solution to a constrained minimization problem : it minimizes u − u − over C q , the polar cone of N q , which is the closed half space As u − • G ≤ 0, the solution to that problem is clearly which we have to identify with u + to prove that the collision law (30) holds at time t 0 .As we know that u + is u − + λ G with condition (39) on λ, we just have to check that We extend notations of Step 3 by introducing U ∈ L 1 (which can be seen as the velocity in the moving G direction) as the limit of u h • G(q h ).Now we recall the property we want to establish : at any t 0 such that D(q(t 0 )) = 0, equation (40) is verified.As t → G(q(t)) is continuous, it can be formulated in terms of U : at any contact time, U + ≤ 0. We propose actually to establish the property (which is stronger as U − ≤ 0 at any contact time) : Proof of (41) : We first recall an important property of U h = u h • G(q h ) which was established in Step 3 : there exists sequences g n h and ρ n h such that with ρ n h ∈ [0, 1] and |g h | ≤ g for some g ∈ L 1 .By summing up over time steps between t and t + τ , it comes and ρ(h, t, τ ) is in [0, 1].Let t 0 be such that U − (t 0 ) ≤ 0.Then, for any ε > 0, there exists η 1 such that Now, as g ∈ L 1 (I ), from (43) there exists η 2 ≤ η 1 such that, for any t ∈ (t 0 −η 2 , t 0 ), for some measurable set A h ⊂ (t 0 −η 2 , t 0 ) with |A h | → 0. Using (42) to "translate" approximately inequality (44), it comes The last inequality and L 1 convergence of U h to U (now considered in the translated interval (t 0 , t 0 + η 2 )) imply so that U + (t 0 ) ≤ 3ε, for any ε > 0. Therefore we have This completes the proof.

Numerical non-uniqueness
This first set of results is somewhat distant from the original purpose of this work, but it illustrates an interesting property of Scheme (S), which is the capability to compute multiple solutions.We consider the situation of a single material point moving on the real line, subject to the constraint q(t) ≥ 0, with an inelastic impact law at 0 + .The system reads : where N q is {0} whenever q > 0, and R − as soon as q = 0. We shall consider a force field f which does not depend explicitly on q, and which is defined in the time interval I = (0, 4) by (see Figure 2) We chose α = 1.8.When the time step h is a negative power of 2, the right-hand side involved in the time-stepping process can be computed exactly.We present two solutions which have been computed with h = 2 −8 and h = 2 −9 , respectively.Those approximations correspond to distinct piecewise polynomial functions, which can be shown to be both solutions of the problem.One can observe that for time steps of the type 2 −2n , we have convergence to a solution (Figure 3), and for h = 2 −2(n+1) , we have convergence to the other solution (which corresponds to Figure 4).Other sequences of time steps were not investigated (f -integral can no longer be computed exactly).

Sticky particles
This second set of computations illustrates the good behaviour of the scheme when the time step is large.The bodies we consider are punctual masses moving on a line, with no external force, so that when two particles meet, they collide and behave like a single particle.We consider the following situation: the initial condition is q 0 = (0, a, 2a, . . ., 1 − a, 1) with a = 1 N − 1 , the initial velocity vector u 0 is chosen arbitrarily in [−1, 1] N , and particles move according to system (8)- (12).
Remark 12 This one-dimensional model is now commonly referred to as sticky particle model (although particles do not stick in a strict sense, as they could be pulled apart with infinitesimal forces).It is used in Brenier [4] to establish existence of solutions to the pressureless gas model same time, so that the individual behaviour of particles is poorly described.Nevertheless, all these computations give exactly the same final state (up to errors due to the constrained minimization procedure).This behaviour is actually a direct consequence of Proposition 4 and Remark 11.Indeed, in this one-dimensional model, the final state (when all the masses are stuck together) is completely determined by the position of the aggregate and its velocity.As the center of mass moves with a constant velocity, the scheme is exact for its trajectory so that, no matter how large the time step is (as soon as all collisions have been taken into account), the approximated final state is exact.

Many-body computation
Numerical tests have been performed satisfactorily for large numbers of spheres (up to 50 000).We present here its application to a polydisperse collection of 2000 spheres.For readability reasons, we chose a situation such that the spheres remain Fig. 8 Step 0 Fig. 9 Step 3 Fig. 10 Step 6 Fig. 11 Step 9 in a fixed plane of R 3 , so that they behave like discs in R 2 , but the algorithm can be applied indifferently to genuine 3D problems.Radii are distributed in the interval [0.8r 0 , 1.2r 0 ], with r 0 = 1.5 × 10 −2 , and the mass m i of sphere i is proportional to its volume.
As initial condition, we consider a cluster of spheres with uniform downward velocity U 0 = −e z .We simulate the impact of this cluster with a rigid obstacle defined as the half space R × R × R − .To illustrate the stability of the scheme, we run a case where the time step h = 0.1 is such that the displacement of a sphere during one time step is typically 3 times its own size.Figures 8, 9, 10, and 11 represent the 2000 spheres at steps 0, 3, 6, and 9.

Conclusion
We presented a scheme to compute the motion of rigid bodies with nonelastic impact law.Because of its stability properties, its robustness (it produces feasible configurations only, even for large time steps), this scheme is particularly suitable to control the minimal distance between rigid particles in the context of fluid-particle simulations, with a controlled influence of the perturbation on the energy.
Beside this particular purpose, the fact that it can be used to handle situations where uniqueness does not hold, and its special behaviour in the case of simultaneous or quasi-simultaneous collisions even for large time steps, should make it an efficient tool to model more general granular flows.

Figures 5 ,Fig. 5 Fig. 6 Fig. 7
Figures 5,6 and 7 represent the computed pathlines of the particles in the time interval [0, 3], for a "small" time step h = 0.01, and two larger time steps h = 0.5 and h = 1.5.For h = 0.01, all collisions are distinctly captured by the scheme.For larger time steps, several collisions are taken into account by the scheme at the