Matthieu Constant
email: matthieu.constant@uclouvain.be

Nathan Coppin
email: nathan.coppin@uclouvain.be

Frédéric Dubois
email: frederic.dubois@umontpellier.fr

Riccardo Artoni
email: riccardo.artoni@univ-eiffel.fr

Jonathan Lambrechts
email: jonathan.lambrechts@uclouvain.be

Vincent Legat
email: vincent.legat@uclouvain.be

Numerical investigation of the density sorting of grains using water jigging

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Immersed granular flows are intensively used in industrial processes and constitute a widespread research topic. The interaction of a fluid with solid grains gives rise to complex systems for which understanding the effect of microscopic properties on the macroscopic behaviour is still limited. Sorting of grains and particles is a common challenge across many research fields yet executed in very different ways depending on the application. In chemistry and biology, classification of cells ( [START_REF] Yamada | Separation of the dead cell fraction from x-irradiated rat thymocyte suspensions by density gradient centrifugation[END_REF]) or bacteria ( [START_REF] Chesnot | Quantitative and qualitative comparison of density-based purification methods for detection of cryptosporidium oocysts in turbid environmental matrices[END_REF]) requires advanced separation techniques adapted to the microscopic size of the particles. Based on a microfluidic system, pinched-flow fractionation ( [START_REF] Yamada | Pinched flow fractionation: continuous size separation of particles utilizing a laminar flow profile in a pinched microchannel[END_REF]) used with an appropriate sedimentation force ( [START_REF] Morijiri | Sedimentation pinched-flow fractionation for size-and density-based particle sorting in microchannels[END_REF]) constitute a passive method of size and density sorting while active methods make use of an external force field like surface acoustic waves ( [START_REF] Jo | Active density-based separation using standing surface acoustic waves[END_REF]). In sedimentology, [START_REF] Steidtmann | Size-density sorting of sand-size spheres during deposition from bedload transport and implications concerning hydraulic equivalence[END_REF] carried out experiments to evaluate the effects of grain properties on the hydraulic sorting in sediment flows. Size and density sorting mechanisms have been classified in four categories ( [START_REF] Slingerland | Occurrence and formation of water-laid placers[END_REF][START_REF] Hughes | Hydraulic sorting of heavy-mineral grains by swash on a medium-sand beach[END_REF]) and affect the sediment transport in braided rivers ( [START_REF] Bridge | The interaction between channel geometry, water flow, sediment transport and deposition in braided rivers[END_REF]). In civil engineering, the environmental impact of the construction and demolition waste is a major issue ( [START_REF] Neto | An economic analysis of the processing technologies in cdw recycling platforms[END_REF]). Aggregates produced by the construction industry through demolition or maintenance are not easily reused. Recycling waste materials can be enhanced by sorting and separation systems able to improve the homogeneity of aggregate collections having the desired requirements ( [START_REF] Paranhos | A sorting method to value recycled concrete[END_REF]). Among the sorting processes, jigging is a mineral processing technique relying on gravity to separate different types of grains. Grain collections are repeatedly suspended before settling by a pulsating fluid motion ( [START_REF] King | Modeling and simulation of mineral processing systems[END_REF]). The sorting efficiency of jigs has been studied for coarse aggregates but the identification of influencing parameters is still challenging ( [START_REF] Sampaio | Improve the quality of recycled aggregate concrete by sorting in air jig[END_REF]). [START_REF] Lyman | Review of jigging principles and control[END_REF] reviewed the jig theory and discussed the historical and more recent control systems during jigging operations. Laboratory observations are restricted to the external parts of the flow or required state-ofthe-art imaging techniques to provide information about the internal characteristics of the flow such as recirculation or convection cells. Existing models based either on the equilibrium of forces acting on grains, either on thermodynamic considerations ( [START_REF] Mayer | Fundamentals of a potential theory of the jigging process[END_REF]) require to be improved in order to better understand the stratification mechanisms ( [START_REF] Ambrós | Jigging: A review of fundamentals and future directions[END_REF]). Numerical simulations are able to complete laboratory experiments and provide information otherwise unavailable. This work is carried out in the framework of the MigFlow1 project and is devoted to the numerical simulation of water jigging with the aim to prove the efficiency of this free open-source software to compute immersed granular flows. It uses an unresolved model combining Computational Fluid Dynamic (CFD) and Discrete Element Method (DEM). Grains are considered in a Lagrangian way and their collisions are computed using the Nonsmooth Contact Dynamics ( [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF]). The fluid is computed at a larger scale than grain diameter. The fluid point properties are averaged on a representative volume following the methodology developed by [START_REF] Anderson | A fluid mechanical description of fluidized beds: Stability of the uniform state of fluidization[END_REF]. The volume fraction of the solid phase in the representative volume is used to write the fluid equations in a volume-averaged form. Such an unresolved model is able to give insight in the microscale solid dynamics due to the discrete representation of the grains while the representation of the fluid at a coarse scale is computationally convenient. The use of continuum rather than a point equations for the fluid flow requires the introduction of an explicit volume-averaged fluid-grain interaction force ( [START_REF] Xu | Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics[END_REF]) based on empirical formulas ( [START_REF] Van Wachem | Methods for multiphase computational fluid dynamics[END_REF]). Unresolved CFD-DEM models are convenient in applied sciences and have been widely used to describe industrial devices like fluidized beds ( [START_REF] Xu | Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics[END_REF][START_REF] Zhu | Discrete particle simulation of particulate system: theoretical developments[END_REF][START_REF] Kafui | Discrete particle-continuum fluid modelling of gas-solid fluidised beds[END_REF][START_REF] Li | Gas-particle interactions in dense gas-fluidized beds[END_REF]) or dense particulate reaction system ( [START_REF] Zhong | CFD simulation of dense particulate reaction system: Approaches, recent advances and applications[END_REF]). The experiments of density sorting achieved using the water jig are described in Section 2. The apparatus is presented and the sorting rate is discussed for two types of binary granular mixtures with two different initial configurations of the pile. In section 3, the numerical model and its Finite Element formulation are developed. A special attention is paid to the treatment of the boundary conditions. Finally, the simulation results are presented in Section 4. The numerical parameters of influence relevant to the computation of the laboratory experiments are varied to observe the sensitivity of the process.

Experimental observations

The experimental setup (Figure 1) is a laboratory scale water jig (Allmineral Alljig) consisting of a tank of size 33 × 33 × 35 cm. The tank is connected to a pump able to inject and drain water through a wire grid at its bottom. Initially, a granular pile of height h g is placed in the tank without water covering the pile. The pump injects progressively water into the tank until a prescribed rate maintained during a prescribed duration before cutting off the inflow rate allowing the level of water going down. Under the action of the water flow, grains are resuspended before settling. This cyclic process induces the segregation of the grains with respect to their density: the grains with the smaller density end up on the top of the pile. Drainage effect is possible during the raising and lowering of water level. The walls of the tank are vertically divided in sections. Six drawers of height 2.5 cm and four drawers of height 5 cm are stacked to form the tank walls. These drawers can be used to scrape the grains layer by layer allowing to estimate the composition of the grains' mixture in each layer. In this work, the jig was operated at a fixed frequency and the thrust was adjusted (by manually tuning the water level in the reservoir) to yield a significant lifting of the grains. Experiments were made with model granular materials, i.e. spherical beads of different materials but similar size (d s ∼ 10mm): ceramic (ρ s = 3670 kg/m 3 ), glass (ρ s = 2500 kg/m 3 ), plastic (POM, ρ s = 1410 kg/m 3 ). In particular, two binary mixtures were studied, glass/plastic (density ratio 1.77) and ceramic/glass (density ratio 1.47). The initial, volumetric, proportion of the two species was kept equal to 50/50 for both mixtures. The effect of the initial configuration was also investigated: in some cases, the mixtures were initially mixed, while in other runs, the grains were initially inversely sorted (lightest grains below heaviest grains) as presented in Figure 2. In a typical run, after the grains were loaded into the tank, the water jig was activated and the grains subjected to a given number of injection-draining cycles, after which the motion was stopped. Then the material was carefully retrieved layer by layer (by means of the drawers described above) and separated manually to determine the composition for each layer and the given number of cycles. During each run, the movement of the grains was also recorded at a sidewall by means of a CMOS colour camera (Basler acA2440-75uc) in order to characterize qualitatively the solids flow pattern. Pictures were also used to estimate the water inflow dynamics. Careful inspection of the videos revealed that the jig had three operating phases. As it is switched on, the apparatus starts pushing fluid into the tank. A larger quantity of fluid is injected in the tank during the water injection than the quantity that is drained during the water withdrawal. In this initiation phase, the maximum height reached by the grains during the lifting part is reinforced after each period due to their gradual immersion in the fluid. Then, the process reaches its operating regime in which grains are lifted periodically at the same level. Finally, when the apparatus is switched off, the elevation of the grains is smoothly damped until the water is completely drained from the tank. The initiation and final phases correspond respectively to the filling of the tank by increasing water pulses and the emptying of the tank by decreasing water pulses such that the granular pile is totally covered only during the operating regime. An uncertainty surrounding the total number of injection-draining cycles achieved when analysing the experimental results therefore comes from such initiation and final phases during which the signal has a smaller amplitude. For this reason the evolution of the separation with the number of cycles given below must be taken as a qualitative indication, more than a quantitative measure. The first interesting experimental observations were made concerning the effect of the initial configuration of the materials. It was evident that the two configurations (random or inversely sorted) create different sorting dynamics. When the grains are inversely sorted at the beginning of the experiment, the lightest grains gather at the center of the tank while the heaviest grains slide along the boundary and the corners of the tank (Figure 2, top). After the cluster of lightest grains has pierced the layer of heaviest grains, the sorting is rapidly achieved. At the opposite, when the grains are randomly placed in the initial pile, there is no possibility to create a global cluster. The sorting is vertically uniform as the lightest grains diffuse in the pile until the top (Figure 2, bottom). The sorting process is slower in this case than in the previous one. As it was said above, the average composition at different heights was determined by recovering the content of each drawer, manually separating the grains and weighing the separated samples. If, for a given binary mixture, we refer to the masses of light and heavy particles as m L and m H (and to the respective densities as ρ L and ρ H ), considering that the mean particle diameters were similar, the number fraction of, say, the light particles, can be estimated as m L /ρ L /(m L /ρ L + m H /ρ H ). Figure 3 resumes the evolution of the spatial distribution of the light component based on this layer-averaged measure. Clearly, for the initially segregated mixtures, this layer-averaged index does not account for the three dimensional structure of the bed which was observed from direct inspection. However, it still gives important hints on the bed evolution in relation to the density difference between the materials. For the mixture composed of the plastic and glass beads (density ratio of 1.77) the data confirm that the system first evolves smoothly due to the slow rising of the light grains in the center and the descent of heavy grains in the corners. Then an abrupt inversion takes place and is visible from the difference between the curves corresponding to 10 and 13 cycles in Figure 3a. Figure 3b shows that, for the mixture composed of glass and ceramic beads, the evolution is slower, and the inversion of the profile less abrupt. This is due to the lower density difference between the particles (density ratio of 1.47), which reduces the segregation driving force. Finally, Figure 3c displays the evolution of the plastic / glass mixture starting from an initially mixed configuration. As discussed with respect to the visual inspection, the behavior of this configuration is completely different from the initially segregated one. Particles start to separate near the bottom and the top: after 8 cycles, the separation is limited to the top and bottom drawers. Then, the evolution is smooth, eventually passing through a linear concentration profile. Data also confirm that the evolution is longer than for the initially segregated mixture, because a steady state is obtained after roughly 25 cycles. Collective phenomena play therefore an important role for the spatiotemporal evolution of the segregation. Resuming, the water jigging experiments show that both the initial configuration of the grains and the density difference between the materials have a strong effect on the sorting dynamics. In the following, we will discuss how a numerical model based on an unresolved CFD-DEM strategy can be used to reproduce the experimental observations.

Numerical Model

Unresolved CFD-DEM models are multiscale models using a fluid spatial discretization coarser than the grains considered as discrete bodies. Grains trajectories and forces are computed in a Lagrangian way while the fluid computation is based on volume-averaged equations solved using an Eulerian method. These methods are able to represent small scale effects due to the grains configuration like clogging while the coarse fluid representation scale highly reduces the computational time compared to methods that are fully resolved on both phases. Simulations presented in this work are achieved using the open source software MigFlow. The Finite Element Method (FEM) is used to solve the fluid equations. Linear equal order interpolation functions are used for simplicity and efficiency. It provides smooth fluid fields that are suitable for the stability of the coupling with the grains dynamics compared to high order methods. Grain trajectories and solid interactions are computed using the Nonsmooth Contact Dynamics (NSCD) method [START_REF] Jean | The non-smooth contact dynamics method[END_REF]. It is a time-stepping method to solve solid interactions in which all the contacts happening in a given time-step are solved in order to prevent any overlapping of the rigid bodies. This method is emphasized in order to prevent the oscillations a pile of grains undergo when using a contact force proportional to overlapping. Moreover, it gives a completely incompressible mixture which is more adequate to fluid-grains flows. In the following simulations, a perfectly inelastic collision law is considered without restitution coefficient. More details about the model and its validation are presented in [START_REF] Constant | Implementation of an unresolved stabilised fem-dem model to solve immersed granular flows[END_REF][START_REF] Constant | Simulation of air invasion in immersed granular beds with an unresolved fem-dem model[END_REF].

Numerical Setup

Let us consider the density sorting of glass and plastic beads that are initially mixed described above. When the system reached its operating regime, the mean water elevation h at each repetition of the injection-draining process can be estimated based on the experiment, along with the period of this process. The water elevation is computed at each period as the difference between the highest level reached during the injection part and the lowest level reached during the draining part. Levels are measured at the center of the front wall and the water elevations are averaged over all the injection-draining cycles to compute h (Figure 4). Parameters measured from the corresponding laboratory experiment are:

h = 7.0 ± 0.1 cm T = 0.776 ± 0.032 s
The inflow velocity condition is approximated as a sine function of period T . The sine amplitude is determined based on the water height reached during half a period considering a uniform inflow velocity over the bottom area. The vertical fluid velocity at the bottom of the box is modelled as:

u b = hπ T sin 2π T t (1)
All the injection-draining cycles are achieved with full amplitude pulsations during the entire computation so that we only consider the operating regime of the device. It results that the sorting is some steps ahead in the simulation. Grains are then more influenced by the fluid in the simulations than in the experiments for an equal number of injection-draining cycles. The computational box is considered full of water and grains are completely immersed since the beginning of the process. At the top boundary, a force is imposed to the fluid to maintain the pressure at a constant level while a no-slip boundary condition is applied between the fluid and the lateral walls of the box.

Fluid Problem

Computing the fluid flow at a coarse scale gives rise to averaged fluid velocity and pressure fields. An explicit coupling is considered between the grain and fluid equations. When computing the grain trajectories, the fluid flow and the force it applies on the grains have already been computed. The grain positions are used, in the same way, to compute a smooth fluid volume fraction φ that is introduced in the volume-averaged Navier-Stokes equations [START_REF] Anderson | A fluid mechanical description of fluidized beds: Stability of the uniform state of fluidization[END_REF] describing the fluid dynamics on the computational domain Ω:

∂φ ∂t + ∇ • u = 0, in Ω, (2)
∂ρu ∂t + ∇ • ρuu φ = ∇ • (2µφd) -φ∇p • I + f + φρg, in Ω, (3)
where u is the fluid velocity multiplied by φ, p is the pressure, I is the identity tensor, ρ is the density, µ is the dynamic viscosity, f is the fluid-grain interaction force, g is the gravity and d is the symmetric gradient of u/φ. The problem solved to compute the fluid flow inside the box during the sorting is constituted by these equations associated to the boundary conditions:

pI • n = p re f I • n, on Γ t , (4)
u = v, on Γ l , ( 5
) u • n = u b , on Γ b , (6)
where Γ t , Γ l , Γ b are respectively the top, lateral, bottom walls of the box and n is the unit outward normal to Ω. As we consider a no-slip boundary condition on the lateral walls, v is the zero vector.

Finite Element Formulation

Considering a proper decomposition of the domain Ω in n e elements Ω e , the finite element approximation spaces for the velocity and the pressure are

U h = {u|u ∈ C 0 (Ω) ∪ P 1 (Ω e ) ∀e = 1, ..., n e }, P h = {p|p ∈ C 0 (Ω) ∪ P 1 (Ω e ) ∀e = 1, ..., n e },
where C 0 (Ω) is the set of continuous functions of order 0 on Ω and P 1 (Ω e ) the set of polynomial functions of order 1 on each Ω e . Defining Ûh and Ph the corresponding suitable test function spaces, the stabilized finite element formulation of the averaged Navier-Stokes equations is to find

(u h , p h ) ∈ (U h × P h ) such that ∀( ûh , ph ) ∈ ( Ûh × Ph ): < ∂φ ∂t , ph > -< u h , ∇ ph >=A + S p (7)
< ∂ρu h ∂t , ûh > + < ∇ • ρu h u h φ + φ∇ • p h I, ûh > -< 2µφd h , ∇ • ûh >=-< f h + φρ g, ûh > + B + E + G + S u (8)
where the notation < •, • > is used for the L 2 -inner product on the domain Ω while S p and S u hold for the stabilization terms needed to use equal order interpolation functions. More information about the numerical implementation of the model are described in [START_REF] Constant | Simulation of air invasion in immersed granular beds with an unresolved fem-dem model[END_REF]. The right-hand terms A, B, E, G are the boundary terms obtained from the integration by part of the gradient and divergence terms for which the notation •, • Γ is used to specify the L 2 -inner product on the boundary Γ = Γ l ∪ Γ b ∪ Γ t . Boundary conditions of the problem are then weakly satisfied through the surface integrals.

1. The impermeability of the lateral walls and the inflow rate at the bottom of the tank are enforced through the surface term of the divergence integration:

A = -u h • n, ph Γ t -u b , ph Γ b (9)
2. Integration by part is performed two times in both ways on the convection term to complete the previous condition and set the convection flux through the boundary. The first integration by part gives:

< ∇ • ρu h u h φ , ûh >= -< ρu h u h φ , ∇ ûh > + i=b,l,t ρu h φ u i • n, ûh Γ i , (10)
while the second one performed the other way around on the first term on the right-hand side considering the upwind convection velocity gives:

< ∇ • ρu h u h φ , ûh > =< ∇ • ρu h u h φ , ûh > - i=b,l,t a i ρu i φ u n i , ûh Γ i - ρu h φ u n i , ûh Γ i , (11)
where the convection velocity on each boundary

u n i = u i + u h 2 • n (12)
and the a i coefficient expresses a stability condition depending on the convection velocity characteristics:

a i = 0 if u n i < 0, 1 otherwise. ( 13
)
The normal velocity u l • n is zero on the lateral walls while the normal velocity u b • n is equal to u b at the bottom of the tank. Including the above development in the momentum equation ( 8), the boundary term associated to the convection term is:

B = i=b,l,t a i ρ φ u n i (u i -u h ) Γ i . (14)
3. A similar integration by part in both ways is achieved on the pressure gradient term to enforce the pressure condition on the top boundary:

E = -(p h -p re f )I • n, c ûh Γ t (15)
4. Finally, the surface integrals of the viscous terms are used to enforce the no-slip boundary conditions on the lateral walls:

G = -2µφd h • n, ûh Γ + Cµ h u h -v, ûh Γ l + u h -v • nn, ûh • n n Γ l (16)
where the second term is a jump on the velocity required to correct the normal gradient of the velocity computed with the inner value of u h and the third term is the similar correction applied to the transpose of the normal gradient of the velocity. The penalization coefficient C can be increased to enforce a strict satisfaction of the boundary condition. Its value will be discussed later.

Solid Problem

Solid interactions are modeled using the non-smooth contact dynamics (NSCD). The set of equations stating the dynamics of each grains is derived from Newton's second law of motion:

m s du s dt = m s g + F s + j∈{S}\s R c s j (17)
where m s is the mass and u s is the velocity of the grain s in the set of grains S, F s is the force applied by the fluid on the grain and R c s j is the contact reaction force exerted by the grain j on the grain s. Reactions between grains and boundaries are not stated in the above expressions but are treated in an analogous way to grain-grain reaction forces. The NSCD method consists in finding, at each time step, the reaction forces so that when the grains are moved, there is no interpenetration. It is based on an implicit scheme allowing a larger time step, reducing the number of solid sub-time steps achieved at each fluid step. The MigFlow software iteratively solves the contacts at each time step using a perfectly inelastic contact law. The algorithm is based on a queue listing the contacts that still have to be solved [START_REF] Constant | Implementation of an unresolved stabilised fem-dem model to solve immersed granular flows[END_REF] to improve computational efficiency compared to the usual Gauss-Seidel like algorithm [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF][START_REF] Jean | Non-smooth contact dynamics approach of cohesive materials[END_REF]. In Equation [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF], the discrete fluid-grain interaction force F s is used. Applied at the center of the grains s, the force F s exerted by the fluid on the grains s is:

F s = -V s ∇p| x s -γ(φ, x s , u s , u) u s - u φ x s ( 18
)
where V s is the volume of the grain, x s is the position of the grain, u s is the velocity of the grain and γ(φ, x s , u s , u) is the drag coefficient. This force is related to the continuous fluid-grain interaction force f in Equation ( 3) by the expression:

f = - s∈S F s + V s ∇p| x s δ s ( 19
)
where δ s is the Dirac function at the position of the grain s.

Computational result

Let us compute the density sorting of glass (ρ s = 2500 kg/m 3 ) and ceramic (ρ s = 3670 kg/m 3 ) beads. A pile of glass beads of 14.392 kg is initially placed under a pile of the ceramic beads of 19.300 kg so that the sorting is achieved when the configuration is inverted. In the simulations, granular piles consist of about 22000 spherical grains having a diameter uniformly distributed in the range:

d s ∈ [9, 11] • 10 -3 m
to represent the slight polydispersity due to the manufacturing process. In order to respect the total mass of each material, the number of grains varies from one deposit to another because of the random distribution of the grains diameter. The particle Reynolds number in the simulations is:

Re s = ρ u b,max d s µ = 1000 • 0.283 • 11 • 10 -3 10 -3 = 3113. ( 20
)
Figure 5 shows the configuration of the grains at different stages of the process. The computation provide the same dynamics than the one observed in the experiments. The least dense grains (the glass beads) progressively gather at the center of the box while the densest grains (the ceramic beads) slide along the corners of the box. As the injectiondraining cycles continue, the least dense grains pierce the layer of densest grains through its center (Figures 7). At the end of the process, the granular pile reaches an inverted configuration compared to the initial one.

Comparisons show that simulations and laboratory experiments are qualitatively comparable. Achieving a similar analysis of the layer-averaged number fraction of light particles as a function of the height for the simulations than for the laboratory experiments (Figure 3) gives a quantitative comparison. Figures 6 (a) and (b) compare the layeraveraged number fraction of light grains with respect to the height for a glass-plastic binary pile initially inversely segregated and mixed respectively for both simulations and experiments. It can be seen that similar vertical profiles are found but for different numbers of periods. The main difference between simulations and laboratory experiments is that grains are completely immersed since the beginning in simulations which is not the case in the laboratory experiments. Moreover, as explained above, due to the progressive water injection, there exists some uncertainty on the number of periods in the laboratory experiments. It can explain the difference between the number of periods required to obtain the same vertical profiles in the simulations and laboratory experiments.

Sensitivity analysis

Some characteristics influencing the stratification such as particle density, size or shape as well as bed and tank dimensions are listed by [START_REF] Ambrós | Jigging: A review of fundamentals and future directions[END_REF]. This paper does not pretend to give insight in all the influencing parameters of the jigging process but only the numerical and physical parameters relevant to the full description of the specific laboratory experiments described above.

In the following, the simulations of the density sorting allow to estimate the requirements to obtain the sorting dynamics observed in the laboratory experiments and the optimal parameters speeding the process up. First, the effect of the drag force is described. It is important because jigging relies on the different suspension and settling rates to separate grains with respect to their density. Jigging creates successive expansions and contractions of the pile so that the sorting is mainly achieved when the bed is fluidized. We will estimate if the friction between rigid bodies is still able to disturb the dynamics. The density ratio of the two types of grains and the initial configuration of the pile are then considered. Laboratory experiments achieved with initially mixed binary piles show a very different sorting dynamics than the ones achieved with initially inversely segregated binary piles. Numerical simulations are used to give insight in the sorting rate for both cases. Finally, the boundary treatment is considered through the variation of the penalization coefficient. It will be shown that the sorting dynamics is really dependant on the no-slip boundary condition enforced on the lateral walls.

Drag coefficient

The fluid-grain interaction force appears in the momentum Equation [START_REF] Yamada | Pinched flow fractionation: continuous size separation of particles utilizing a laminar flow profile in a pinched microchannel[END_REF]. Neglecting the Basset force, virtual mass effect and the particle lift forces [START_REF] Crowe | Multiphase flows with droplets and particles[END_REF], the fluid-grains interaction force is constituted by the pressure gradient force and the drag force. There is no agreement in the scientific community about an expression of the drag force applied by a fluid on a set of grains for all flow regimes and the formula we use is empirically determined. The coefficient γ in Equation ( 18) depends on the grain Reynolds number and parametrizes the drag force with respect to the flow regime [START_REF] Di Felice | Fluid-particle drag force in binary-solid suspensions[END_REF]. [START_REF] Richardson | The sedimentation of a suspension of uniform spheres under conditions of viscous flow[END_REF] make the general hypothesis that the drag force for a granular assembly is computed as the drag force acting on a single grain multiplied by a function of the porosity. A simple power law of the porosity multiplies the drag coefficient given by [START_REF] Dallavalle | Micromeritics: the technology of the particles[END_REF] representing a smooth transition between the Stokes and Newton regimes so that:

γ(φ, x s , u s , u, a) = φ -β x s A s ρ 2         0.63 u s - u φ x s 1 2 + 4.8 µ d s ρ φ| x s 1 2         2 (21)
where the value of β = 1.8 is given by [START_REF] Wen | Chemical engineering progress symposium series[END_REF]. Due to the empirical nature of the drag force, its value can be easily questioned as it potentially biased the numerical results. To determine the influence of the drag on the numerical results, the γ parametrization is multiplied by a parameter C d . Let us consider the case of glass (ρ s = 2500 kg/m 3 ) and ceramic (ρ s = 3670 kg/m 3 ) beads that are initially inversely segregated. The density sorting is achieved due to the fluid-grain interaction and particularly due to the different drag forces exerted on the two types of grains. Figure 8 shows the evolution of the mean density ρsup in the top half of the box. It can be seen that decreasing the drag force applied on the grains smooths the sorting process so that it takes more time to be completed. Not only the beginning of the process is delayed but the sorting rate (slope) is also changed. Figures 9 (a) and (b) show the total kinetic energy of the light and heavy grains respectively during the simulation. It can be seen that the solid kinetic energy decreases with the drag coefficient while the fluid energy remains unchanged.

Friction

Solid grains interact through collisions constituted by an inelastic restitution in the normal direction and a friction in the tangential direction that depends on the roughness of the grain surface. Using manufactured spherical grains and in presence of an interstitial fluid, the friction coefficient of the materials is difficult to estimate. However, it is possible to investigate numerically the effect of the friction on the sorting dynamics. Considering binary mixtures of grains, there are five friction coefficients susceptible to change the sorting dynamics. Each grain type can interact with grains of the same type, grains of the other type and walls. We still consider the sorting of ceramic and glass beads that are initially inversely segregated so that the five friction coefficients to investigate are the friction coefficients between glass-glass beads = µ gg , glass beads-plexiglas walls = µ gw , glass-ceramic beads = µ gc , ceramic-ceramic beads = µ cc , ceramic beads-plexiglas walls = µ cw .

As a first investigation, let us consider that all the friction coefficients are the same and let us note it µ s . Figure 10 shows the sensitivity of the process to the friction coefficient µ s . It can be observed that the sorting process is slower for small values of µ s and we conclude that friction effects play a crucial role in the process dynamics. It is interesting now to determine if this sensitivity is due to one particular coefficient or if it results from a combined effect of several friction coefficients. To investigate the individual effect of each friction coefficient, different values of each coefficient are used keeping the values of all the other coefficients at a value of 0.5. For example, treating the glass-glass friction coefficient, we test the values µ gg = 0.1, 0.3, 0.5, 0.7, 0.9 with µ gw = µ gc = µ cc = µ cw = 0.5. In this way, we obtain the sensitivity of the process to each coefficient individually. The presence of a fluid creates lubrication between grains so that friction is highly decreased ( [START_REF] Brenner | The slow motion of a sphere through a viscous fluid towards a plane surface[END_REF]). That is why the friction coefficient is only varied in the range between 0.1 and 0.9. Figure 11 shows the sensitivity of the time evolution of the density sorting to the variation of the friction coefficients. The shaded area presents the variation of the time evolution curves when varying µ gg , µ gw , µ gc and µ cc . Not only the sensitivity to these parameters is small but there is no clear tendency in the variation which can be explained due to the initial configuration. Due to the initial configuration of the bed, grains of each type remain agglomerated during all the sorting process. It results that glass-ceramic contacts are limited to the interface between the agglomerates. Moreover, the motion of glass beads relatively to ceramic beads mainly take place during the suspension and deposition phases i.e. during the fluidization of the pile so that contacts hardly occur. That is why µ gc does not affect significatively the sorting process. Analysing the sorting dynamics for initially inversely segregated piles, we might think that increasing the friction between grains of the same specie will result in a slower sorting process. Indeed, during the sorting we observe that glass beads gather at the center of the tank before piercing the layer of the densest grains. As a result, the densest grains are compacted near the corners of the box. Increasing the friction coefficient between the grains of the same specie might then results in a more difficult gathering of the grains. However, this slow-down is not so clear and varying the friction between grains of the same specie provide quite unpredictable results. The least dense grains moving at the center of the tank, it was expected that their friction with the walls is not relevant so that µ gw does not impact the results. However, the densest grains are compacted against the walls during the sorting process. By increasing their friction with the walls, motions of the densest grains near the wall are prevented. As a result, the inflow rate being the same, the motion at the center of the tank is increased making the sorting faster. The coloured lines on Figure 11 show that varying µ cw keeping the other friction coefficients constant has a huge impact on the sorting duration (much more important than the variation induced by the other coefficient represented by the shaded area). Regarding the kinetic energy of the light and heavy grains during the beginning of the simulation (Figures 12 (a) and (b)), it is clear that increasing µ cw tends to favour motion of the light grains. Their kinetic energy increases with µ cw while the kinetic energy of heavy grains is decreased. Varying µ cw has then an impact on the relative motion of the two species during the sorting process. Not only the sorting rate is affected by the friction between heavy grains and the walls but also the sorting dynamics. This can be reflected by the grain velocity distribution of each specie. Let us compute a continuous representation of the time-averaged grain velocity projected in the horizontal plane. This continuous solid velocity at each point {i, j} of a regular grid in the horizontal plane is obtained as

V i, j = ∆t T end t s 1 πh 2 exp - (X i -x s ) 2 + (Z j -z s ) 2 h 2 u s πd 3 s 6 ( 22
)
where X i and Z j are the coordinates of the grid point, x s and z s are the coordinates of the grains, T end is the total simulation time and h is the smoothing length of the discrete field. Figures 13 (a)-(e) show the continuous representation of the velocity for different µ cw . These figures reveal that for a value of µ cw much lower than the other friction coefficients, the motion of the light grains is not favoured at the center of the box anymore. Light grains are able to cross the bed near the corners of the box and the velocity field is more spread, showing a Poisson like distribution. It is another proof that increasing µ cw with respect to the other friction coefficients tends to favour the motion of light grains at the center of the box, increasing the sorting rate. Finally, the contact network is compared for each value of µ cw . Figures 5 (a)-(e) show the contact network at the early stage of the negative part of the inflow signal at 3.6 s. It can be seen that the contact network is reinforced near the walls by the increase of µ cw but also that the network shapes are different. Contact network shows that the top of the bed is more rounded while increasing µ cw due to light grains crossing the bed at the centre of the box.

Initial pile properties

The sorting dynamics is very influenced by the initial state of the pile. Presenting the experimental results, we briefly compared the sorting of an inversely segregated deposit with the sorting of a mixed deposit. It has been observed that the sorting is slower for the case of an initial mixed pile of grains. The density ratio of the pile constituents (defined as the ratio of the larger density over the smaller one) has also been shown to influence the sorting rate.

Figure 19 shows the time evolution of the mean density in the upper half of the pile for the sorting of ceramic-glass and plastic-glass piles. Beads diameter in the plastic-glass piles are similar to beads diameter in ceramic-glass piles and consist of 14.392kg of glass beads and 8.117kg of plastic beads.

As it has been assumed in the laboratory observations, the sorting is smoother for the initially mixed pile. As the grains of the same type are scattered in the pile, there is no possibility to form large clusters. The grains have to diffuse alone in the pile till the top which gives a slower and more continuous sorting dynamics. Figures 5 (a) and (b) present the trajectories of the grains projected in the z=x plane (y is the vertical coordinate) for both initial configurations. In the case of an initially inverted configuration, two convection cells are observed near the boundaries constituting a symmetric pattern in the box while trajectories are more chaotic in the case of an initially mixed configuration. The total kinetic energy of the grains during the simulation for light and heavy grains in both configurations are also observed in Figures 16 (a) and (b) respectively. For an initially inverted configuration, an increase of the peaks in the solid kinetic energy is observed when the sorting is suddenly achieved due to the collective phenomenon. At the opposite, the solid kinetic energy for an initially mixed binary granular mixture oscillates between almost uniform values which confirmed a smooth sorting.

The density ratio also appear to change the sorting rate. The larger the density ratio, the faster the sorting. This observation is in agreement with what it is expected. Jigging is a sorting process based on gravity concentration [START_REF] Sampaio | Improve the quality of recycled aggregate concrete by sorting in air jig[END_REF]. Using water jigging, the sorting depends on the gravity and hydrodynamic forces during the expansions and contractions of the pile. A larger difference of density between the different types of grain result in a larger difference of the forces applied on the grains. The separation of the grains is then faster. The continuous representation of the solid velocity also confirmed the different patterns for the two initial configurations. The continuous velocity fields of the heavy and light grains for both initial configurations are shown in Figures 18 and17. It can be seen that for an inverted configuration, the maximum velocity of the light grains is located near the middle of the box with rounded levels while the maximum velocity of heavy grains are located near the corners. The velocity reflects then exactly the paths chosen by the two types of grains. For the mixed configuration, the velocity exhibits a Poisson profile expected for the velocity of a continuun medium in a square section with partial slip walls. It tends to favour the hypothesis of a smooth process without dominant flow area and collective phenomena.

Partial slip coefficient at the lateral boundaries

The boundary condition ( 6) is a Dirichlet condition corresponding to a no-slip wall boundary condition on the lateral wall of the box. Usually, such a condition is strictly enforced on the boundary. According to the finite element formulation developed above, it requires to set a large penalization coefficient for the jump terms. Computing the density sorting with a penalization coefficient C = 10 6 gives results presented in Figure 20. It can be seen that the evolution dynamics at short times corresponds to the laboratory results. The least dense grains go through the pile of densest grains at the centre of the box to achieve the sorting. However, it can be seen that if the injection-draining cycles repeat several times after the sorting is complete, the least dense grains slide near the boundaries of the box while the centre part of the pack stay sorted.

In order to comply with the requirement of the model, the fluid discretization scale has to be larger than the grain diameter. As a result, it is not possible to have a fully resolved boundary layer of the fluid flow near the lateral walls. By strictly enforcing the Dirichlet boundary condition, the prescribed value affects the flow much farther from the wall than it should be. The velocity is more important at the centre of the box. When grains of the centre of the pile are lifted by the fluid, grains at the top of the pile slip downwards in the near wall area. In order to avoid the slipping of the grains along the lateral walls, it is required to reduce the region of influence of the Dirichlet boundary condition. Due to the poor resolution of the boundary layer, it is more appropriate to allow a small deviation between the computed solution on the boundary and the Dirichlet boundary condition by specifying a smaller penalization coefficient.

Figure 21 shows the computation of the density sorting obtained with a coefficient C = 1. It can be seen that the sorting dynamics is inverted compared to the laboratory results. Due to the small penalization coefficient, the no-slip condition is not sufficiently enforced. The velocity near the boundaries is too high and especially near the corners where the volume fraction of fluid is larger. As a result, the least dense grains are carried by the fluid in the areas where the flow is the most important.

Varying the penalization coefficient gives the opportunity to estimate the transition point between the two dynamics.

Let us consider the division of the domain in two sub-volumes Ω ext and Ω int (Figure 22) so that we define

ρ int (t) = s∈Ω int m s s∈Ω int V s (23)
ρ ext (t) = s∈Ω ext m s s∈Ω ext V s (24) (25)
that are the mean grain densities in Ω int and Ω ext with m j and V j the mass and the volume of the grain j respectively. These mean grain densities vary in time due to the motion of the grains. If the least dense grains gather at the centre of the tank before piercing the layer of densest grains like in the case of a strict satisfaction of the Dirichlet boundary condition, the time-integrated ratio of the mean densities computed between the start time t 0 and end time t end of the simulation:

r = 1 t end -t 0 t end t 0 ρ int ρ ext dt (26)
is greater than one and lower than one otherwise. Figure 23 shows the evolution of the average value of r using five different deposits that are initially inversely segregated. It can be seen that the transition point is around a value of C = 350.

Conclusion

The density sorting by water jigging of binary mixtures constituted of manufactured spherical grains has been studied both experimentally and numerically. Laboratory experiments are used to study the density sorting of granular piles using two different initial configurations. The separation of binary piles initially mixed is smooth while the collective phenomena play an important role for pile initially inversely sorted. Simulations are used to determine the sensitivity of the process to key parameters relevant for the chosen experimental setup. The results show an important dependence of the sorting rate on the initial configuration of the pile. The density sorting is faster when the least dense grains are able to gather and form agglomerates rising to the top of the pile. At the opposite, dispersing grains in the initial pile, they have to diffuse upwards alone in the pile. The pile being fluidized by the water pulses, this can be understand as the rising of an air bubble in fluid. The more grains there are in clusters, the larger the cluster diameter is and the larger the rising velocity is. The sorting rate benefits then from the formation of aggregates.

Water jigging relies on the difference in suspension and settling velocities of grains with different densities to separate granular mixtures. As expected, the sorting rate has been showed dependent on the drag force intensity and the density ratio between the two types of grains. Increasing the density ratio or the drag force intensity result in a faster density sorting.

More surprising was the results obtained using different friction coefficients. It has been shown that, considering an initially inversely segregated ceramic-glass pile, the friction coefficients implying the least dense (glass) grains have an insignificant influence on the sorting rate as the friction coefficient between ceramic and ceramic. However, the friction coefficient between ceramic and the walls has been shown influencing the sorting rate. It can be assumed that increasing the friction between the ceramic beads and the walls prevents motion of the pile near the walls. The water flow is then favoured at the centre of the tank with a larger velocity. It results a larger drag force at the centre of the tank because it is proportional to the fluid-grain relative velocity so that the sorting is faster. This hypothesis has been confirmed by the analysis of the kinetic energy, solid velocity field and contact network during the sorting process. It reveals that increasing the ceramic-wall friction coefficient with respect to the other friction coefficients tends to favour motion of the light grains at the center of the box while decreasing it is able to change the sorting dynamics by allowing motion of the light grains near the center of the box. Finally, the numerical treatment of the boundary conditions has been studied. Boundary conditions enforcement may seem trivial but due to the coarse resolution of the fluid flow, small scale phenomena have to be modeled explicitly.

In particular, the boundary layer is underresolved. It results that strongly enforcing the no-slip boundary condition create a too large boundary layer while weakly enforcing the no-slip boundary condition create a too small boundary layer. This gives rise to completely different sorting dynamics. It confirms that boundary effects due to solid or fluid friction completely trigger the sorting pattern. Further researches should probably focus on the value of the penalization coefficient used to enforce the no-slip boundary condition and compare it to the wall shear stress. The sorting dynamics might then be related to the ratio between the drag force at the centre of the tank and the wall shear stress. Comparisons with models able to resolve the flow near the boundaries at a fine scale should also be considered to estimate the shear stress between the wall and the fluid. At this end, semi-resolved models based on a volumeaveraging of the discrete properties at a scale which is independent from the fluid discretization [START_REF] Wang | A semi-resolved cfd-dem approach for particulate flows with kernel based approximation and hilbert curve based searching strategy[END_REF][START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF][START_REF] Esteghamatian | A stochastic formulation for the drag force based on multiscale numerical simulation of fluidized beds[END_REF] seem to be very promising approaches to consider in the future. (e) µ cw = 0.9

Figure 1 :

 1 Figure 1: The experimental apparatus: a laboratory scale water jig (Allmineral Alljig).

Figure 2 :

 2 Figure 2: Evolution of the density sorting of glass and plastic beads that were (top) initially inversely sorted and (bottom) mixed. A cycle is constituted of a resuspension phase and a settling phase.

Figure 3 :Figure 4 :

 34 Figure 3: Evolution of the layer-averaged number fraction of light grains for the different mixtures studied in this work: (a) plastic / glass, initially inversely segregated, (b) glass / ceramic, initially inversely segregated, (c) plastic / glass, initially mixed. The legends provide the number of cycles performed.

Figure 5 :

 5 Figure 5: Density sorting of glass (cyan) and ceramic (grey) beads in water at different stages of the computation achieved with a penalization coefficient C = 600 and without friction.

Figure 6 :

 6 Figure6: Layer-averaged number fraction of light grains during the sorting process. Solid lines correspond to simulations while dashed lines correspond to laboratory experiments. The Pearson correlation r x 1 ,x 2 is given for comparable curve to asses the agreement between the vertical profiles.

Figure 7 :Figure 8 :

 78 Figure 7: Cuts of the pile to show the grain configuration after 8 injection-draining cycles

Figure 9 : 9 Figure 10 : 9 Figure 11 :

 9910911 Figure 9: Total kinetic energy of the (a) light and (b) heavy grains for different drag coefficients during the simulation

Figure 12 :

 12 Figure 12: Kinetic energy of (a) the glass (light) and (b) the ceramic (heavy) grains during the first 10 s of the simulation using different values of µ cw using an inverted initial configuration

Figure 13 :

 13 Figure 13: Time-averaged continuous representation of the velocity field [m/s] of glass (light) beads during the sorting process of a binary pile initially inversely sorted for different friction coefficients between ceramic (heavy) beads and the walls

Figure 14 :Figure 15 :

 1415 Figure 14: Contact network in the ceramic-glass mixture at the early stage of the negative part of the inflow signal at 3.6 s for different values of the ceramic-wall friction coefficients.

Figure 16 :

 16 Figure 16: Total kinetic energy of the (a) light and (b) heavy grains during the simulation for an initially inversely sorted and an initially mixed configurations

Figure 17 :

 17 Figure 17: Time-averaged continuous velocity of the solid phase during the sorting simulation of a binary granular plastic-glass bed initially mixed

Figure 18 :

 18 Figure 18: Time-averaged continuous velocity of the solid phase during the sorting simulation of a binary granular bed initially inversely sorted

Figure 19 :

 19 Figure19: Time evolution of the mean density in the upper half of the pile for initial piles constituted of ceramic-glass beads and plastic-glass beads. Orange lines correspond to a pile with an initial mixed configuration while the blue lines correspond to an initially inversely sorted pile. The different curves for each case hold for different initial piles of a given configuration.

Figure 20 :

 20 Figure 20: Density sorting of glass (cyan) and ceramic (grey) beads at different stages of the computation using C = 10 6 .

Figure 21 :Figure 22 :Figure 23 :

 212223 Figure 21: Density sorting of glass (cyan) and ceramic (grey) beads at different stages of the computation using C = 1.

www.migflow.be

Acknowledgements

Matthieu Constant is a Research Fellow with the Belgium Fund for Research in Industry and Agriculture (FRIA).