V Loisel

M Abbas

) O Masbernat

E Climent

Inertia-driven particle migration and mixing in a wall-bounded laminar suspension flow

Laminar pressure-driven suspension flows are studied in the situation of neutrally buoyant particles at finite Reynolds number. The numerical method is validated for homogeneous particle distribution (no lateral migration across the channel): the increase of particle slip velocities and particle stress with inertia and concentration is in agreement with former works in the literature. In the case of a two-phase channel flow with freely moving particles, migration towards the channel walls due to the Segré-Silberberg effect is observed, leading to the development of a non-uniform concentration profile in the wall-normal direction (the concentration peaks in the wall region and tends towards zero in the channel core). The particle accumulation in the region of highest shear favors the shear-induced particle interactions and agitation, the profile of which appears to be correlated to the concentration profile. A 1D model predicting particle agitation, based on the kinetic theory of granular flows in the quenched state regime when Stokes number St = O(1) and from numerical simulations when St < 1, fails to reproduce the agitation profile in the wall normal direction. Instead, the existence of secondary flows is clearly evidenced by long time simulations. These are composed of a succession of contra-rotating structures, correlated with the development of concentration waves in the transverse direction. The mechanism proposed to explain the onset of this transverse instability is based on the development of a lift force induced by spanwise gradient of the axial velocity fluctuations. The establishment of the concentration profile in the wall-normal direction therefore results from the combination of the mean flow Segré-Silberberg induced migration, which tends to stratify the suspension and secondary flows which tend to mix the particles over the channel cross section.

I. INTRODUCTION

The behavior of finite size non-Brownian particles suspended in a Newtonian fluid raises many questions from non-Newtonian response [START_REF] Barnes | Shear-thickening (dilatancy) in suspensions of nonaggregating solid particles dispersed in Newtonian liquids[END_REF] to laminar-turbulent transition modulation. [START_REF] Matas | Transition to turbulence in particulate pipe flow[END_REF] Such effects are of prime importance in two-phase transport applications, separation, oil recovery, or fracking processes. At vanishing particle Reynolds numbers, channel flows laden with neutrally buoyant particles remain homogeneous, unless the concentration is high enough to cause shearinduced migration (whether they are Brownian or not), towards the channel center. [START_REF] Frank | Particle migration in pressure-driven flow of a Brownian suspension[END_REF][START_REF] Leighton | The shear-induced migration of particles in concentrated suspensions[END_REF][START_REF] Nott | Pressure-driven flow of suspensions: Simulation and theory[END_REF] At finite Reynolds numbers, the finite size particles migrate under dilute conditions in pressure-driven flows, a phenomenon that was first observed by Ref. 35 in pipe flows. The single particle inertial effect has been used in some applications like particle sorting and separation in microfluidic channel flows (see the review of Ref. 11). Despite the small scales involved in these applications (channels with typical dimensions 100 µm and particles of diameter 10 µm), typical axial velocities of a) micheline.abbas@ensiacet.fr O(10-100 cm/s) yield high shear rates such that the channel Reynolds number Re is of O(10-100) and the particle Reynolds number Re p is O(0.1-1). While the micro-hydrodynamic phenomena for inertia-driven particle migration are well documented in dilute conditions, important issues related to collective effects at moderate concentration remain unclear. This work addresses the behavior of moderately concentrated suspensions of finite size particles in pressure-driven flow.

Early experiments of Ref. 35 showing preferential equilibrium positions of particles in pipe flows have given rise to a number of studies on isolated spherical particles in both bounded and unbounded shear flows. Theoretical predictions (at low but finite particle Reynolds number) have shown that the particle radial equilibrium position in a Poiseuille flow results from a balance between inertial lift force induced by the interaction of the particle stresslet with the shear gradient flow, and a wall-induced hydrodynamic repulsion. [START_REF] Asmolov | The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number[END_REF][START_REF] Ho | Inertial migration of rigid spheres in two-dimensional unidirectional flows[END_REF][START_REF] Schonberg | Inertial migration of a sphere in Poiseuille flow[END_REF][START_REF] Vasseur | The lateral migration of spherical particles sedimenting in a stagnant bounded fluid[END_REF] The equilibrium position moves towards the tube wall when the Reynolds number and/or the particle radius are increasing. More recently, the experimental study of Ref. 25 carried out in the range Re p = O(1-10) (where Re p = γa 2 /ν, γ being the local shear rate experienced by a particle of radius a, and ν the fluid kinematic viscosity) has shown that more complex mechanisms arise at intermediate particle Reynolds numbers, where a second equilibrium position emerges closer to the tube centerline. Lift force intensity and particle equilibrium position as a function of flow Reynolds number have been numerically investigated in a pipe geometry [START_REF] Shao | Inertial migration of spherical particles in circular Poiseuille flow at moderately high Reynolds numbers[END_REF] and theoretically [START_REF] Asmolov | The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number[END_REF][START_REF] Matas | Lateral force on a rigid sphere in large-inertia laminar pipe flow[END_REF] in both pipe and channel geometries. It was clearly demonstrated that the lift force magnitude increases while the Reynolds number increases, and that it is smaller in pipe than in channel geometry. At a given Reynolds number, it results in an equilibrium position closer to the wall of a plane channel than of a cylindrical pipe. The existence of a second equilibrium position in pipe geometry was shown to be due to finite size effect of particles.

In dilute conditions (0.1% concentration), particle pair interactions at finite inertia lead to the structuring of particles in chain-like patterns that develop in the streamwise direction. [START_REF] Matas | Trains of particles in finite-Reynolds number pipe flow[END_REF][START_REF] Shao | Inertial migration of spherical particles in circular Poiseuille flow at moderately high Reynolds numbers[END_REF] The characteristic separation length within this particle alignment, which occurs around the radial equilibrium position, probably results from a combination of the narrow band of equilibrium position and the development of reverse streamlines upstream and downstream of a particle at finite particle Reynolds number. The separation length between these particles decreases as the particle Reynolds number is increased, following the evolution of the distance between the nearest reverse streamline and the leading particle. [START_REF] Matas | Trains of particles in finite-Reynolds number pipe flow[END_REF] The rheological properties of a dilute unbounded sheared suspension at finite inertia have been investigated, based upon theoretical calculations. [START_REF] Lin | Simple shear flow round a rigid sphere: Inertial effects and suspension rheology[END_REF][START_REF] Subramanian | The influence of the inertially dominated outer region on the rheology of a dilute dispersion of low-Reynolds-number drops or rigid particles[END_REF] Expressions of shear viscosity and normal stress differences have been obtained to order O(φRe 3/2 p ), by dividing the velocity perturbation field into an inner region around the particle where Stokes equations provide the leading order terms, and an outer region where inertial effects are dominant. A shear-thickening behavior of the suspension is observed, the shear viscosity being a growing function of the shear rate (or particle Reynolds number) even in the infinitely dilute regime. Normal stress differences are finite, the first and second differences being, respectively, negative and positive. In the limit of infinitely dilute suspension, the effect of inertia therefore induces a non-Newtonian response.

The rheology at finite concentration is more complex to describe since it depends on the suspension microstructure. The signature of the microstructural anisotropy is a deficit of particles in the extensional quadrant, i.e., downstream side of a reference particle. [START_REF] Brady | Microstructure of strongly sheared suspensions and its impact on rheology and diffusion[END_REF] Therefore, normal stresses develop, leading to particle migration towards the low shear regions in curvilinear sheared flows. It is shown that Stokesian suspensions develop heterogeneity under the influence of hydrodynamic and non-hydrodynamic (like contact forces) interactions, with or without Brownian motion (this subject being addressed among many other works by Refs. 29, 30, 36, and 41). Nevertheless, the microstructural anisotropy and thus the normal stress differences increase with the particle Reynolds number (as shown by the numerical simulations of Refs. 17 and 42). Consequently, the particle shear-induced migration is expected to compete with the inertia-driven migration (Segré-Silberberg effect) in a complex way.

The study of spatially developing suspension flow in pipe or channel geometry at finite Reynolds number and moderate concentration has been poorly addressed. Combination of inertial lift forces and finite concentration leads to collective motion and heterogeneity in the flow cross section, which have been evidenced a few times in the literature. At finite concentration, Ref. 15 has shown that in a laminar suspension flow in a pipe, established concentration profiles are strongly dependent upon the particle Reynolds number, for Re p > 0.1, in both dilute or concentrated regimes. Under these conditions, low Reynolds number shear-induced migration or suspension balance models fail to reproduce their experimental data.

In channel flow geometry and at finite concentration, the wall-normal inertial migration will induce a concentration gradient across the channel, since particles tend to migrate towards a single equilibrium position close to the channel walls where the shear rate is the highest. It is then expected that shear-induced particle interactions will lead to particle velocity fluctuations, and consequently diffusion of particles towards the channel core (both the shear-induced velocity fluctuations and diffusion being weakly dependent on the particle Reynolds number as shown by Ref. 42). At steady state and for a statistically steady flow in the channel, a momentum balance in the wall normal direction corresponds to a balance between the inertial lift force and the shear-induced stress flux, resulting in a concentration profile, which is so far unexplored. In this work, we use numerical simulations based on the Force Coupling Method (FCM) to study the development of a suspension flow, at finite Reynolds number and moderate concentration. We show that the suspension exhibits unexpected flow features, which clearly put into question the relevance of a one-dimensional momentum balance approach.

We start the paper in Section II by briefly explaining the force coupling method used for the two-phase flow simulations and the tests used to validate the case of a single particle in shear flow at finite particle Reynolds number. In Section III, we present a set of short time simulations at finite concentration using Couette and Poiseuille configurations (wall-bounded flows). At short time scales, the suspension homogeneity is relatively preserved due to small values of particle-fluid slip velocities. Concentration gradients are therefore, negligible and the particle stress and slip velocity profiles can be obtained in quasi-homogeneous conditions. By comparing the flow response in both Couette and Poiseuille flows, we were able to analyze the effect of the linear or quadratic velocity profile on the suspension flow properties. In Section IV, the pressure-driven suspension flow is simulated for a long time and its transient evolution is carefully described and discussed. After particle migration towards the channel wall is completed, the suspension flow becomes stratified into clear and particle-laden zones, resulting in concentration profile peaking near the wall. Assuming the concentration profile is established and steady, a 1D model of particle fluctuating kinetic energy in the wall-normal direction is derived and compared to the numerical predictions. It is shown that such a model cannot reproduce the simulated concentration and agitation profiles in the wall normal direction. The development of secondary flows at longer simulation times is shown to be responsible of remixing the particles in the channel cross section. The origin of these secondary flows is discussed and an interpretation is proposed based on inertia-driven mechanism.

All the simulations were performed with one particle-to-channel size ratio H/a = 32, where H is the channel height, and a moderate concentration φ b = [1%-5%]. Along this paper, φ b refers to the average volumetric solid fraction, in the computation domain, of a suspension simulated numerically, and φ denotes local solid volume fraction or concentration in general. The channel Reynolds number Re = U H/ν is in the range [1000-1400], which is large enough to promote the evolution of the suspension structure, and just below the laminar-turbulent transition ( U stands for the flow average velocity) which has been studied in a previous work recently published by the same authors. [START_REF] Loisel | The effect of neutrally-buoyant finite-size particles on channel flows in the laminar-turbulent transition regime[END_REF] The particle size and concentration were chosen such that the particles significantly reduce the transition threshold of the suspension flow (see Ref. 24).

II. SIMULATION METHOD AND VALIDATION FOR A SINGLE PARTICLE

The numerical scheme has been presented in Ref. 20. Its main features are recalled in this work for the sake of completeness. The numerical simulation of particle trajectories and suspension flow field is based on the FCM (which is a multipole expansion of momentum source terms added to the Navier-Stokes equations). Flow equations are dynamically coupled to Lagrangian tracking of particles (see details in Refs. 27, 21, and 10). The fluid is assumed to fill the entire simulation domain, including the volume occupied by the particles. The fluid velocity and pressure fields are solutions of the mass and momentum conservation equations (Eqs. ( 1) and ( 2)),

∇ • u = 0, (1) ρ Du Dt = -∇p + µ∇ 2 u + f (x,t) . (2)
ρ and µ are, respectively, the density and dynamic viscosity. u represents the velocity vector. The presence of the dispersed phase in the fluid is then represented by a body force distribution f(x,t) written as a multipole expansion (Eq. ( 3)) where the first term is the monopole F (n) representing the force that the nth particle applies on the fluid (due to an external forcing or particle-to-particle contact forces). The second term is the dipole G (n) tensor. Its anti-symmetric part is related to external torques applied on the particle. The symmetric part is set through an iterative procedure to ensure that the strain-rate within the fluid volume occupied by the dispersed phase is zero (enforcing solid body response),

f i (x,t) = N n=1 F (n) i ∆ x -Y (n) (t) + G (n) i j ∂ ∂ x j ∆ ′ x -Y (n) (t) . (3)
Y (n) is the nth particle center position. ∆(x -Y (n) ) and ∆ ′ (x -Y (n) )
, the spatial distributions of the monopole and dipole forcing, respectively, are generated using finite size envelopes that are adapted to the shape of the spherical particle. The particle translation and rotation velocities V (n) and Ω (n) are obtained from a local weighted average of the volumetric fluid velocity (respectively, rotational velocity) field over the region occupied by the particle in (respectively, Eq. ( 5))

V (n) = u(x,t) ∆(x -Y (n) ) d 3 x, (4)
Ω (n) = 1 2 [∇ × u(x,t)] ∆ ′ (x -Y (n) ) d 3 x. (5)
Isotropic Gaussian envelopes are used for spherical particles. The respective width of these Gaussian envelopes (σ and σ ′ ) is set with respect to the particle radius in such a way that the Stokes settling velocity and the hydrodynamic perturbation generated by a particle in a shear flow are both accurately reproduced at low Reynolds numbers. [START_REF] Lomholt | Force-coupling method for particulate two-phase flow: Stokes flow[END_REF] This ensures, among other results, that the particle velocity matches the Stokes settling velocity under low particulate Reynolds number and captures the Faxén correction for the motion in a non-uniform flow. [START_REF] Maxey | Localized force representations for particles sedimenting in Stokes flow[END_REF] Particle trajectories are then obtained from numerical integration of dY

(n) dt = V (n) . (6)
The key point of this modeling approach is that it allows the calculation of the hydrodynamic interactions with a moderate computational cost. In order to capture correctly the dynamics of dilute suspension flows, four grid points per particle radius are usually used when the monopole force is not zero, and in the case where only dipole forcing is relevant, three grid points per particle radius are sufficient.

When the particle volume fraction reaches a few percent, pairwise short-range hydrodynamic interactions and rigid body contact forces need to be accounted for. This is done by adding a monopole term following Ref. 2, which prevents particles from overlapping. The repulsion contact force between a pair of particles (α) and ( β) is written as a function of the relative position vector x α β = Y (α) -Y (β) and the distance r α β = x α β . If r α β < R ref (cut-off length scale of the repulsive barrier), then the contact force acting on particle α due to particle β, scaled by F ref , reads

F α β c = F ref 2a        R 2 ref -r α β 2 R 2 ref -4a 2        2 x α β . (7)
Otherwise the contact force is zero. In the present simulations, the value of the scaling factor F ref is chosen so that the number of overlapping particles is less than one percent of the total particle number. The dynamics of the suspension flows shown in this work is not sensitive to the value of the force scale in a rather large range of variation of F ref .

In the absence of external forces (no gravity effect for neutrally buoyant particles), the coupling between the particles and the carrier flow occurs exclusively from the force dipole term which is mainly related to the local flow strain rate. Accordingly, the method has been validated in finite Reynolds number flow situations where the force dipole is the only contribution to hydrodynamic interactions. Different tests have shown that the method is accurate under these conditions when the average concentration φ b < 20% and the particle Reynolds number Re p < 10.

A. Particle stresslet in a Couette flow

A first test consists in placing a freely rotating particle at the center of a Couette flow. The computation domain is a cubic box of size (20a) [START_REF] Asmolov | The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number[END_REF] where a is the particle radius. The mesh is uniform in all directions. The mesh resolution is three grid points per particle radius (which is the general recommendation for FCM dipole). The shear is imposed thanks to two walls moving with equal velocities in opposite directions (no-slip boundary conditions on the walls and other boundary conditions are periodic). Figure 1 shows the time evolution of the shear stresslet component scaled by µγa 3 for two particle Reynolds numbers Re p = γa 2 /ν = 10 -5 and 1. The time is scaled with the diffusion time t * = tν/a 2 = tγ/Re p . The curves of the shear stresslet component plotted for both Reynolds numbers collapse at small times t < 1. Then, each evolution converges to its specific value under steady state. The characteristic time scale a 2 /ν corresponding to the diffusion of the velocity perturbation due to the presence of the particle increases with Re p . This yields an increase of the convergence time. Typically, dimensionless time required to reach the complete convergence of the stresslet terms is around 2 for Stokes flow. When Re p varies from 10 -5 to 1, this corresponds to 100 and 2000 time steps, respectively, to achieve convergence.

The stresslet terms were computed for particle Reynolds numbers up to 10. Results are in very good agreement with those obtained by Ref. 28 from boundary fitted direct numerical simulations (see Figure 2). The amplitude of normal stress differences N 1 = G 11 -G 22 and N 2 = G 22 -G 33 is an increasing function of the particle Reynolds number. Therefore, it is suggesting that inertia promotes the non-Newtonian character of sheared suspension flows. The agreement of our FCM simulations is also good for those quantities.

B. Particle equilibrium position in a channel flow

The experiments of Ref. 35 highlighted the migration of neutrally buoyant spherical particles towards equilibrium positions in laminar pipe flows. Their study showed that the equilibrium position was located at 0.62R relatively to the tube center (where R is the tube radius) for Re < 30. When the Reynolds number increases, the equilibrium position shifts towards the wall. The inertial migration was later explained by different theories that quantified the lift force responsible of the particle cross-streamline motion. It will be discussed with further details in Section III.

The Segré-Silberberg effect is an appropriate test for the numerical method with a freely moving particle in channel flow. The test consists of initially introducing a particle in the flow, anywhere except in the channel central plane which is an equilibrium position unstable though, and calculating its trajectory in time within a wide range of Reynolds numbers [100-1500]. The time evolution of the trajectories is displayed in Figure 3. In this figure, the particle is located at the beginning of the simulation at a distance H/4 from the closest wall. Then, the particle migrates towards the channel wall and reaches an equilibrium position. The larger the Reynolds number, the closer to the wall is the equilibrium position. Moreover, the migration velocity is an increasing function of the Reynolds number.

The evolution of the dimensionless particle equilibrium position 2 y eq /H with the Reynolds number is plotted in Figure 3. Most of the numerical results relative to this test are obtained with the channel-to-particle size ratio H/a = 32 and with 3 grid points per particle radius. They are in good agreement with the prediction of Ref. 3 (with a relative error lower than 1.5%). Note that the situation of Re p < 1 (used in the work of Asmolov) corresponds to a channel Reynolds number less than 350. The agreement of the numerical results with the theoretical prediction of Asmolov suggests that, for the considered range of Reynolds numbers, the particle Reynolds number has small effect on the migration mechanism. Similar results obtained with smaller particles (H/a = 40) and a spatial resolution of three points per particle radius did not reveal any influence on the particle equilibrium position. Moreover, we tested the influence of the domain size and grid resolution on the numerical results. Increasing the streamwise domain length from 28.8 to 57.6a (or equivalently from 0.9H to 1.8H) did not lead to any significant modification on the particle migration velocity or trajectory (at Re = 100 and 1500). When the grid resolution is increased from 3 to 4 grid points per particle radius, the particle equilibrium position is slightly changed. The largest discrepancy is observed at Re = 1500 where the distance between the channel wall and the particle center is 1.3% larger for the finest mesh grid.

III. HOMOGENEOUS SUSPENSION FLOWS AT SHORT TIME SCALE

In order to analyze the macroscopic properties of a laminar homogeneous suspension flow in a channel at finite Reynolds numbers, we performed short time simulations that ensure suspension homogeneity. An initial random draw of particle positions is set and statistics are computed when stresslet components are converged. The simulation is stopped before lateral migration has led to significant displacement of particles. Ensemble averages are taken over several tests run at different random positions of the particles. However, the number of time iterations should be carefully chosen. On one hand, the average quantities in the present simulations are representative of a homogeneous suspension only when the number of time iterations n realized for each particle distribution is smaller than the time needed for the particles to migrate. On another hand, n should also insure the convergence of the second order terms (dipole) of the force coupling method (see Figure 1). This type of simulation allows to uncouple the effect of particle Reynolds number and of segregation due to particle inertial migration on the effective properties of the suspension flow.

Because particles undergo cross-streamline migration at finite Reynolds numbers, we were constrained to fix the number of iterations to 100. In Figure 1, symbols indicate the value of the shear stresslet terms reached after n = 100 time steps. With this choice of n, the stresslet is slightly overestimated at Re p = 1 (around 1%) whereas the resulting error at Re p = 10 -5 is negligible.

No-slip boundary conditions are imposed at the walls ( y = 0 and 32a). Periodic boundary conditions are imposed in the flow (x) and spanwise (z) directions. The channel size in these directions is carefully chosen to limit the interaction of the particles with their images through the periodic boundaries. The dimensions of the computational domain are (20a × 32a × 20a). In the Couette flow configuration, the walls move in opposite directions. In Poiseuille flow configuration, the flow is driven by a constant pressure-drop. In both cases, the average flow shear rate is noted γ . The time step and the number of iterations are, respectively, set to 0.01 γ -1 and 100. This simulation time allows capturing fairly well the particle-flow interactions. Also it corresponds on average to 4 and 7 particle encounters in the Couette flow for φ b = 1% and 5%, respectively.

A. Flow parameters and statistical quantities

The parameters used for the homogeneous suspension flow simulations are summarized in Table I. The macroscopic behavior of the suspension flow is characterized through two major quantities. First, the stresslet terms obtained from the FCM are used to characterize the suspension rheology because they are the main contribution to the suspension stress. Following Ref. 5, in the absence of any external torque applied on the rigid particles, the solid phase stress tensor of a suspension is obtained from the integral of the stresslet, variation of particles momentum, and Reynolds stresses over the suspension volume (first, second, and third terms on the right hand side of

Σ (p) i j = 1 V A 0 1 2 σ i j x j + σ j k x i n k dA - 1 V V 0 ρ f ′ i x j dV - 1 V ρu ′ i u ′ j dV, ( 8
)
where V is the volume of the fluid domain containing particles of external surfaces A 0 and volume V 0 . n k is a unit vector normal to A 0 pointing outwards. f ′ i is the local particle acceleration relative to the average value of the acceleration in V . It includes the non-hydrodynamic forces such as interparticle forces occurring when particles are in contact. u ′ i is the particle-induced velocity perturbation, and the calculation of the last contribution is straightforward. The suspension properties are homogeneous in the flow and spanwise directions. However, the shear gradient and the walls induce stress gradients in the wall-normal direction. Stresslet profiles in that direction will be therefore investigated. To calculate the average stresslet contribution to the particle-induced stress, averages are performed on slabs of volume V sl parallel to the walls, using the following summation:

S sl i j = N sl n=1 ζ (n) sl G (n) i j V sl . (9)
N sl stands for the number of particles present in the slab weighted by the volumetric percentage ζ (n) sl of the particle (n) belonging to this slab. The most telling example to illustrate the increase of the suspension stress due to the particles is the suspension viscosity. It is well known that a collection of particles enhances the local shear viscosity as a function of the solid volume fraction. For dilute suspensions of non-Brownian rigid spheres under unconfined Stokes flow, the Einstein's expression [START_REF] Einstein | Investigations on the Theory of the Brownian Movement[END_REF] of the effective shear viscosity of the suspension varies linearly with the concentration µ eff = µ(1 + 5 2 φ) in the dilute regime, while the extra term in the viscosity expression comes from the particle contribution. In Stokes flow, the viscosity enhancement comes exclusively from the stresslet. In order to obtain the increase of suspension viscosity due to inertia, the particulate contribution to the shear stress S sl 12 is calculated using the term G (n) 12 obtained by FCM simulations at finite inertia. Scaling S sl 12 by the fluid viscosity µ and the local shear rate γ sl leads directly to the increase of the suspension viscosity µ eff due to the stresslet contribution

µ eff -µ µ = S sl 12 µ γ sl . (10)
Second, the Eulerian slip velocity profile of the solid phase with respect to the fluid phase is defined by

u s = u p -U x, z χ p x, z χ p x, z , (11)
where χ p is the solid phase indicator (equal to 1 if the Eulerian grid point is located inside the particle volume and to 0 otherwise). The brackets indicate the average in time and in the homogeneous directions (x, z). Hence, the slip velocity is defined as the difference between the particle velocity field defined as

u p (x) = V (n) + (x -Y (n) ) × Ω (n) if x -Y (n)
< a and 0 otherwise, and the average suspension flow velocity U x, z .

B. Plane Couette flow

The increase of the flow viscosity induced by the particles is displayed in Figure 4 at different particle Reynolds numbers and concentration using S sl 12 /(µγ sl φ b ). Compared to the case of an isolated sphere, the simulations performed with the FCM account for wall effects and multi-body hydrodynamic interactions. Results for small volume fraction (φ b = 1%) are represented with dashed and solid lines for Re p = 10 -6 (which will be noted here Re p = 0) and Re p = 1, respectively. Results for larger fraction (φ b = 5%) are represented in Figure 4 The two vertical thin solid lines are linearly extrapolated from the DNS simulation of a single particle in the center of a Couette flow, [START_REF] Mikulencak | Stationary shear flow around fixed and free bodies at finite Reynolds number[END_REF] at Re p = 0 and 1 from left to right. The horizontal thin lines define a distance of one particle radius from the walls.

In the central region of the domain, the particle stresslet at φ b = 1%, Re p = 0 almost matches that for an isolate particle in unbounded Stokes flow (the relative difference is 2%). At φ b = 5% and Re p = 0, the multi-body hydrodynamic interactions add a quadratic term in the suspension viscosity µ eff = µ(1

+ 5 2 φ + K φ 2 )
, where K is a coefficient between 5 and 7.5 depending on the local microstructure. [START_REF] Batchelor | Stress system in a suspension of force-free particles[END_REF] The FCM results lead to K ≈ 5. One can note that the stresslet profiles are flat in the center of the channel. Near the walls, the stress due to the rigidity constraint in the particle volume is increased due to the interaction with the rigid wall, independently from the concentration, and yields modification of the stress distribution on the particle surface (this holds true for a single particle in a linear flow near a wall at small and finite Re p ). [START_REF] Magnaudet | Small inertial effects on a spherical bubble, drop or particle moving near a wall in a time-dependent linear flow[END_REF] At Re p = 1, the stress is increased for both concentrations (φ b = 1% and 5%) compared to the case with negligible flow inertia, and the wall effect on the particle stress acts on a wider region. The stresslet is found constant only in a narrow band in the channel core (12 < y/a < 18). Considering the value of the stresslet only in that band, we find that at Re p = 1 the contribution of the stresslet (divided by the concentration) is increased with respect to Re p = 0 by a factor of 1.036 and 1.023, respectively, for the concentrations 1% and 5%. This increase falls in the range of other predictions found in the recent literature although there are to some differences in the cases reported. Reference 32 found that in 2D (suspension of discs) the stresslet at Re p = 1 increases by 1.15 with respect to the case Re p = 0 in a dilute suspension. Reference 38 predicted, by considering the total shear stress, that the increase of viscosity is (µ eff -µ)/µ = (2.5 + 0.6Re 3/2 p )φ in an unbounded dilute sheared suspension flow laden with rigid spheres, leading to an increase of a factor 1.24 when Re p is increased from 0 to 1, independently of the concentration. As for Ref. 17 who used the Lattice-Boltzmann method to simulate a suspension flow in a plane Couette, they obtained a factor 1.2 at φ b = 5% when Re p is increased from 0 to 0.1.

In pure shear flow, the torque-free particles are experiencing a constant shear rate that makes them rotate in the spanwise direction. When inertia is neglected (Re p ≪ 1), the rotation rate is equal to the mean flow vorticity, Ω 0 = -γ/2. When flow inertia increases, the rotation rate deviates from the local flow rotation rate and (Ω -Ω 0 )/γ scales with Re 3/2 p . [START_REF] Lin | Simple shear flow round a rigid sphere: Inertial effects and suspension rheology[END_REF] In the presence of a wall, the velocity perturbation is modified by the no-slip boundary condition. Thus, the particle exhibits a slip velocity which is similar to a Faxén correction because of the symmetry breaking of the velocity on the particle surface (the slip velocity is the highest when a particle is very close to the wall [START_REF] Ho | Inertial migration of rigid spheres in two-dimensional unidirectional flows[END_REF] ). In the center of the Couette flow, this interaction is symmetric and particles exhibit no slip in the flow direction.

Profiles of the slip velocity components in the flow and wall-normal directions are displayed in Figure 5 are negative in the upper part (where the fluid velocity is positive). Note that all the slip components are linear across the Couette gap, and they are almost equal suggesting an average slip motion angle of -π/4. As for the multi-body hydrodynamic interactions, they seem to have weak impact in the range of concentrations considered. Indeed the velocity profiles closely match the migration velocity of a single particle when placed at different positions across the Couette gap (filled squares).

The hydrodynamic repulsion by the walls at finite inertia induces a finite lift force that drives the neutrally buoyant particles towards the center. The theoretical predictions available in the literature on Couette flows with neutrally buoyant particles assume that the wall lies within the inner region of the flow perturbation around the sphere, [START_REF] Ho | Inertial migration of rigid spheres in two-dimensional unidirectional flows[END_REF][START_REF] Vasseur | The lateral migration of spherical particles sedimenting in a stagnant bounded fluid[END_REF] so that the dominant inertial contribution comes from a regular perturbation of the Stokes solution. Thus, they are rather valid at small particle Reynolds numbers, and seem to over-estimate the lift force at Re p = O(1) (open circles and plus symbols in Figure 5(b)). When inertia becomes important at the particle scale, the leading order inertial effects result from the outer expansion of the perturbation and the wall likely falls in the outer region. The evaluation of these terms requires the solution of a singular perturbation problem which can be obtained using the matched asymptotic expansions technique (see, for instance, Ref. 34 Let us first consider the particle contribution to the shear stress in the suspension. In a steady Poiseuille flow, the profile of the shear stress is linear. Also the strain rate (∂U/∂ y) is linear across the channel when the viscosity is constant. Figure 6 presents the plot of the stresslet contribution to the shear stress (divided by the fluid viscosity) as a function of the local shear rate γ sl . The curves are linear for the four test cases considered here (linear interpolation in solid and dashed lines). Therefore, the suspension response to the shear resembles to a Newtonian fluid up to Re p = 3 (maximum particle Reynolds number at Re = 500).

The slope of the stress-strain curves is clearly increased with the concentration. The increase of the suspension viscosity (µ eff -µ)/µφ can be calculated using this slope divided by the cross section average shear rate, and not the local shear rate to avoid divergence near the channel center, i.e., S sl 12 y /(µ γ sl y φ b ), where y indicates the average across the channel height. These values are reported in Table II. At low inertia and low concentration, the viscosity increase is close to the Einstein's limit, whereas it is slightly higher at φ b = 5% due to hydrodynamic interactions. At finite inertia, the viscosity increases slightly less than the value predicted in Ref. 38 (3.6 for Re p = 1.5), as in the case of Couette flow.

Second, concerning the particle motion in the flow direction, whatever the flow regime, a neutrally buoyant finite-size spherical particle placed in a quadratic flow lags the local flow. In Stokes flow (Re ≪ 1), Faxén law relates the particle velocity to the force and torque and accounts for the non-uniform flow velocity gradient at the particle scale. The slip velocity calculated from Faxén law, without accounting for the presence of the channel walls, is equal to -4 3 U m (a/H) 2 (where U m is the velocity at the channel center). One can find, for example, in Ref. 16, the correction due to the walls in Stokes flow which is dependent on particle position in the channel. Profiles of the slip velocity in the flow direction, normalized by 2U m a/H (which is the average shear velocity at the particle scale), are plotted in Figure 7(a). At Re = 10 -3 , the slip velocity is constant across the channel width (except near the walls). At low Reynolds number, the profiles are close to the Faxén estimate in the limit of unbounded parabolic flow, and to its extension to bounded channel flows as written in Ref. 16. The small discrepancy between the FCM results and the theoretical predictions is due to the accuracy of the numerical representation of the multipole expansion, where higher order terms than the dipole are not accounted for. Terms like the quadrupole start to be important when the gradient of the strain rate in the particle volume is not negligible (corresponding to larger particle-to-channel size ratio). At Re = 500, the magnitude of the slip velocity of the particulate phase increases. For both Stokes and inertial regimes, the streamwise slip velocity appears to be independent of the concentration when φ b < 5%. The profiles of the cross-streamline migration velocity (normal to the channel walls) scaled by 2U m a/H are shown in Figure 7(b). Only half of the channel height is shown for clarity, knowing that the profiles are perfectly symmetric. At negligible inertia, the migration velocity is vanishingly small for both concentrations. Similarly to Couette flow, when fluid inertia is important, the neutrally buoyant particles experience a lift force normal to the walls. The difference with a linear shear flow is that, away from the wall, the lift force is induced by the interaction between the stresslet and the curvature of the bulk velocity profile for finite particle Reynolds number. It is directed towards the wall for particles located around the channel center, and away from it for the particles close to the channel wall (Segré-Silberberg effect).

In order to validate the numerical computation of the migration velocity in the channel flow, results obtained with a single particle located at different positions across the channel at Re = 133 (open squares in Figure 7(b)) are compared with the theoretical prediction (triangles) taken from Ref. 26 at the same Re (i.e., Re = 200 in their paper). One can note that the FCM results agree reasonably well with the theoretical prediction valid in the limit Re = O(10-1000). Also the location of the equilibrium position, where the wall-normal velocity cancels, is correctly captured. Deviation between numerical and theoretical values is observed only very close to the wall.

In Figure 7(b), the wall-normal velocity profiles obtained for φ b = 1% and 5% at Re = 500 well match that of a single particle at same Re (filled squares). The velocity of a single particle cancels at a distance of 2.4a from the wall, in agreement with the theoretical prediction of Ref. 3. At φ b = 1%, the velocity cancels at the same position but at φ b = 5%, the lateral migration velocity decreases without reaching zero (thick solid line). In this case, hydrodynamic interactions between particles are probably hindering the wall repulsion. Also the hydrodynamic interaction of two particles almost touching the wall is probably not accurately calculated by the FCM.

IV. DEVELOPING TWO-PHASE FLOW

In this section, the long-time behavior of the suspension flow is analyzed, using the same channel-to-particle aspect ratio (H/a = 32). The domain dimensions are (100 × 32 × 80)a 3 and contain 320 × 100 × 256 cubic cells. It is large enough to minimize the effect of particle interaction through the periodic boundaries. Initially, the particles are randomly seeded in the domain corresponding to a homogeneous distribution. Simulations with different solid fractions φ b = [0.5%-5%] and different Reynolds numbers Re = [1000-1400] are performed at constant average velocity by adjusting the longitudinal pressure drop. The flow regime remains laminar at all times and no transition to turbulence is observed.

A. Stratification of the suspension

Over time, the freely moving particles migrate towards the channel wall leading to a stratification of the suspension, as illustrated in Figure 8. At φ b = 1%, a single layer of particles forms near each wall of the channel, very close to the equilibrium position of a single particle at the same Reynolds number. Particle alignment in the flow direction is observed (see Figure 8(a)). This behavior is similar to the formation of particle trains detected in pipe flows in the experiments of Ref. 25 and later in the simulations of Ref. 37. The particle alignment depends on Re p according to Ref. 25 and is likely due to hydrodynamic interactions of the finite-inertia particle-induced disturbance in the wall shear flow. At φ b = 2.5% and 5%, the suspension splits into three layers: the center of the channel becomes depleted from particles while they migrate towards the walls due to Segré-Sielberberg effect. Wall regions become concentrated while particle layers develop (Figure 8(b)).

The profiles of the solid volume fraction across the channel are displayed in Figure 9 after the migration is completed (when t is approximately equal to 50H/ U ). These profiles correspond to the volume fraction of particles averaged in space along the streamwise (x) and vorticity direction (z), and in time over 15 time units (H/ U ). The peak of concentration is located at approximately one particle radius from the channel wall. The local minimum next to the peak is due to excluded volume effect (because of the finite size of the particles). The height of the peaks and the level of the concentration in the depleted core zone are increasing function of the global concentration. Note that even though there is enough room for all particles to be trapped at the single particle equilibrium position (as it is the case at φ b = 1%), most of the particles remain freely suspended throughout the entire channel. The development of these profiles in the wall normal direction suggests that the inertial migration of the particles is balanced by concentration gradient and shear-induced particle diffusion oriented towards the channel center. The mass balance equation applied on the solid phase, at steady state, leads to equilibrium between inertial migration and shear-induced diffusion fluxes, the equilibrium being independent of the Reynolds number as suggested by the concentration profiles in Figure 9(b). However the wall-normal particle migration velocity estimated from homogeneous suspension simulations (Figure 7(b) extended to higher Re) is orders of magnitude higher than the opposite shear-induced particle diffusion velocity based on a low Reynolds number diffusion coefficient. [START_REF] Abbas | Dynamics of bidisperse suspensions under Stokes flows: Linear shear flow and sedimentation[END_REF][START_REF] Drazer | Microstructure and velocity fluctuations in sheared suspensions[END_REF][START_REF] Yeo | Dynamics and rheology of concentrated finite-Reynolds-number suspensions in a homogeneous shear flow[END_REF] This leads us to the conclusion that a simple FIG. 9. Solid volume fraction profiles of a laminar channel flow after the particle migration towards the wall is completed. (a) Re = 1400 and φ b = 1 (thin line), 2.5 (dashed), and 5% (solid thick line). (b) φ b = 5% and Re = 1400 (dashed-dotted), 1200 (light gray -green online), and 1000 (dark gray -blue online). mass balance on the particulate phase is not sufficient to predict the particle concentration profiles established in the channel flow.

B. Agitation of the particulate phase

Information on particle fluctuating motion is particularly interesting for particle-phase stress prediction across the channel flow (this will be explained further in this section). Particle fluctuating motion is not isotropic. Velocity fluctuations in the wall-normal and spanwise directions are very close whereas streamwise component is larger than both by a factor close to 4. We call particle agitation the isotropic part of particle velocity fluctuations taken here as 0.5(u ′2 p, y + u ′2 p, z ). The bar symbol denotes the time and phase average operation on the particles in the x and z directions, i.e., ψ p = ψ p χ p / χ p = ψ p χ p /φ, where χ p is the local solid volume fraction and χ p = φ is the average wall-normal dependent concentration. The explicit y dependence of φ is omitted in this section for brevity.

The particle agitation profiles in the wall-normal direction are plotted in Figure 10(a) for φ b = 2.5% and 5% and Re = 1400. They peak at a distance from the wall slightly larger than one particle diameter and decrease towards the channel center, following the trend of the concentration profiles. The agitation profiles at φ b = 5% seem, however, to be slightly dependent upon flow Reynolds number (corresponding to Re = 1200 and Re = 1000, respectively). In the channel center, the particle agitation and concentration are both non-zero. Therefore, it is suggested that the diffusive flux term plays a non-negligible influence in the budget of the particle phase agitation. The role of the diffusive flux was already pointed by Ref. 30 who first performed the simulation of (quasi 2D) concentrated suspension flows in a channel at zero Reynolds number. However, concentration FIG. 10. Profiles of the dimensionless particle agitation T * = T /(γ m a) 2 . Black lines are results obtained from the simulations. The thick and thin line styles refer to φ b = 5% and 2.5%, respectively. In plot (a), solid, dashed, and dashed-dotted lines are for Re = 1400, 1200, and 1000, respectively. In plot (b), the dashed (blue online) lines correspond to Eq. ( 20) with α ≈ 2000. The dotted (red online) lines result from a prediction based on pure shear flow prediction (Eq. ( 21)).

profiles of suspension channel flows at low Reynolds number are quite different since they peak at the channel center, [START_REF] Frank | Particle migration in pressure-driven flow of a Brownian suspension[END_REF][START_REF] Nott | Pressure-driven flow of suspensions: Simulation and theory[END_REF] whereas the concentration profile peaks near the wall in the present case.

We examine now the transport equations of the momentum (in the wall-normal direction) and agitation of the particle phase at steady state. The mean wall-normal component of particle velocity being negligible, the momentum balance in that direction reduces to a balance between the divergence of the particle normal stress and the lift force due to Segré-Silberberg effect (noted F ss y ). The momentum balance on the solid phase can be written following Eulerian two-phase formulation as

ρ p d dy φΣ p, y y + φF ss y = 0, (12)
where Σ p = P p -µ p ∇U p + ∇ t U p is the particle phase stress tensor, which can be split into two contributions, namely, the particle phase pressure and the viscous shear stress. This is true in the situation where the particle phase velocity field is divergence free. The incompressibility condition is satisfied in the channel flow, thanks to the negligible particle mean motion in the wall-normal and spanwise directions (which implies that ∇. φU p = 0 reduces to ∇.U p = 0). The general form of this constitutive relation for the particle stress tensor is valid for both non-Brownian Stokesian suspensions (see Ref. 30) and gas-solid suspensions, [START_REF] Boelle | Second-order prediction of the particle-phase stress tensor of inelastic spheres in simple shear dense suspensions[END_REF] where P p is the granular pressure and µ p the granular viscosity. In the frame of kinetic theory of granular flows (KT-GFs) at large Stokes number, P p and µ p both result from a kinetic contribution (normal Reynolds stress) and a collisional contribution (see Ref. 7). The granular pressure P p is related to the particle Reynolds stress, i.e., in average notation, P p = ρ p φT where T is the granular temperature equal to one third of the mean particle fluctuating energy T = [START_REF] Abbas | Shear-induced self-diffusion of inertial particles in a viscous fluid[END_REF] 3 T ii , where φT i j = u ′ p,i u ′ p, j χ p . In the present study, particle Reynolds and Stokes numbers are of the order of unity and maximum solid volumetric fraction near the wall does not exceed 20%. The Stokes number, being the ratio of particle relaxation to fluid flow characteristic time scale, is defined as St = γ m τ p = 2 9 ρ p /ρ f Re p , where γ m is the average shear rate over the channel flow. An alternative definition of the Stokes number would be St = 1 3 ρ p /ρ f Re p accounting for the added-mass effect (but this would not lead to any significant variation for neutrally buoyant particles). In that regime (low Reynolds and Stokes numbers and moderate concentration), particle soft collisions are induced by the shear (T < (γa) 2 ), corresponding to the so-called quenched state in KTGF as first introduced by Ref. 39, and later extended by Ref. 31. In the KTGF formulation, granular viscosity is due to particle agitation and can be calculated with collision integrals. An expression of particle phase viscosity as a function of concentration can be deduced from the calculation of particle stress tensor in homogeneous pure shear flow in the quenched state regime, following Ref. 39,

µ p = µ f φ 2 b(St), (13)
where b(St) = 64 35π St 3 1 + 9π 16 St -1 . Although there is no simple analytical expression for F ss, y in Eq. ( 12), it is clearly the only term that may allow at steady state a strong correlation between the concentration and the agitation profiles of the particle phase. Indeed, in the limit of zero Reynolds number, this term would vanish and Eq. ( 12) would naturally reduce to an inverse relation between these two variables. [START_REF] Frank | Particle migration in pressure-driven flow of a Brownian suspension[END_REF][START_REF] Nott | Pressure-driven flow of suspensions: Simulation and theory[END_REF] At steady state, the macroscopic balance of particle agitation

T is 0 = -∇.Q p -Σ p : ∇U p -ε p . (14)
The first term of this equation is the divergence of the flux of fluctuating kinetic energy, which is usually closed as a linear function of ∇T,

Q p = -κ p (φ)∇T, (15)
where κ p is the conductivity, which in rarefied gas (i.e., dilute hard sphere systems) has the same order of magnitude as particle phase viscosity, µ p . Then, the conductivity is assumed to be proportional to the granular viscosity (the pre-factor is the inverse of a Prandtl number) κ p = α µ p . The second term of Eq. ( 14) is the production by the mean velocity gradient. The third one is the dissipation term. Note that this term also includes the contribution of the fluctuation of the power of Segré-Silberberg force. As the dispersed phase fraction and the Stokes number considered here are moderate, we assume that the characteristic time scale of the dissipation rate is mainly controlled by the drag rather than the collisions,

ε p = 3ρ p φλ(St) T τ p , (16)
where τ p is the response time scale of the particles (τ p = 1 3 ρ p µ f a 2 ). λ(St) is a function of Stokes number exclusively, accounting for the contribution of shear-induced collisions to the agitation. In homogeneous shear flow, λ(St) can be calculated from the equilibrium between the production (second term of Eq. ( 14)) and dissipation (third term), the granular temperature given from the quenched state theory, [START_REF] Parmentier | Transition models from the quenched to ignited states for flows of inertial particles suspended in a simple sheared viscous fluid[END_REF][START_REF] Tsao | Simple shear flows of dilute gas-solid suspensions[END_REF]

T * = T (γa) 2 = φ f (St), (17)
where

f (St) = 128 945π St 3 1 + 9π 16 St -1 + 9 2 St -2 . Consequently, λ(St) = 1 + 9π 16 St -1 1 + 9π 16 St -1 + 9 2 St -2 . ( 18
)
The particle agitation in the channel flow being mainly driven by the shear, we verify that the quenched state closure laws ensure the observed level of T, assuming constant average shear rate across the channel. When St = 1, Eq. ( 17) gives T * ≈ φ π ≈ 0.31φ, which is a very close result to that predicted by Lagrangian simulations of Ref. 1 at same Stokes (see also Ref. 31). However, when St < 1, the prefactor f (St) in Eq. ( 17) provides much smaller values than those predicted by the simulations. Based on FCM simulations, Ref. 42 has calculated the particle Reynolds stress normal components in homogeneous concentrated sheared suspensions, at different concentrations and particle Reynolds number. Compared to shear flow at low Re p , particle agitation slightly decreases as the particle Reynolds number increases. With neutrally buoyant particles at φ = 20%, they obtained T * ≈ 0.24φ when Re p = 1 (St = 0.33), and T * ≈ 0.3φ when Re p < 10 -2 (St ≈ 3 × 10 -3 ). In the limit of dilute suspension and zero Re p , Stokesian dynamics based on simulations from Ref. 12 gives T * ≈ 0.34φ. Instead Eq. ( 17) leads to T * ≈ 4 × 10 -3 φ when St = 0.2 and T * ≈ 2 × 10 -3 φ when St = 10 -2 , which is two orders of magnitude smaller. Clearly, quenched state closure laws (Eq. ( 13)) are not applicable when St < 1, mainly because it does not account for the hydrodynamic interactions occurring during particle encounters.

As a consequence, given the range of Stokes and particle Reynolds number (0 < St, Re p < O(1)) and particle phase fraction (0 < φ < 0.2) investigated in the channel flow, a good approximation of particle phase viscosity in Eq. ( 19) will be µ p = 4µ f φ 2 and λ(St) = 1 ensuring that T * ≈ 0.3φ in pure shear flow case.

If we now write the macroscopic balance (Eq. ( 14)) in the wall-normal direction y, we obtain

d d y κ p (φ, St) dT dy + µ p (φ, St) dU p, y dy 2 -3ρ p φλ(St) T τ p = 0. ( 19
)
In this equation, the gradient of the particles velocity is the same as that of the fluid and is known in the whole channel height since the flow is parabolic (except at one particle diameter from the walls). In dimensionless form, Eq. ( 19) writes

d 2 T * d y * 2 - 27 32αǫ 2 T * φ = - 1 αǫ 2 (1 -y * ) 2 . ( 20
)
The star indicates dimensionless quantities such that y = y * (H/2), T = T * (γ m a) 2 where U p, y = U * f , y (γ m a). This model can be tested along the channel width y * when the calculated profile of the particle concentration φ( y * ) is introduced in Eq. (20). Neglecting the diffusive term, the temperature profile T * ( y * ) is given by

T * = 32 27 (1 -y * ) 2 φ( y * ). (21)
These profiles have been reported in Figure 10(b) at particle concentration of φ b = 2.5% and 5%. It can be seen that for a given concentration, Eq. ( 21) does not reproduce well the shape of the calculated temperature profiles neither the influence of particle concentration. The prediction can be improved with the diffusive term but in order to match the two profiles at 5%, the coefficient α between the viscosity and the conductivity must be taken equal to O(ǫ -2 ) ≈ 10 3 , which is unrealistic. Even when this high ratio is considered, a unique set of parameters (µ p , κ p , λ) does not allow capturing the evolution of the temperature profile for all bulk concentrations. As a conclusion, the temperature profiles in the wall normal direction do not correspond to a steady developed profile that can be described by the energy budget (Eq. ( 14)). The temperature and concentration profiles in y-direction probably result from other contributions involving interactions in the vorticity direction (z). These interactions have been highlighted, thanks to longer simulation times, and they are presented in Sec. IV C.

C. Secondary flows and interface instability

At φ b = 5%, after the suspension has been segregated into particle-laden layers, and a central fluid core, the separating boundary does not stay permanently parallel to the channel walls. Unexpectedly, the development of wavy-like patterns for the particle concentration in the spanwise direction is observed approximately 70 time units after the separation occurred. This pattern is clearly illustrated in Figure 11(a), where the particle positions projected onto the ( y, z) plane exhibit the shape of the transverse concentration waves (Figure 13(b)). The order of magnitude of the wavelength of this interface is half of the channel height. Note that the average velocity profile of the suspension flow does not deviate from the Poiseuille flow profile (see Figure 11(c)).

Reference 9 demonstrated that when two stratified fluid layers of different rheological properties are sheared, the interface experiences an instability in the vorticity direction. Considering the particular case of a plane Couette flow of two stratified fluids (Newtonian fluid and concentrated suspension at very small Re), Ref. 9 has shown that the instability of the interface is due to a jump in the second normal stress difference across the two layers. According to their analysis, a necessary condition for the instability to develop is that the jump (N 2N -N 2n N ) in the second normal stress difference between the Newtonian N 2N and the non-Newtonian N 2n N fluid at the interface be positive. The normal stress difference of the Newtonian fluid being null (N 2N = 0), the instability onset criterion reduces to N 2n N < 0. This configuration of two stratified layers of different rheological properties is comparable to the present case of study, even though the shear is not constant along the channel height.

While the flow free of particles in the central region is Newtonian, the layers laden with particles are likely non-Newtonian. The anisotropy of the normal stresses occurring in sheared suspensions at finite Reynolds numbers has been reported in many studies (see Refs. 17, 32, and 42). Therefore, the particle-laden zones exhibit non-zero normal stress differences. The main contributions to the second normal stress difference N 2 are due to the stresslet components (S sl 22 -S sl 33 ) and to the transverse velocity fluctuations difference (v ′2 -w ′2 ) (second and third term of Eq. ( 8)). The profiles of these contributions averaged over (x) and (z) are shown in Figure 12. Both the stresslet and inertial contributions are almost zero in the region free of particles, and they are of same order of magnitude in the particle-laden regions. These contributions sum up in the second normal stress difference, which turns out to be positive. It is not surprising to find positive N 2 in so far as the particle Reynolds number is Re p = O(1) and the maximum concentration reached near the wall does not exceed 20%. Reference 42 has shown that N 2 is positive in this range of parameters, and that the sign switch takes place at 20% ≤ φ ≤ 30%. Back to the definition of Ref. 9, the jump in the second normal stress difference ∆N 2 = (N 2N -N 2n N ) is negative in the present case and does not meet their criterion for the onset of the interface instability. vorticity, these authors have shown that the second normal stress difference N 2 (which is non-zero in a concentrated suspension even in Stokes flow) generates secondary flows whenever the flow cross section is not axisymmetric. In the particular case of a plane Poiseuille flow, an initially imposed concentration gradient in the channel cross section triggers the occurrence of dune-like patterns similar to those observed in the present work. The authors ascribe the origin of the secondary circulating flows to spanwise gradients in the curvature of the flow stream-surfaces, which is not the case since the flow is periodic in the spanwise direction.

The contour plot of the streamwise velocity fluctuations (averaged over x) is displayed in Figure 13(a). It can be observed that in the zones laden with particles, the mixture flows with negative velocity with respect to the average local Poiseuille flow whereas the fluid velocity in the depleted zones is positive. This can be explained by momentum conservation arguments based upon suspension effective viscosity fluctuations, where the more (respectively, less) viscous or concentrated regions have smaller (respectively, larger) velocities when driven by the same pressure drop along the longitudinal direction. The contours of the spanwise shear rate shown in Figure 13(a) reveal that particles are submitted to a non-uniform shear field in the plane (x, z) parallel to the channel wall combined with that of the streamwise plane (x, y), and therefore they migrate, in the (x, z) plane direction, from faster flow regions to slower ones. It can be seen in this figure that the amplitude of variation of this gradient fluctuation in a plane parallel to the walls is of the order of 1 over a distance z/a which is of order of 15. This exactly compares with the variation of the mean streamwise velocity gradient in the wall normal direction. The particle Reynolds number based on this spanwise shear rate varies between 0 and 2 in Figure 13(a). Therefore, the existence of an inertial lift force in the flow-spanwise plane, from faster (unladen) to slower (particle laden) flow regions, can be thought of the driving force of the instability. As the origin of the velocity gradient fluctuation ∂ (∂u ′ /∂z) /∂z is due to concentration fluctuations in z direction, as ∂u ′ /∂z develops, the migration of particles from lower to higher velocity gradient zones is increasing these concentration fluctuations and consequently ∂u ′ /∂z. Therefore, the development of the spanwise instability is a self-sustained process (the character of spanwise concentration fluctuation is reinforcing). This is confirmed in Figure 13(b) where the contour plot of the streamwise flow vorticity is reported. It is averaged in the streamwise direction and over time ([130-200]H/ U ) after the dune-like patterns took place. In Figure 13(c), this map is superimposed to Figure 11(a) representing the projected positions of particles in the channel cross section. It can be seen that each individual dune pattern is associated with a single pair of counter-rotating vortices, the size of which nearly fits that of the dune. This figure clearly illustrates the feeding mechanism of the dune from both sides of it (which correspond to the depleted zones). These secondary flows remix the particles in the channel cross section. In turbulent flows, the development of secondary flows is known to originate from the anisotropy of the Reynolds stress tensor in a pipe cross section (see, for instance, the recent work of Ref. 6). However, in the present case, their amplitude relatively to the bulk flow velocity is much smaller than in turbulent flows. Here, the anisotropy of the Reynolds normal stress components arises from the local concentration gradients, combined to localized forcing at the particle positions.

The wall normal migration induced by each pair of vortical structures is pointing to the opposite direction of the Segré-Silberberg induced migration in the (x, y) plane. Therefore, the development of a steady concentration profile along the channel height is resulting from an equilibrium between Segré-Silberberg migration due to the interaction between the particles and the Poiseuille flow on one hand, and on the other hand the same mechanism due to the development of concentration waves in the spanwise plane (such as illustrated by the scheme of Figure 13(c)). Note that according to such a mechanism, the development of both concentration and temperature profiles is entirely governed by inertial effects.

The observation of these instabilities has never been reported and it is likely to be limited in terms of concentration and particle Reynolds number. At φ b = 1%, the particles form a single particle layer close to the wall. At 5%, the dune height is of the order of half the channel height. Increasing the concentration of particles in the channel will tend to make the two wavy layers interact and damp the concentration gradients in the wall normal direction. It is likely that above 10%-15% of particle concentration, such patterns are not visible and a more homogeneous profile will be observed. Increasing the particle Reynolds number is also limited by the onset of turbulence which occurs at Reynolds number close to 1600 with H/a = 32. 20

V. CONCLUSION

We have studied the dynamics and migration of a suspension in a channel flow under laminar regime. Finite-size neutrally buoyant particles were modeled by the force coupling method coupled to direct numerical simulations at finite Reynolds numbers.

In the first part, we have analyzed homogeneous suspensions under plane Couette and Poiseuille flow configurations. In Couette flow, the particle induced stress increases with the Reynolds number. The effect of increasing the particulate concentration on the suspension stress is less important at finite than at low Reynolds numbers. This can be explained by a reduction of the length scale of the particle-induced flow perturbation when flow inertia is increased. The cross-stream and streamwise slip velocities increase with the Reynolds number and they are of the same order of magnitude at Re p = O(1). The main difference between the Couette and Poiseuille plane flows is that the particles tend to migrate towards the center because of wall repulsion in a Couette flow and towards the walls in the Poiseuille flow configuration.

In the second part, the developing two-phase flow structure in the channel is analyzed. First, the particles migrate towards the channel walls due to the Segré-Silberberg effect. Then, the suspension becomes stratified forming two concentrated layers close to the channel walls separated by a nearly pure fluid region in the core of the channel. Concentration and agitation profiles in the wall-normal direction are first analyzed through a steady two-phase flow model in which suspension transport properties as a function of particle concentration are derived from the kinetic theory of granular flows in the quenched state regime when St > 1 and from numerical simulations when St < 1. The model clearly fails to reproduce the agitation profile in the wall normal direction, suggesting the existence of secondary flows. This is confirmed by the simulations of flow dynamics at longer times which show the development of a particle concentration wavy structure in the spanwise plane (x, z). These structures are associated to counter-rotating vortices, which are thought to be driven by a secondary migration due to inertial lift force in this plane from depleted (less viscous and with higher velocity) zones towards laden (more viscous and with lower velocity) zones. In the laden zones, the balance between the outward particle wall-normal flux raising from the secondary flows and the flux towards the wall induced by the mean flow shear gradient (mainly due to Segré-Silberberg effect) leads to the establishment of the observed concentration and agitation profiles. Such a flow configuration is probably only observable in a narrow range of particle Reynolds number and moderate concentration (<10%). To the best of our knowledge, we have reported the first observations of such a peculiar suspension flow instability. Experimental evidences would be highly desirable to confirm our findings.

FIG. 1 .FIG. 2 .

 12 FIG. 1. Time evolution of the shear stresslet G 12 (scaled by µ γ a 3 ) of a particle in a pure shear flow. The time is normalized by the diffusion time scale a 2 /ν (µ and ν being the fluid dynamic and kinematic viscosity). The solid and dashed lines correspond to different flow inertia Re p = 1 and 10 -5 , respectively. The corresponding values obtained by the DNS of Ref. 28 are plotted as short horizontal lines. The blue plus and blue triangle symbols indicate the value of the stresslet after 100 time steps.

FIG. 3 .

 3 FIG. 3. Time evolution of the dimensionless lateral position of a single particle during the cross-streamline migration. The time is scaled by a/v 0 s where v 0 s is the velocity at the initial migration stage. v 0 s /(2U m a/H ) = 1200Re 0.535 is a good estimate of the initial particle velocity, where U m is the flow velocity at the channel center. The particle reaches different equilibrium positions at different Re ranging from 100 (bottom) to 1500 (top line). Inset: equilibrium lateral position and comparison with Asmolov's theoretical prediction (solid line).

  by symbols.

FIG. 4 .

 4 FIG.4. Particle contribution to the increase of the Couette flow viscosity. The thick dashed and solid lines correspond to φ b = 1% with Re p = 0 and Re p = 1, respectively. The cross and plus symbols correspond to φ b = 5% with Re p = 0 and Re p = 1, respectively. The two vertical thin solid lines are linearly extrapolated from the DNS simulation of a single particle in the center of a Couette flow,[START_REF] Mikulencak | Stationary shear flow around fixed and free bodies at finite Reynolds number[END_REF] at Re p = 0 and 1 from left to right. The horizontal thin lines define a distance of one particle radius from the walls.

  FIG. 5. Profiles of the solid phase slip velocity at Re p = 1 scaled by γa. (a) and (b) show the streamwise and wall-normal components, respectively. The line style stands for the particle concentration (solid line for φ b = 5%, dashed line for φ b = 1%). The plus and circle symbols are the analytical prediction of the wall-normal component obtained by Refs. 40 and 16, respectively. Filled square symbols correspond to the simulation with a single particle. Oscillations in the streamwise slip profile at φ b = 1% are due to the limited number of random particle realizations (1000). Grey (blue online) layers define a distance of one particle radius from the walls.

  for Re = O(1-100) and Ref. 3 for Re = O(1-1000) in channel flow).C. Pressure-driven channel flowSimilar simulations were carried on homogeneous suspensions subject to Poiseuille flow at Re = 10 -3 and Re = 500, with two different concentrations φ b = 1% and 5%.

FIG. 6 .

 6 FIG.6. Stresslet-shear rate diagram in Poiseuille flow. Square and plus symbols are, respectively, for Re = 500 and 10 -3 . Dark and light gray (online blue and red) colors are for φ b = 1% and 5%, respectively. The lines are obtained from linear fitting for each case.

FIG. 7 .

 7 FIG. 7. (a) Slip and (b) migration velocity profiles in Poiseuille flow scaled by 2U m (a/H ). Solid and dashed lines are for φ b = 5% and 1%, respectively. Thick and thin lines are for Re = 500 and 10 -3 , respectively. The vertical line in (a) is calculated from the Faxén law (4U m /3)(a/H ) 2 , and the stars are from Ref. 16 at Re = 0. In (b), the filled squares obtained using a single particle placed at different y positions at Re = 500. The open squares are the numerical results for Re = 133, and the open triangles represent the prediction based on the matched asymptotic expansion (taken from Ref. 26 in channel flow). Grey (blue online) layers define a distance of one particle radius from the walls.

FIG. 8 .

 8 FIG. 8. Particle positions in the channel at the initial stage of lateral inertial migration. x is the flow direction in both figures. (a) At φ b = 1% and t = 200H / U particle trains are formed parallel to the streamlines indicated by the lines oriented in the flow direction (from smaller towards larger x). Only the bottom half of the channel is shown. (b) At φ b = 5% and t = 70H / U , stratification of the two-phase flow is observed. The colors are related to the particle cross-stream position in the channel gap.

  FIG. 11. Instability of the interface between the particle-laden and the particle-free zones. (a) Front view of the whole channel. (b) Three clips (of thickness 2d) perpendicular to the streamwise direction. (c) Mixture velocity profile (solid line) compared against the plane Poiseuille profile (plus symbols).

FIG. 12 .

 12 FIG. 12. Profiles of the contributions to the second normal stress difference scaled by the wall shear stress at Re = 1400 and t = 70H / U . Thick solid line: v ′2w ′2 . Thin solid line: S sl 22 -S sl 33 . Dashed line: prediction of S sl 22 -S sl 33 from the single particle in unbounded pure shear flow considering local Re p and shear rate.

FIG. 13 .

 13 FIG. 13. (a) Contours of the streamwise velocity fluctuations scaled by U followed by the spanwise velocity gradient scaled by the average flow shear rate. (b) Contours of streamwise vorticity scaled by the average shear rate. All these contours correspond to averages in the streamwise direction at the same time t = 200H / U at Re = 1400. In (c), the secondary flows are superposed with the particle front view of panel (a).

TABLE I .

 I Parameters of simulated homogeneous suspension flows.

	Configuration	Re	Re p	Re sim	Re p,sim
	Couette	γ H 2 /ν	Re(a/H ) 2	2 × 10 -3 1000	2 × 10 -6 1
	Poiseuille	U H /ν	3Re(a/H ) 2	1 × 10 -3 500	3 × 10 -6 1.5

TABLE II .

 II Increase of the viscosity due to the combined effect of concentration and inertia.

	φ b (%)	Re	S sl 12 y µ γ sl y φ b
	1	10 -3	2.58
	1	500	3.23
	5	10 -3	2.84
	5	500	3.22

ACKNOWLEDGMENTS

The authors would like to acknowledge the anonymous referee 1 for his/her valuable comments that improved the introduction section. We thank A. Pedrono for technical support, and P. Floquet for solving the differential equation of the particle agitation. This work was granted access to the HPC resources of CALMIP and GENCI under the Allocation Nos. 2013-P1002 and x20132a6942, respectively.