A Guiot

S Pattofatto

C Tournier

L Mathieu

Modeling of a polishing tool to simulate material removal

In plastic injection mould and prosthesis industries, "mirror-effect" polished surfaces are required for obtaining transparent parts or surfaces without scratches. Traditionally done manually, we have proposed to automate polishing on 5-axis machining centre using a passive elastomeric carrier. One of the main advantages of automatic polishing is the repeatability of the machine movements in order to achieve restricted form deviations. However, the material removal rate (MRR) during polishing depends on parameters such as contact pressure, relative velocity and tool wear. We have thus developed a model dedicated to our process to compute the effective MRR along the polishing tool path regarding the contact area and the contact pressure between the tool and the part.

INTRODUCTION

The development of High Speed Machining (HSM) has dramatically modified the organization of plastic injection molds and tooling manufacturers. HSM in particular reduces mold manufacturing cycle times by replacing spark machining but does not enable to remove the polishing operations from the process. As the usual polishing process is expensive in terms of price and downtime of the mold, we have developed an automatic polishing process on 5-axis machine tool [START_REF] Pessoles | Automatic polishing process of plastic injection molds on a 5-axis milling center[END_REF]. Our approach consists in using the same tools as those used in manual polishing for pre-polishing and finishing polishing (Figure 1). Pre-polishing is performed with abrasive discs mounted on a carrier. This carrier is a deformable part made in an elastomer material fixed on a steel shaft that allows mounting in the spindle. Hence, we do not have a force feedback control but we manage the polishing force with the length of the tool introduced in the numerical controller. The polishing strategies are for the most part issued from previous experiences as for trochoidal tool paths (figure 1) computed on fractal trajectories coming from robotized polishing or cycloidal weaving paths representative of manual polishing [START_REF] Mizugaki | Fractal path generation for a metal-mold polishing robot system and its evaluation by the operability[END_REF], [START_REF] Tam | Robotic polishing of free-form surfaces using scanning paths[END_REF], [START_REF] Tsai | Efficient automatic polishing process with a new compliant abrasive tool[END_REF]. The main advantage of this process is the repeatability of the machine movements to respect the specified form deviation. Indeed, in manual polishing the material removal rate (MRR) depends on the polishing tool velocity and on the polishing pressure applied by the operator. An adequate polishing force facilitates the removal of cusps and stripes left on the part during milling or previous polishing operations. Nevertheless, the contact pressure has to be as constant as possible to avoid over-polishing and respect form deviation tolerances. Furthermore, even if the effective MRR is quit constant, the polishing time spent in every area of the part and consequently the polishing velocity along the tool path has to be constant and the distribution of the tool path has to be consistent. The aim of the paper is thus to propose a tool to simulate the material removal during polishing. The final objective is to provide a tool to optimize polishing parameters along the tool path (spindle speed, feedrate, imposed dis-placement) and the tool path itself in order to manage geometrical deviations.

In the next section several models of the material removal rate from the literature are presented. Then in section 3 we will present our simulation tool and the tool wear evaluation. Section 4 is dedicated to the modeling of the behavior of the disc carrier by using Finite Element Analysis. Simulations are analyzed in section 5 through experimental investigations before concluding remarks in section 6.

RELATED WORKS

The measurement or quantification of abrasion is most frequently given by the material removal rate (MRR), that is to say, the thickness of material removed per time unit. We can find two different approaches to model the MRR, the analytical models and the experimental models [START_REF] Nasri | A process model for robotic disc grinding[END_REF].

The analytical models are based on modeling the interaction between the tool and the workpiece at the level of the abrasive particle [START_REF] Bai | A mathematical model for material removal and chemical-mechanical synergy in chemical-mechanical polishing at molecular scale[END_REF] [7] [START_REF] Oh | An integrated material removal model for silicon dioxide layers in chemical mechanical polishing processes[END_REF]. We find these approaches in the electronics industry in particular for wafer polishing. The proposed models give an estimate of the material removal rate from all material characteristics (Young's modulus, hardness, Poisson's ratio), shape of the surface (type, size, distribution of asperities and particles) and the operating conditions (pressure, relative velocity, the effect of a chemical or fluid).

There is no comparable model for the realization of mechanical parts. The models used for polishing of these materials are rather experimental models. The experimental models are derived from the analysis of many polishing trials. The model of Preston [START_REF] Preston | The theory and design of plate glass polishing machines[END_REF], developed within the context of glass polishing, is probably the one that remains today the most used (eq.1). It states that the material removal rate is proportional to the average pressure of contact, P, and to the tool/part relative velocity V:

MRR = dz dt = k ⋅ P ⋅V (1)
The coefficient of proportionality also called "Preston coefficient" K is determined experimentally and depends on process parameters (part material, abrasive, lubrification, etc.).

In their work, Luo et al. [START_REF] Luo | Modification of the Preston equation for the chemical-mechanical polishing of copper[END_REF] studied the influence of the fluid for copper polishing. In addition to mechanical mechanisms, they highlight the fact that there are chemical mechanisms involved in abrasion. The proposed model is a slight evolution of the Preston model:

MRR = (K 1 . P + B) V + R c (2)
Where K 1 , B and R c are three constants to be determined experimentally. R c is the rate of material removal due to the chemical mechanism. Wang et al. [START_REF] Wang | A material removal model for polishing glass-ceramic and aluminium magnesium storage disk[END_REF] also proposed an extension of the Preston model for the manufacture of storage disk made of aluminum and magnesium alloy or glass-ceramic material:

MRR = K 2 . P n/2 . V 1-n/2 (3)
K 2 and n are two constants to be determined experimentally. This model seems to give good results for different materials.

Finally, analytical models have the advantage of predicting the removal material rate from all process data, which can also be a drawback because the knowing of all the physicochemical properties of the polishing components is mandatory. The experimental models are used for diverse materials. Preston's equation was originally used for polishing glass, but is now used for metal or semiconductor polishing. Whatever the model, the MRR always depends on pressure and relative velocity between the tool and the part.

However, literature models do not take into account the wear of the abrasive tool. Furthermore, these models have been developed for planar polishing on automatic machines. Thus, the contact area between the abrasive disk and the part is always a plane so that the contact pressure is uniform over the part. The relative velocity between the part and the tool is a combination of the tool and part rotation speed. In automatic polishing, the cutting speed Vc of the tool is one hundred larger than the tool feedrate Vf, and Vf is not constant when polishing on a 5-axis machine tool [START_REF] Lavernhe | Kinematical performance prediction in multi-axis machining for process planning optimization[END_REF]. The contact area and the pressure distribution could be determined by finite element computation, taking into account the disk carrier material behavior, inclination and position of the tool. This will be developed in section 4. In a first approach, we consider that the pressure distribution in the contact area is uniform and the form and size of the contact area are supposed constant. This allows the definition of our simulation tool presented hereafter.

PROPOSED MODEL FOR MATERIAL REMOVAL RATE SIMULATION Our simulation tool is based on the Preston model to compute the cutting depth per time unit. It takes into account effective cutting speed at each point of the contact area as well as the effective feedrate.

A global model is first proposed and in a second step, the tool wear, the geometry of the contact area and the pressure distribution are investigated and included in the model.

Global approach

As mentioned in the previous section, there are many abrasion models, but in mechanical and particularly for the polishing of molds and dies, the easiest model to use is based on the Preston equation (eq.1). Using this equation for a little element of surface that can be assumed to be a plane of normal z, it is equivalent to write during a little time Δt:

Δz = K . P . Vc . Δt (3)
The cutting depth Δz is proportional to the polishing time spent on an area of the part, which is in contact with the disk. If the pressure and the velocity are kept constant during polishing, one can think that the MRR will be constant. Actually, the disk efficiency is different over the time because of its wear. Hence, the Preston coefficient K depends on time and it will be decomposed by the product of the initial Preston coefficient K 0 =K (t=0) and a decreasing time function representative of the efficiency of the abrasive disk (wear law) w(t).

The Preston coefficient K(t) can be expressed as follow:

K(t) = K 0 . w(t) (5)
Hence the Preston equation becomes:

Δz = K 0 . w(t). P . Vc . Δt (6)
If we consider a constant cutting speed, we can define a nominal material removal rate MRR 0 :

MRR 0 = K 0 . P. Vc = cst (7) Δz = MRR 0 . w(t) . Δt (8)
In order to go further in the definition of the model, we need to define the contact area. In a first approach, the contact area is modeled by an ellipse with a minor axis 2a and a major axis 2b (figure 2) and the pressure distribution is supposed to be constant.

The first step of the method consists in meshing the surface to polish and sampling the tool path according to a maximum length between points l_max, initially fixed by the user.

The result is a list of n points M(i) such as:

|M(i+1)-M(i)| < l_max (9)
Figure 2: Tool/surface contact.

Then for each point M(i) of the tool path, the contact area is positioned such as the point O (tool location) corresponds to the point M(i) and oriented along the inclination direction θ (figure 3). When the contact area is in place, any node j among the m nodes of the mesh is marked in a (n,m) dimension matrix Ins(i,j). Any point included in the contact area (red points) takes the logical value 1 and all others take the logical value 0. Along all segment i included between points M i and M i +1 , the feedrate Vf and wear law w(t) are supposed to be constant. Thus the segment travel time can be approximated by:

Δt i = M i +1 -M i Vf i (10)
Taking l_max small enough allows making the assumption that the points in contact for the position M i still in contact during the entire segment i. Now it is simple to obtain the abrasion time T(j) for each node of the mesh:

T( j ) = Δt i i =1 n ∑ ⋅ Ins( i , j ) (11)
We define the local abrasion time for each node T w ( j ) weighted by the tool wear:

T w ( j ) = w( t i )⋅ Δt i i =1 n ∑ ⋅ Ins( i , j ) (12)
Finally, the cutting depth D(j) at the node j is obtained by the relation:

D( j ) = MRR 0 ⋅ w( t i ) i =1 n ∑ ⋅ Δt i ⋅ Ins( i , j ) (13)
We thus have exposed a model to predict the cutting depth based on a theoretical contact area between the tool and the part modeled has an ellipse with constant pressure distribution. The next section focuses on the definition of the tool wear.

Tool wear identification

The tool wear has been identified experimentally by polishing plane stripes and measuring the polishing depth along the stripes. Since the feedrate has been kept constant during the test, the correlation between efficiency and time can be established. The tool wear evolution plotted in Figure 4 The next section focuses on the definition of the effective contact area by using finite element analysis.

ELASTOMERIC CARRIER MODELING

The purpose of this section is to model and simulate the mechanical behavior of the disc carrier when pressed on the workpiece. The geometry of the contact area and pressure distribution at the interface with the workpiece are computed regarding the imposed displacement and the angle between the tool axis and workpiece. The steel tool shaft is modeled as a virtual part of infinite rigidity in the simulation. In this section, the abrasive disk is not taken into account. Simulation results are compared to experimental investigations.

Numerical simulation set up

The disc carrier is made of elastomeric material the composition of which is unknown exactly.

Traction tests allowed us to determine the value of Young's modulus E (9 MPa). The Poisson coefficient has been set to 0.49 as much of the elastomeric materials. We consider a purely elastic behavior of the disc carrier. The connection between the disc carrier and the shaft is modeled as a rigid joint and the part is considered as infinitely rigid. The desired displacement is imposed along the shaft axis. The disc carrier has a diameter of 18mm and is meshed with tetrahedral elements of maximum size equal to 0.5 mm (Figure 5).

Experimental validation

To validate the numerical simulation, the tool has been photographed during its deformation on a plate of polycarbonate. The following pictures show the comparison for imposed displacements equal to 0.6, 1.4 and 1.8 mm (Figure 6). We see that there is a very good correlation between testing and simulation. It is interesting to observe how the contact area evolves according to the imposed displace- ment. For low compression, the contact zone is located at the front of the disk on a thin strip. When the compression increases, the contact area grows and divides into two peripheral areas.

Finally, for larger compression, contact area as well as the maximum pressure converge to the inside edge of the disc carrier. There is also a loss of contact in the area at the front of the tool that does not appear clearly in the experimental results.

We therefore have a modeling tool of the compression of the disk that gives satisfactory results regarding the geometry of the contact area. However, the distribution of pressure field at the contact has not been confirmed yet.

EXPERIMENTAL INVESTIGATIONS

The experiment allows us to show correlations and differences obtained between simulation and geometrical deviation. The test surface is a planar face the section of which is 50 mm x 50 mm machined in a stock made of X35CrMoV5 (Hrc 53) steel.

• The plane is milled using a 12 mm diameter hemispherical tool along parallel planes so that the scallop height is inferior to 2 µm.

• Pre-polishing operation is performed with a 18 mm diameter abrasive disk made of silicon carbide (grade 600). Tool paths are trochoidal trajectories based on parallel planes (Figure 7). The distance between planes is set to 13 mm, the trochoid diameter and step are set to 8 mm and 0.8 mm respectively. The imposed displacement is 0.5 mm leading to a contact force of 5 N [START_REF] Pessoles | Automatic polishing process of plastic injection molds on a 5-axis milling center[END_REF].

The programmed feedrate (2 m/min) is reach during most of the operation. width of the polished area for a straight-line tool path. A constant and uniform pressure distribution P has been used. Thus, the result of this first simulation depicted on Figure 9 shows the mapping of the local polishing time weighted by the tool wear Tw(j) (eq. 12). We are able to distinguish the polished areas and where the tool path overlaps. It seems that the efficiency of the abrasive disk along the tool path is under estimated.

A second simulation has been computed by using the contact area issued from the finite element analysis without the abrasive paper (Figure 10). This simulation represents the mapping of Tw(j) (eq. 12). The contact area is included in the term Ins(i,j). There is a significant difference between simulation and measure. The pseudo time (s) spent on each point of the surface is half as large as in the previous simulation. The contact area seems to be thinner than the real one. We are able to distinguish each arc of the trochoid tool path.

A third simulation including the abrasive paper has been conducted (Figure 11). The abrasive paper is introduced as an additional part with a purely elastic behavior in the finite element analysis. An approximate Young modulus has been identified through a compression test of several superposed abrasive disks (E=22MPa). The thickness of the contact area is greater leading to a superior local polishing time. The abrasion topography is much more satisfying and is closer to the real part and to the elliptic simulation.

CONCLUSION AND FUTURE WORKS

In this paper we have presented a simulation tool to predict the geometrical deviation due to polishing operation with abrasive paper. The model is based on Preston equation and includes tool wear and contact area modeling. The simulation results are quite satisfactory, although they are not perfect. The introduction of FEA including the abrasive paper in the simulation provides good results. Although the simulation of the deformation of the disc carrier is satisfying, the modeling of the contact area must be improved. The tool wear is established for the same set of tool and workpiece material and for the same polishing parameters. However, the abrasion efficiency seems to be under estimated.

The simulation tool has been tested for planar polishing but we intend to extend our model on convex and concave surfaces. Then we will be able to tunes the polishing strategy (imposed displacement, feedrate and tool path) in order to manage the geometrical deviations on free form surfaces.

Figure 1 :

 1 Figure 1: Polishing tool and trochoidal tool paths.

Figure 3 :

 3 Figure 3: Simulation model.

Figure 4 :

 4 Figure 4: experimental tool wear evolution.

Figure 5 :

 5 Figure 5: Numerical simulation.

Figure 6 :

 6 Figure 6: Simulation and experiment (imposed displacements: 0.6; 1.4; 1.8mm).

Figure 7 :

 7 Figure 7: Polishing trajectory.

Figure 8 :

 8 Figure 8: Measurement.

Figure 10 :

 10 Figure 10: "FEA" simulation.

Figure 9 :

 9 Figure 9: "Elliptic" simulation.

Figure 11 :

 11 Figure 11: "FEA" with abrasive paper.

ACKNOWLEDGMENTS

This work is founded by the French Ministry of Research. It was carried out within the Manufacturing 21 working group, which comprises 16 French research laboratories. The topics covered are: modeling of the manufacturing process, virtual machining and emergence of new manufacturing methods.


B Maury

A Preux

Keywords: Pressureless Euler equations, congestion, optimal transportation

Pressureless Euler equations with maximal density constraint : a time-splitting scheme

Introduction

We are interested in the pressureless gas dynamics system, that describes the free motion of inertial particles. The system simply expresses mass conservation and momentum conservation:

∂ t ρ + ∂ x (ρu) = 0 ∂ t (ρu) + ∂ x (ρu 2 ) = 0. (1)
Here, ρ represents the density of particles and u is the velocity field. The fact that a single velocity can be locally defined is a strong implicit assumption here. The real free transport equation for non-interacting particles would be of the kinetic type (namely the Collisionless Boltzmann equation), allowing for various velocities to coexist at the same place. In the previous system [START_REF] Bouchut | On zero pressure gas dynamics[END_REF], even for smooth initial data, transport characteristics are likely to cross, leading to an incompatibility in terms of velocity. Preserving the monokinetic character of the representation calls for considering interactions between particles, although they do not explicitly appear in the equation. As described on [START_REF] Bouchut | On zero pressure gas dynamics[END_REF] or [START_REF] Brenier | Minimal geodesics on groups of volume preserving maps and generalized solutions of the Euler equations[END_REF], this model can create Dirac Masses even if initial data are smooth. It is therefore necessary to define measure-valued solutions for this system. The existence of such solutions of (1) has been proven constructively in [START_REF] Grenier | Existence globale pour le systeme des gaz sans pression[END_REF] and [START_REF] Rykov | Generalized variational principles, global weak solutions and behavior with random initial data for systems of conservation laws arising in adhesion particle dynamics[END_REF], by approximating the initial measure by Dirac masses (sticky particles). The motion of the corresponding collection of particles is then computed, following an event-driven approach to handle binary collisions, and the corresponding sticky particle solutions are shown to converge as the number of particles tends to infinity. This approach has been introduced in [START_REF] Ya | Gravitational instability: An approximate theory for large density perturbations[END_REF] as a model describing the formation of galaxies in the early stage of the universe. This system has been studied in the framework of hyperbolic systems of conservation laws (see for instance [START_REF] Brenier | Sticky particles and scalar conservation laws[END_REF] and [START_REF] Nguyen | Pressureless Euler/Euler-Poisson systems via adhesion dynamics and scalar conservation laws[END_REF]), although it lacks hyperbolicity. A proof of uniqueness has been given in [START_REF] Huang | Well Posedness for Pressureless Flow[END_REF] with appropriate conditions on entropy and energy. This system has also been formulated and analyzed as a first order differential inclusion, see [START_REF] Natile | A Wasserstein approach to the one-dimensional sticky particle system[END_REF][START_REF] Cavalletti | A Simple Proof of Global Existence for the 1d Pressureless Gas Dynamics Equations[END_REF]. Some other approaches have been proposed, for instance, in [START_REF] Cavalletti | A time discretization for the pressureless gas dynamics equations[END_REF] with an optimal transportation approach, with a viscosity regularization in [START_REF] Sobolevskii | The small viscosity method for a one-dimensional system of equations of gas dynamic type without pressure[END_REF] and [START_REF] Boudin | A solution with bounded expansion rate to the model of viscous pressureless gases[END_REF], with finite size sticky particles in [START_REF] Wolansky | Dynamics of a system of sticky particles of finite size on the line[END_REF].

We are interested here in the situation where the density is subject to remain below a threshold value, say 1; it corresponds to the pressureless gas dynamics with congestion constraint :

                 ∂ t ρ + ∂ x (ρu) = 0 ∂ t (ρu) + ∂ x (ρu 2 + π) = 0 ρ ≤ 1 (1 -ρ)π = 0 π ≥ 0.
(2) A first approach and numerical algorithm have been given in [START_REF] Bouchut | A Hierarchy of Models for Two-Phase Flows[END_REF], and the existence of solutions is proved in [START_REF] Berthelin | Existence and weak stability for a pressureless model with unilateral constraint[END_REF]. The approach is constructive, and dedicated to the one-dimensional setting: it is based on a generalization of the sticky particles model. Particles have been replaced by macroscopic congested blocks, but the dynamic is similar, with a pressure π that is active in congested zone only (a typical collision between two identical blocks is represented in Fig. 1

(top)).

In [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF] a similar system is obtained and analysed as a limit of the Aw-Rascle model.

Both systems (1) and (2) admit infinitely many solutions, some of them being highly non-physical. Consider for example the case without congestion (1), and ρ defined as follows: ρ is a Dirac mass at 0 during [0, 1], and splits at time 1 into two half Dirac masses, one going to the right at constant velocity u, and the other one at velocity -u, where u > 0 is arbitrary. This measure path, with its obvious associated velocity, verifies the system although it does not correspond to any physical reality. The fact that such solutions are compatible with (1), although no interaction force explicitly appears, can be explained: a Dirac mass can be considered as the sum of two half Dirac masses at the same point. Splitting is achieved by exerting an impulsion of one of the halves, and the opposite impulsion on the other. Since these impulsions sum up to zero, it can be done without any real action on the system, and the obtained path verifies the weak formulation of the system. The splitting time is of course arbitrary, and each of the subparticles can be further split in lighter particles, at any time, as soon as local conservation of mass and momentum are respected. Similar bizarre solutions can be built for the congested case. An initially steady block can be suddenly splits in two halves, with a pressure field that is singular in time, but regular in space at the instant of the splitting (piecewise affine function vanishing at the ends of the block, maximal at the center). One may hope to recover uniqueness by demanding that the kinetic energy decreases. This is not enough, as the following example shows: consider two Dirac masses heading to 0 with opposite velocities u and -u. They may collide at 0 into a steady Dirac mass with weight 2, but they may also bounce against each other, and head back to -∞ and +∞ at velocities -eu and eu, respectively. For any e ∈ [0, 1] (restitution coefficient), the solution is energy decreasing (with conservation for the purely elastic case e = 1).

To sum up, both systems allow for multiple solutions to exist, some of them being non-physical (energy increasing). But still infinitely many energydecreasing solutions exist. A collision law is obviously lacking. Before writing the multidimensional version of (2), let us describe now the two main ingredients that we will use to write the model (with a collision law) and to elaborate a time-splitting strategy.

First ingredient: microscopic granular model

At the microscopic level, the motion of rigid particles can be modelled as follows. Consider N rigid spheres in R d , with common radius r > 0, unit mass, and centers q 1 , . . . ,q N . Denote by q ∈ R dN the vector of positions, and let D ij = |q j -q i | -r denote the inter-grain distance. Prescribing grain rigidity amounts to define a set of feasible configurations, that is

K = q = (q 1 , . . . , q N ) ∈ R dN , D ij ≥ 0 ∀i = j .
The cone of feasible velocities (i.e. velocities that do not lead to overlapping between grains) is defined for any q ∈ K as

C K (q) = {v ∈ R dN , D ij = 0 =⇒ G ij • v ≥ 0},
where G ij ∈ R dN is the gradient of D ij . The inertial motion of this system of grains submitted to a force field f (possibly including interactions, i.e. f may depend on q) with nonelastic collision can be written (see e.g. [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF][START_REF] Maury | A time-stepping scheme for inelastic collisions[END_REF])

                           dq/dt = u du/dt = f + i<j π ij G ij q ∈ K π ij D ij = 0 ∀i, j π ij ≥ 0 ∀i, j u + = P C K (q) (u -). (3)
The π ij 's play the role of a discrete pressure, that is a Lagrange multiplier associated to the non-overlapping constraint q ∈ K. Velocities u + and u -are pre-and post-collisional velocities (as detailed in [START_REF] Maury | A time-stepping scheme for inelastic collisions[END_REF], the velocity has bounded variation, so that left and right limits are well-defined). Note that, when there is no contact, all the distances D ij are positive, so that π ij ≡ 0, C K (q) = R dN , and we recover the usual Newton law for a N-particle system. The interaction forces π ij are actually measures in time, and they are indeed singular as soon as a collision happens. Existence of solutions to this system is proved in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF][START_REF] Maury | A time-stepping scheme for inelastic collisions[END_REF]. Note that, although the system is formally well-posed (uniqueness holds when the number of collision is bounded), uniqueness requires analiticity of the forcing term f (see in particular [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF]).

Remark 1.

Let us remark that alternative collisions laws may be accounted for (we shall not extend this generalization to the macroscopic setting) by introducing the so called outward normal cone to K, that is

N K (q) = C • q = v ∈ R dN , v • w ≤ 0 ∀w ∈ C K (q) =    v = - i<j π ij G ij , π ij ≥ 0 , D ij π ij = 0    .
Given a restitution coefficient e ∈ [0, 1] (e = 0 for non elastic shocks and e = 1 for purely elastic shocks), a more general collision law is (see e.g. [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF])

u + = u --(1 + e)P N K (q) (u -). (4)

Second ingredient: first order crowd motion model

The handling of congestion in evolution equations has been studied for the macroscopic crowd motion introduced in [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF] and [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF]. The model is based on a so-called desired velocity, that is the velocity U (x) someone at x would spontaneously have, if he/she were alone. We aim at accounting for a congestion constraint ρ ≤ 1, we denote the corresponding set of densities by K. To that purpose, we consider that the desired velocity is projected (in a L 2 sense) on the cone of feasible velocities. The latter set is defined as the set of all those velocities that do not lead to a violation of the constraint. More formally (see [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] for details), it is defined in a dual way :

C K (ρ) = v ∈ L 2 , v • ∇q ≤ 0 , ∀q ∈ H 1 , q ≥ 0 , q(1 -ρ) = 0 . ( 5
)
The model writes

∂ t ρ + ∇ • (ρu) = 0 u = P CK (ρ) (U ). ( 6
)
In the spirit of the Catching Up algorithm initially proposed in [START_REF] Moreau | Evolution problem associated with a moving convex set in a hilbert space[END_REF], and applied to a microscopic crowd motion model in [START_REF] Maury | A discrete contact model for crowd motion[END_REF], a time-stepping algorithm has been proposed in [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF] to handle this type of evolution problem. It simply consists in transporting the density during one time step with the desired velocity, possibly violating the constraint, and then project in the Wasserstein sense in the set K of admissible densities. It reads as follows :

ρn+1 = (Id + τ U ) # ρ n ρ n+1 = P W2 K (ρ n+1 ). (7)
Now let us expand System (2) in higher space dimension, with the continuous counterpart of a non-elastic collision law:

                     ∂ t ρ + ∇ • (ρu) = 0 ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇π = 0 ρ ≤ 1 (1 -ρ)π = 0 π ≥ 0 u + = P CK (ρ) (u -). ( 8
)
In the previous system, π is a nonnegative scalar field that is the continuous counterpart of the interaction pressures π ij ′ s accounting for interactions between neighboring grains in the microscopic model. As previously said, in the onedimensional setting, it is natural and easy to build exact solutions to this system. Let us stress that sticky block solutions described in [START_REF] Berthelin | Existence and weak stability for a pressureless model with unilateral constraint[END_REF] verify the non-elastic collision law that we added to the system. Such a collision is represented in Fig. 1 (top). See in particular the pressure field π (middle) that is activated at the very instant of impact; its opposite gradient instantaneously anihilates the pre-collisional velocity field. In higher dimension, such analytical solutions are not available. Consider for instance the two-dimensional counterpart of the 1D sticky blocks: two rectangular blocks colliding with opposite velocities. At the instant of collision, it can be checked that the pressure field π solves a Poisson equation over the two-blocks domain, with a right-hand-side that is the opposite of the divergence of the precollisional velocity (it is a singular distribution supported by the interface), and with homogeneous Dirichlet boundary conditions on the exterior boundary. The corresponding pressure field and the corrected velocity are represented in Fig. 1 (bottom), on the middle and on the right, respectively.

In the one-dimensional setting, existence of solutions can be obtained in various ways, as already mentioned ([15, 16, 17]). In higher dimension, system (8) (without the non-elastic collision law) can be obtained as a formal limit of compressible Euler equations with a singular pressure term (see [START_REF] Degond | Numerical simulations of the Euler system with congestion constraint[END_REF]). Note also that solutions to a viscous version of (1) can be built (see [START_REF] Lions | On a free boundary barotropic model[END_REF][START_REF] Bresch | Singular limit of a Navier-Stokes system leading to a free/congested zones two-phase model[END_REF][START_REF] Perrin | Free/Congested Two-Phase Model from Weak Solutions to Multi-Dimensional Compressible Navier-Stokes Equations[END_REF]) as limits of solutions to compressible Navier-Stokes equations with a stiff pressure term.

Yet, up to our knowledge, there is no well-posedness result (in terms of existence or uniqueness) for the congested Euler system (8), for dimensions higher than 1.

We propose here a new time-discretization scheme for [START_REF] Huang | Well Posedness for Pressureless Flow[END_REF] in the spirit of the prediction-correction scheme used in the crowd motion case. This scheme will allow us to build sequences of densities and velocities, that may be expected to converge to weak solutions of our problem. Let us stress that we are not able to establish such a convergence, but we provide partial results that assert a certain consistency of the scheme with respect to the equations (see Section 2.2). Beside, a numerical scheme can be built on this time-splitting strategy, and we shall present preliminary results in Section 3.

Compared to the first order in time model, the present context of inertial particles raises two deep additional difficulties: 1. The velocity field is no longer prescribed as in the crowd motion case, it is attached to the mass and it is highly varying in time. In particular, it is not possible to control its regularity in space, so that the transport step in [START_REF] Nguyen | Pressureless Euler/Euler-Poisson systems via adhesion dynamics and scalar conservation laws[END_REF] may not be one-to-one, which means that particle trajectories cross each other. This is likely to lead to non-physical behaviour: particles crossing without seeing each other, or even artefactual energy increases during the subsequent projection phase. This calls for an additional step in the time-stepping procedure, by prescribing some sort of monotonicity in the transport phase. 2. For the very same reason (the velocity field is not prescribed once and for all), velocities have to be updated at the end of each time step. In particular, the effect of congestion (i.e. the projection phase) changes the velocity field and reduces energy, since we consider non-elastic collisions. This necessitates to follow the motion in a Lagrangian way during the whole scheme, i.e. to map the initial density to the final one respectfully of the actual motion of particles.

Time-stepping scheme

We first present here the issues raised by the attempt to propose a relevant time-stepping strategy for the congested pressureless Euler equations. We shall in particular explain, in a one-dimensional setting, how the two constraints (particles do not cross and the density remains below 1) can be accounted for in a mutually consistent way. We then present the time-stepping scheme in any space dimension, and we detail some of its properties.

Time discretization strategy

We first consider the microscopic model (3) in one dimension (aligned grains), assuming that the grain diameter is zero (which makes clear sense in the onedimensional setting), and that no force is exerted. We assume that, at initial time, the ordering is consistent with the indexing:

q i (0) ≤ q j (0) for 0 ≤ i ≤ j ≤ N.
Since they do not leap over each other, the set of admissible configurations expresses that the order is preserved:

K = {q ∈ R N , q i ≤ q j ∀i ≤ j},
and the set of feasible velocities follows

C q = {v ∈ R N , v i ≤ v j whenever q i = q j }.
The model of 1d granular flow with inelastic collisions may be phrased formally

dq dt = u , u + = P Cq (u -).
For this model, a numerical scheme and convergence results have been obtained in [START_REF] Bernicot | Existence results for non smooth second-order differential inclusions, convergence result for a numerical scheme and application to the modeling of inelastic collisions[END_REF] and [START_REF] Maury | A time-stepping scheme for inelastic collisions[END_REF]. Let τ > 0 be the time step and consider q n , v n ∈ R n the position and velocity vectors at time nτ . We ensure that the updated position vector q n+1 = q n + τ u n is admissible by prescribing ∀i ≤ j, q n i + τ u n i ≤ q n j + τ u n j and we call C τ q n the corresponding set of velocities. The approximated solutions are built using the following scheme :

u n+1 = P C τ q n (u n ); q n+1 = q n + τ u n+1 . ( 9
)
Note that the conditions defining C τ q n can be rewritten

u n i -u n j • q n i -q n j ≥ -1 τ q n i -q n j 2 . ( 10
)
There is an obvious link between this microscopic setting and the macroscopic system [START_REF] Bouchut | On zero pressure gas dynamics[END_REF]. Indeed, to any solution (q(t), u(t)) of the previous model, one can associate a solution to the macroscopic equation ( 1) considering the density

ρ t = n i=1 δ qi(t) and ρ t u t = n i=1 u i (t)δ qi(t)
. It suggests a time-stepping scheme to tackle the unconstrained macroscopic problem in one dimension, by simply requiring that trajectories do not cross during a time step τ , in the spirit of the strategy that is proposed in [START_REF] Cavalletti | A time discretization for the pressureless gas dynamics equations[END_REF]. This scheme reads as follows: consider the density-velocity couple (ρ n , u n ) at time step n, define the velocity ũn as the projection of u n on the cone of feasible velocities (macroscopic counterpart of ( 10))

C τ ρ n = v ∈ L 2 ρ n , (x -y)(v(x) -v(y)) ≥ -1 τ |x -y| 2 for ρ n ⊗ ρ n -a.e.(x, y) ,
and move ρ n by Id + τ ũn to obtain ρ n+1 . The new velocity u n+1 that is defined on the support of ρ n+1 is finally obtained by conservation of momentum between any y in the support of ρ n+1 and the set of its preimages x ∈ (Id + τ ũn ) -1 . Now consider the problem with congestion, i.e. the density is subject to remain below 1, it is natural to treat it in the spirit of Algorithm [START_REF] Nguyen | Pressureless Euler/Euler-Poisson systems via adhesion dynamics and scalar conservation laws[END_REF], that was dedicated to the first order in time model (crowd motion), by considering that the previous velocity u n plays the role of the desired velocity. It consists in projecting the predicted measure (obtained according to the previous considerations) on the set of feasible densities (second step of ( 7)).

The drawback of this scheme is that the intermediate density obtained after handling the non-crossing of trajectories may contain a singular part, typically at points y that can be written x+ τ ũn (x) for x running over a set with non zero ρ n -measure. This feature would not be a problem if the pre-collisional velocities were prescribed, like in the crowd motion problem. But in the case of pressureless equation, it is necessary to build the new velocity field over the support of the new density. It requires to follow the motion of individual particles before and after the two sub-steps of the procedure. This is prevented by the presence of Dirac masses in the intermediate measure: they make it impossible to properly define a map between the previous density field and the next one, in a way that would respect the Lagrangian motion of particles.

For this reason, we propose to reinforce the non-crossing constraint by prescribing that trajectories do not cross in a strict sense, i.e. particles remain, during the whole time step, at a distance that is bounded from below. More precisely, consider two particles initially at x and y, with velocities v x and v y , respectively. We require that the moving particles remain at a distance that is larger than c times their initial distance, where c ∈ (0, 1) is a small parameter:

x -y, (x + τ v x ) -(y + τ v y ) ≥ c |x -y| 2 ⇐⇒ x -y, v x -v y ≥ c -1 τ |x -y| 2 .
This constraint will be integrated in the process through the definition of the set of feasible velocities introduced in the next section (Eq. ( 12)). The positivity of c forces the intermediate density ρn+1 to remain bounded (see Lemma 2.2) and the transport Id + τ ṽn to be invertible. Thus, one can reconstruct the trajectory of particles from the density ρ n+1 to ρ n with the composition of two Lagrangian transports.

The next section will detail the time-stepping scheme issued from the previous one-dimensional consideration. The approach we propose is far from being covered by a full convergence analysis. Yet, beyond the fact that very few methods exist to build solutions for the congested pressureless Euler equation in dimension d ≥ 2, the approach is supported by the following properties, that are detailed in the next section: 1. The scheme is well-defined, in the sense that each time step defines in a unique way a new density, that verifies the congestion constraint, and an associated velocity field; 2. In a one-dimensional setting, in the case of an initial density that is a collection of saturated clusters (sum of characteristic functions of intervals), our scheme essentially behaves like the so-called sticky block approach, meaning that numerical solutions actually converge, in this setting, to solutions of the continuous problem (as proved in [START_REF] Berthelin | Existence and weak stability for a pressureless model with unilateral constraint[END_REF]); 3. In any dimension, the scheme conserves the momentum of material systems during collision. More precisely, for each material subsystem involved in a numerical collision, the momentum is conserved; 4. Since the approach is not based on a relaxation of the congestion constraint, it natively handles empty zones without difficulty, as confirmed by the numerical illustrations that are presented.

The scheme

Let Ω ⊂ R d be a bounded domain, we denote by P( Ω) the space of probability measures on Ω. We shall assume in what follows that the domain is chosen sufficiently large, so that all considered densities are supported in it.

Let τ > 0 be a time step, ρ n ∈ P( Ω) and u n ∈ L 2 ρ n ( Ω) represent the density and velocity at time nτ . The first step consists in projecting the current velocity field on the set of M c ρ n , to prevent particles from crossing. The density is then moved according to this projected velocity field, during τ . The obtained intermediate density is likely to violate the constraint. The third step consists in projecting it (in the Wasserstein sense), on the feasible set K, to account for the saturation constraint. The final step consists in defining the new (post-collision) velocity, using the Lagrange map between the previous density and the new one.

For (ρ n , u n ) given, the scheme reads:

ρ n ρn+1 ρ n+1 s n+1 r n+1 (s n+1 ) -1 • r n+1 K                ũn = P M c ρ n (u n ) ρn+1 = (Id + τ ũn ) # ρ n ρ n+1 = P W2 K (ρ n+1 ) u n+1 = Id -(s n+1 ) -1 • r n+1 τ (11)
with [START_REF] Sobolevskii | The small viscosity method for a one-dimensional system of equations of gas dynamic type without pressure[END_REF] and c ∈ (0, 1) a small parameter. K is the space of all those probability measures supported in Ω that admit a density less than 1. The mappings involved in the last step (reconstruction of the velocity) are

M c ρ n = {v ∈ L 2 ρ n ( Ω) , v(x) -v(y), x -y ≥ c -1 τ |x-y| 2 ρ n ⊗ ρ n -a.e.},
s n+1 = Id + τ ũn ,
and r n+1 , that is the optimal transportation map between ρ n+1 and ρn+1 . Figure 2 illustrates how the time-stepping scheme treats a one-dimensionnal collision between two blocks. The figure reads from bottom to top. The first part of the scheme moves densities with the current velocity, while avoiding contact of pathlines (the velocity is projected on M c ρ n ). This intermediate density violates the congestion constraint, and the second step is meant to enforce it, by projection of the density on K. The 3 circles on the figure correspond to successive positions of a physical particle during the scheme: initially at x, it is moved to y, and then to z by projection. For this very particle, the reconstruction of the after-collision velocity (at the new position z) gives (z-x)/τ . As explained in the proof of Proposition 3, it may happen that a gap remains between the two blocks after this first partial treatment, in that case it may need two more time steps to fully handle the collision.

Let us start by establishing some basic properties pertaining to the proposed scheme. The first lemma asserts that the first step (projection on M c ρ n ) is welldefined, and preserves the total momentum. In a second lemma, we show that the intermediate density ρn+1 is bounded. In the third one, we recall a known result concerning the projection on K.

u -v L 2 ρ ( Ω) ≥ u -ũ L 2 ρ ( Ω) , ∀v ∈ L 2 ρ ( Ω).
Moreover, we have the following properties :

(i) ũ L 2 ρ ( Ω) ≤ u L 2 ρ ( Ω) ; (ii) Ω ũ(x)dρ(x) = Ω u(x)dρ(x).
Proof. Since M c ρ is a closed and convex set, the projection ũ := P M c ρ (u) is well-defined. Inequality (i) comes from the fact M c ρ is also a cone pointed at 0. Let us prove now that the first step of the scheme does not change momentum (identity (ii

)). If u ∈ M c ρ , then u + k is also in M c ρ for any k in R d . If ũ is
the projection of u on M c ρ , we have :

u -ũ 2 L 2 ρ ( Ω) ≤ u -ũ -k 2 L 2 ρ ( Ω)
and thus :

u -ũ, k |k| L 2 ρ ( Ω) ≤ |k| 2 .
Letting |k| tend to 0 yields u -ũ,

k |k| L 2 ρ ( Ω)
≤ 0, for all k ∈ R d and consequently

u -ũ, k |k| L 2 ρ ( Ω) = 0.
Consider the total momenta before and after the projection of velocities:

E tot = Ω u(x)dρ(x) and Ẽtot = Ω ũ(x)dρ(x).
We have

E tot , k = Ẽtot , k for all k ∈ R d and finally Ẽtot = E tot . Lemma 2.2. Let ρ := s # ρ = (Id + τ ũ) # ρ, with ρ ∈ L ∞ ( Ω), and ũ ∈ M c ρ , we have : ρ L ∞ ( Ω) ≤ 1 c d ρ L ∞ ( Ω) .
Proof. We consider here that ũ is regular, and we refer to [START_REF] Ambrosio | Gradient flows: in metric spaces and in the space of probability measures[END_REF], Section 5.5, for the general case. Since ũ is in M c ρ , the eigenvalues of ∇s are larger than c and then det ∇s is larger than

c d . Finally, since ρ = ρ | det ∇s| • s -1 we have ρ L ∞ ( Ω) = ρ | det ∇s| • s -1 L ∞ ( Ω) ≤ 1 c d ρ • s -1 L ∞ ( Ω) ≤ 1 c d ρ L ∞ ( Ω) ,
which ends the proof.

The following lemma concerns the projection on K. It is proven in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF].

Lemma 2.3. Let ρ ∈ P( Ω).

There exists a unique ρ K ∈ K such that

W 2 (ρ, µ) ≥ W 2 (ρ, ρ K ), ∀µ ∈ K,
and there exists π

∈ {p ∈ H 1 | p ≥ 0, (1 -ρ)p = 0} such that ρ = (Id + ∇π) # ρ K .
We may know prove that the scheme defines a unique density-velocity couple.

Proposition 1. Let (ρ n , u n ) ∈ K × L 2 ρ n ( Ω). The scheme (11) defines a unique (ρ n+1 , u n+1 ) ∈ K × L 2 ρ n+1 ( Ω).
Proof. We can check that s n+1 is strictly monotone and consequently injective.

We can define an inverse on s n+1 (spt(ρ n )) = spt(ρ n+1 ). Thanks to Lemma 2.3, we have existence and uniqueness of r n+1 : spt(ρ n+1 ) → spt(ρ n+1 ) and consequently the uniqueness of u n+1 and ρ n+1 . Moreover, we have

Ω (u n+1 (x)) 2 ρ n+1 (x)dx = Ω x -(s n+1 ) -1 • r n+1 (x) τ 2 ρ n+1 (x)dx = Ω (r n+1 ) -1 (x) -(s n+1 ) -1 (x) τ 2 ρn+1 (x)dx ≤ 2 Ω (r n+1 ) -1 (x) -x τ 2 + x -(s n+1 ) -1 (x) τ 2 ρn+1 (x)dx ≤ 2 Ω x -r n+1 (x) τ 2 ρ n+1 (x)dx + 2 Ω s n+1 (x) -x τ 2 ρ n (x)dx ≤ 2 Ω (∇π n+1 (x)) 2 ρ n+1 (x)dx + 2 Ω ũn (x) 2 ρ n (x)dx ≤ 4 Ω u n (x) 2 ρ n (x)dx since Ω ∇π n+1 (x) 2 ρ n+1 (x)dx = W 2 (ρ n+1 , ρn+1 ) ≤ W 2 (ρ n , ρn+1 )
and

W 2 (ρ n , ρn+1 ) ≤ Ω u n (x) 2 ρ n (x)dx,
which ends the proof.

Proposition 2. Total mass and momentum are conserved by the scheme [START_REF] Cavalletti | A time discretization for the pressureless gas dynamics equations[END_REF].

Proof. Conservation of mass is inherent to the approach, which is a succession of transports and projections in the Wasserstein sense. As for the momentum, we have already seen that Ω E n (x)dx = Ω Ẽn (x)dx, where Ẽn = ρ n ṽn .

Ω Ẽn (x)φ(x)dx = Ω ρ n (x)ũ n (x)φ(x)dx = Ω ((s n+1 ) -1 • r n+1 ) # ρ n+1 (x)ũ n (x)φ(x)dx = Ω ((s n+1 ) -1 • r n+1 ) # ρ n+1 (x)( s n+1 (x) -x τ )φ(x)dx = Ω ρ n+1 (x)( r n+1 -(s n+1 ) -1 • r n+1 τ )φ((s n+1 ) -1 • r n+1 (x))dx = Ω ρ n+1 (x)(∇π n+1 + u n+1 )φ((s n+1 ) -1 • r n+1 (x))dx So Ẽn τ = (s n+1 ) -1 • r n+1 ) # (ρ n+1 ∇π n+1 + E n+1
) and consequently :

Ω Ẽn (x)dx = Ω ρ n+1 (x)(∇π n+1 (x) + u n+1 (x))dx = Ω ρ n+1 (x)∇π n+1 (x) + Ω ρ n+1 (x)u n+1 (x)dx = Ω E n+1 (x)dx.

Remark 2. A direct consequence of the previous proposition is that the center of mass moves at a constant velocity

U = Ω u n (x)ρ n (x)dx Ω ρ n (x)dx
.

The next proposition pertains to the one-dimensional setting. We consider the exact solution that corresponds to the collision of two saturated blocks. We establish that, except for the few times steps (no more than 3) that are needed to handle the collision at the discrete level, the discretization scheme recovers the exact solution.

Proposition 3. Let ρ t be the density representing two blocks colliding at time t 0 and η > 0. If c ≤ √ 2 -1 then there exists τ o > 0 such that for all time step τ ≤ τ o the solution given by the scheme is exact, for all step n ∈ {n ∈ N, nτ / ∈ [t 0 -η, t 0 + η]}.

Proof. Let us start by two basic remarks. Firstly, from remark 2, the scheme computes in an exact way the motion of the mass center. Secondly, an isolated block moving at uniform velocity is preserved by the scheme (blocks stay blocks). As a consequence, we just have to check that the scheme makes the blocks stick together in a finite number of steps (no more than 3, as we shall see).

Let ∆U n ≥ 0 be the velocity difference between the two blocks at time nτ and e n their distance. One can check that the updated velocity difference ∆U n+1 is equal to e n -e n+1 τ . The velocity field u n remains in M c ρ n if and only if the constraint is verified in the ends of blocks, so that

∆U n (-e n ) ≥ c-1 τ (e n ) 2 ⇔ τ ∆U n ≤ (1 -c)e n
Now, if τ ∆U n > (1 -c)e n and if the length of the smaller block L is sufficiently large (take τ such that L ≥ 2τ ||v|| ∞ ) the density ρn+1 exceeds the congestion constraint (see Figure 2) over two blocks of total size cτ ∆U n 1-c -ce n . In this case, the mass M projected in the gap between the two blocks is larger than the half of over-congested mass,

M ≥ 1 2 1 c -1 cτ ∆U n 1 -c -ce n ≥ 1 2 (τ ∆U n + (c -1)e n ).
Two cases have to be considered. If τ ∆U n ≥ (1 + c)e n , the mass M is larger than the mass necessary to refill the gap, this implies e n+1 = 0 and τ ∆U n+1 = e n . The following step allows blocks to have the same velocity, e n+2 = 0 and ∆U n+2 = 0. So, the scheme handles the collision with 2 steps. Now if τ ∆U n ∈ ((1 -c)e n , (1 + c)e n ), the projection on the set of feasible velocities imposes that ẽn+1 = ce n (the distance between the ends of blocks is exactly c times the previous one), then e n+1 ≤ ce n and consequently

τ ∆U n+1 = e n -e n+1 ≥ (1 -c)e n .
Thus,

(1 + c)e n+1 ≤ (1 + c)ce n ≤ (1 -c)e n ≤ τ ∆U n+1
since c ≤ √ 2 -1, conditions of the first case are obtained, the scheme handles the collision with 3 steps.

Finally, taking τ

o = min{ L 2 v ∞ , 2η
5 } allows us to prove the property.

Numerical illustrations

We briefly describe here a space discretization strategy that can be carried out to perfom actual simulations, and we present some numerical tests to illustrate the behavior of the numerical scheme in various situations.

Space discretization scheme

We take Ω = [0, 1] 2 and we discretize the domain in a collection of cells (C i ) i∈I with a uniform size h × h and the center located at x i . The quantities used previously are approximated by constant functions in each cell :

f ≈ i∈I f (x i )1 Ci := i∈I f i 1 Ci ,
and we denote by X α h the corresponding space (α = 1 for a scalar field, α = 2 for a vector field). From now on, we will assimilate f to the vector (f i ) ∈ (R d ) I .We denote by I ρ the set of indices corresponding to non-zero values of the ρ i 's:

I ρ = {i ∈ I | ρ i > 0}.
We may now describe the 4 steps of the numerical scheme.

Suppose that we have (ρ

n , v n ) ∈ X 1 h × X 2 h . Step 1 : We replace M c n by the space of functions v in X 2 h such that ∀i, j ∈ I ρ n , v i -v j , x i -x j ≥ c -1 τ |x i -x j | 2 .
The projection can be formulated in a saddle-point form: Consider the function

F ij := v → v i -v j , x i -x j -c-1 τ |x i -x j | 2
, the projection onto this space is defined as ṽn = argmin

ṽ∈R I    sup λij ∈R I 2 +   ṽ -v n ℓ 2 ρ n - i,j∈(I ρ n ) 2 λ ij F ij (v n )      .
Step 2 : We approximate the transport step by a projected Lagrangian scheme (see Remark 3):

ρn+1 := (Id + τ ṽn ) # ρ n ≈ i∈I   Ci (Id + τ ṽn ) # ρ n (x)dx   1 Ci , Ẽn+1 := (Id + τ ṽn ) # E n ≈ i∈I   Ci (Id + τ ṽn ) # E n (x)dx   1 Ci .
Step 3 : In order to project on K, we define the discrete transport plan as the minimizer of a discrete counterpart of the transport cost, among all those plans ζ = (ζ ij ) i,j∈I with ρn+1 as first marginal, and second marginal in K. The cost is defined as

C(ζ) = i∈I ρn+1 j∈I β ij ,
where the cell-to-cell costs β ij are defined (see Remark 4) by

β ij =      ζ 2 ij ρn+1 i |x i -x j | 2 if C i is a neighbor of C j + ∞ otherwise (13)
where cells are considered neighbors if they share an edge or a vertex. The infinite cost for non-neighboring cells prevents the transport from leaping over adjacent cells. This choice is consistent with the fact that very small time steps are used in the presented illustrations.

Step 4 : We reconstruct the velocity by

E n+1 j = i∈I ρn+1 ζ ij Ẽn+1 i ρn+1 i + x j -x i τ , u n+1 j = E n+1 j /ρ n+1 j . Remark 3. Note that the space X h is not stable by the transportation step. If ρ ∈ X 1 h , then (Id+ τ v) # ρ is a piecewise constant function that is not, in general, in X 1 h . A function in X 1 h is obtained by ℓ 2 projection.

Remark 4. The step 3 is undoubtedly the most problematic one. A natural counterpart of the continuous projection problem would consist in minimizing

i,j ζ ij |x i -x j | 2
over all those discrete transport plans (ζ ij ) with ρn+1 as first marginal, and such that the second marginal belongs to K. Yet, this straight discretization of the transport cost does not lead to a satisfying approximation of the projection. This feature is related to the numerical diffusion that is inherent to the space discretization of the Finite Volume type that we chose. This delicate issue will be addressed in a forthcoming paper. Let us simply say here that, if one considers the JKO scheme applied to the gradient flow that corresponds to transport at constant velocity, such a discretization leads to an untenable behaviour. For small time steps, the discretized JKO scheme is static, i.e. the density is not transported. It comes from the fact that, because of numerical diffusion, the expression above highly overestimates the cost of the displacement from a cell to its neighbor, so that immobility is cheaper than cell-to-cell motion. This locking phenomenon can be avoided by changing the cost for short displacements. More precisely, it can be checked in the transport case that changing the cost into

ζ 2 ij ρn+1 i |x i -x j | 2
makes it possible to recover the classical upwind scheme in the one-dimensional case. This explains the choice of the cost of transport between neighboring cells that appears in (13).

Numerical tests

In this part, we illustrate our comments by comparing numerical solutions computed by the SCoPI software for the granular model (see [START_REF] Maury | A time-stepping scheme for inelastic collisions[END_REF][START_REF] Faure | Crowd motion from the granular standpoint[END_REF][START_REF] Lefebvre | Numerical simulation of gluey particles[END_REF] for details on the algorithm) with the solutions computed by our algorithm for the scheme described above. This comparison is not intended to rigorously validate the approach in any way, since the two models are different (see next section), but it asserts a satisfactory behaviour of the macroscopic model to mimic the motion of many-body suspensions. In the microscopic examples we take the number of particles (with a radius of r given below) that corresponds to the density of the macroscopic model.

First numerical example :

In this case, a small cloud of particles (on the left) is set at uniform (rightward) velocity toward another cloud of particles initially at rest. Macroscopic model :

• Ω = [-40, 40] 2 , h = 0.3, τ = 0.0150;

• ρ 0 = 0.21 x 2 +y 2 <25 2 + 0.21 (x+31) 2 +y 2 <5 2 ;

Figures 3 and4 present the evolution of the macroscopic density (left), and the suspension of rigid grains (right), for the first test case. Figures 5 and6 correspond to the second test case. In both examples, the shapes of the saturated zone that is induced by the collision of clouds present similar patterns. Note that symmetry-breaking seems to occur in the microscopic case. It is due to the fact that the grains are randomly distributed on the support of each cloud, so that the initial condition itself is not symmetric.

Conclusive remarks

Let us first stress that the macroscopic model that we propose cannot be obtained as a limit, even formally, of the hard sphere microscopic model. We refer to [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF] for a detailed discussion on these delicate micro-macro issues. The microscopic granular arrangements of rigid spheres (or discs, for the two-dimensionnal setting) induce non local effects that cannot be captured by the crude macroscopic description that we proposed. Even the very notion of maximal density does not have a clear meaning in the microscopic setting, especially in two dimensions: it depends on the local structure of grains.

Let us also mention that the macroscopic model can be extended, at least formally, to more general collision laws (see Remark 1 for the microscopic setting), i.e. with a restitution coefficient e ∈ (0, 1]. It is a straight translation of the microscopic collision (4). Indeed, the dual cone of C K (ρ) (defined by ( 5)) is the outward normal cone N K (ρ) = ∇q , q ∈ H 1 , q ≥ 0 , q(1 -ρ) = 0 which leads to the following macroscopic collision law u + = u --(1 + e)P NK (ρ) u -.

The projection on N K (ρ) is obviously well-defined, since it is a closed convex cone. Yet, although the model seems to make clear sense, and at least formally determines the evolution of the system, very deep stability issues can be expected, even in the one-dimensional setting. Let us consider the case of two identical blocks colliding with opposite velocities, in 1d (see Fig. 1, top). With e = 1 (pure elastic collision), the expected solution can be recovered: the aftercollisional velocities are simply the opposite velocities (pure bouncing between blocks). Now consider that the right hand side block, considered as 1 (0,1) is replaced by the sum of 1 (0,a) and 1 (a+ε,1+ε) , for a fixed a ∈ (0, 1). The behavior can be computed explicitely, and it is highly different from the previous situation (the sub-block at the right hand side will be pushed away at a higher velocity), for arbitrary small values of ε, which leaves little hope that the evolution may depend smoothly on initial conditions.

Fig. 1 .

 1 Fig. 1. Collision of two rectangular blocks: pre-collision velocity, instantaneous pressure and post-collision velocity (top: one-dimensional setting; bottom: two-dimensional setting)

Fig. 2 .Lemma 2 . 1 .

 221 Fig. 2. Time stepping scheme for a 1d collision

Fig. 3 .

 3 Fig.3. First example at times 0, 1, 2 and 5.

Fig. 4 .

 4 Fig. 4. First example at times 10 and 12.

Fig. 5 .

 5 Fig.5. Second example at times 0 and 1.

Fig. 6 .

 6 Fig. 6. Second example at times 1.5 and 2.

• E 0 = 21 (x+31) 2 +y 2 <5 2 .

Microscopic model :

• Ω = [-40, 40] 2 , τ = 0.03, r = 0.1; • 500 particules in {(x+31) 2 +y 2 < 5 2 } at velocity (10, 0) and 12500 particules in {x 2 + y 2 < 25 2 } at velocity (0, 0).

Second numerical example :

In this second situations, two square clouds of particles are thrown against each other. Macroscopic model :