Stany Gallier
email: stany.gallier@herakles.com

Elisabeth Lemaire

Laurent Lobry

François Peters

Franc ¸ois Peters

A fictitious domain approach for the simulation of dense suspensions

Keywords: Suspensions, Fictitious domain, Lubrication, Discrete Element Method

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

Introduction

Suspensions of solid particles embedded in a liquid are a class of two-phase flows that are ubiquitous in industry as well as in biological or natural flows. Fresh concrete or uncured solid rocket fuel are two examples of industrial concentrated suspensions for which the highest particle volume fraction is desired while keeping correct rheological properties and flowing behavior. Such dense suspensions do exhibit an intricate physics which is hitherto far beyond complete understanding. This complexity mainly stems from the wide variety of fluid-particle or particle-particle interactions [START_REF] Coussot | Rheophysical classification of concentrated suspensions and granular pastes[END_REF] so that even the idealized case of hydrodynamically-interacting smooth monodisperse spherical particles in a Newtonian fluid can involve strong non-Newtonian effects such as yield stress, shear-thickening, particle migration or anisotropic microstructures (see [START_REF] Stickel | Fluid mechanics and rheology of dense suspensions[END_REF][START_REF] Morris | A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow[END_REF] for a review).

The development of numerical simulation can help shed light on the complex physics of suspensions. However, due to the importance of flow-particle interactions, only microscale methods -wherein the flow is fully resolved around each particle -are relevant. In contrast, macroscale methods are less computationally demanding but require a local averaging -because computational cells are much larger than particles -which loses the essential details of the flow. For low Reynolds number suspensions, well-suited techniques include Stokesian Dynamics (SD) [START_REF] Bossis | Dynamic simulation of sheared suspensions. I. General method[END_REF][START_REF] Brady | The rheology of concentrated suspensions of spheres in simple shear flow by numerical simulation[END_REF][START_REF] Durlofsky | Dynamic simulation of hydrodynamically interacting particles[END_REF][START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF] and Force-Coupling Method (FCM) [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF][START_REF] Climent | Numerical simulations of random suspensions at finite Reynolds numbers[END_REF][START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF]. Although different, both methods rely on a truncated multipole expansion solution of the Stokes equations. They undoubtedly have numerous upsides and appear to be instrumental in providing among the most influential results in the field of suspension physics. Nevertheless, they are highly specialized and may not tackle any kind of flows. Because they rest on special solutions of the Stokes equations (multipole expansion), they are inherently restricted to non-or weakly-inertial flows. By and large, multipole expansions are mostly available for spheres, thereby preventing from simulating arbitrary-shaped particles. Furthermore, it is not possible to consider non-Newtonian fluids with these techniques. Therefore, direct numerical simulation (DNS) has emerged as an attractive alternative. In this approach, governing equations (Navier-Stokes or Stokes equations) are solved without any further assumptions other than numerical approximations. This makes possible to solve particulate flows with arbitrary particle shape, flow Reynolds number or fluid rheological constitutive equation. The very first class of DNS methods dedicated to particulate flows followed a boundary-fitted approach wherein only the domain occupied by the fluid is meshed [START_REF] Hu | Direct simulation of fluid particle motions[END_REF][START_REF] Feng | Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid Part 1. Sedimentation[END_REF]. As particles move, the constantly evolving fluid domain poses complex remeshing issues. For sheared concentrated suspensions, in which particle separation can be vanishingly small, remeshing becomes extremely involved and makes this approach impractical for more than a few particles. As an exception, let us mention the works of Johnson and Tezduyar [START_REF] Johnson | 3d simulation of fluid-particle interactions with the number of particles reaching 100[END_REF] who performed impressive 3D simulations with up to 100 particles with such body-fitted techniques. In contrast, non-boundary-fitted methods are much more suited for the simulation of suspensions with a large number of particles. The whole domain is mapped onto an Eulerian fixed grid and particles are embedded in this regular non-moving grid. Non-boundary-fitted methods for particulate flows include different techniques, such as immersed boundary methods [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Kim | Immersed boundary method for flow around an arbitrarily moving body[END_REF][START_REF] Feng | Proteus: a direct forcing method in the simulations of particulate flows[END_REF], fictitious domain methods [START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF][START_REF] Patankar | A new formulation of the distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF] or lattice Boltzmann methods [START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation[END_REF][START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2. Numerical results[END_REF][START_REF] Ladd | Lattice-Boltzmann simulations of particle-fluid suspensions[END_REF]. Fictitious domain methods have met considerable interest and have been widely used so far [START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF][START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF][START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF][START_REF] Veeramani | A fictitious domain formulation for flows with rigid particles: A non-Lagrange multiplier version[END_REF][START_REF] Blasco | A Fictitious Domain, parallel numerical method for rigid particulate flows[END_REF]. They have been successfully applied to non-Newtonian flows [START_REF] Patankar | Rheology of a suspension of particles in viscoelastic fluids[END_REF][START_REF] Hwang | Toward the computational rheometry of filled polymeric fluids[END_REF][START_REF] Yu | Viscoelastic mobility problem of a system of particles[END_REF][START_REF] Hwang | Direct simulations of particle suspensions in a viscoelastic fluid in sliding bi-periodic frames[END_REF], heat transfer in suspensions [START_REF] Yu | A fictitious domain method for particulate flows with heat transfer[END_REF] or non-spherical particles such as ellipsoids [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF], cubes [START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF], polygons [START_REF] Wachs | A DEM-DLM/FD method for direct numerical simulation of particulate flows: Sedimentation of polygonal isometric particles in a Newtonian fluid with collisions[END_REF], plates [START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF] or fibers [START_REF] Beaume | Modélisation et simulation numérique directe de l'écoulement d'un fluide complexe[END_REF]. Fictitious domain methods are generally body-force-based methods since solid particles are modeled via a body-force (or momentum forcing) introduced in the momentum equation to enforce a rigid body motion. There exist various methods to handle this body-force and a review is provided in Yu and Shao [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF]. In the original fictitious domain approach by Glowinski et al. [START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF][START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF], the body-force is introduced as a Lagrange multiplier in a weak formulation and solved in an iterative way. However, the method is computationally expansive as a saddle-point problem must be solved.

Patankar's works [START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF][START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF] were the first to develop a non-Lagrange multiplier version of fictitious domain in order to improve computational performance. More recent studies [START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF][START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF][START_REF] Veeramani | A fictitious domain formulation for flows with rigid particles: A non-Lagrange multiplier version[END_REF] have still continued in this way of non-Lagrange multiplier methods by working out an explicit equation on the body-force. These previously cited studies share many similarities and only slightly differ. Surprisingly, despite the availability of these DNS methods, quantitative simulations of dense suspensions are heretofore a privilege of SD or FCM methods. Only scarce studies address the DNS of suspension rheology [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF][START_REF] Patankar | Rheology of a suspension of particles in viscoelastic fluids[END_REF][START_REF] Hwang | Toward the computational rheometry of filled polymeric fluids[END_REF][START_REF] Lefebvre | Apparent viscosity of a mixture of a Newtonian fluid and interacting particles[END_REF] but they remain mostly qualitative and far from actual configurations as they consider 2D particles (cylinders) with low volume fractions. As far as the authors know, there have been no attempts at simulating the rheological behavior of dense suspensions using DNS in the way it is done by SD or FCM. With an eye to demonstrating that DNS is well-adapted for detailed simulations of concentrated suspensions, this study intends to develop a non-Lagrange multiplier version [START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF][START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF] of the fictitious domain approach since it alleviates the computational burden of the original version. However, those methods define body-force values at some Lagrangian nodes attached to the particles moving in the Eulerian background grid. This consequently involves several interpolation steps between Eulerian and Lagrangian nodes, which are known to induce numerical instabilities [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF] and impose a careful choice of the interpolation kernel [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Peskin | The immersed boundary method[END_REF]. It is also attested that the arrangement itself of the Lagrangian particle nodes may alter results. For spheres, the best results are obtained with nodes adequately arranged in concentric shells and with the nodes closest to the particle surface slightly retracted from it [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF][START_REF] Yu | Numerical simulation of particle sedimentation in shear-thinning fluids with a fictitious domain method[END_REF]. A generalization to arbitrary-shaped particles is believed to be tedious. Therefore, our approach intends to keep the simplicity of a fully Eulerian method. Bodyforce quantities are defined at the Eulerian grid points, like the other fluid variables, and an Eulerian advection step is inserted in the method, so that the body-force remains attached to the particle as it moves along.

A peculiar feature of concentrated suspensions is that the average separation distance a⇠ between particles, where a denotes the particle radius, becomes extremely small. The so-called lubrication forces arise between particles in near-contact because of the draining of interstitial fluid in the gap. The magnitude of these forces rises dramatically as particles approach each other and is singular in the limit of touching particles as the normal and tangential forces diverge as ⇠ -1 and ln ⇠, respectively. Consequently, the rheology of suspensions is markedly modified by lubrication forces and qualitatively important effects occur at small separations, typically down to ⇠ ⇠ 10 -2 . For accurate simulations of low Reynolds suspensions, a numerical method should be able to capture those short-range lubrication forces. The required grid spacing should however be smaller than at least 10 -3 a-10 -4 a to resolve lubrication, making long-term simulations of many-particle systems unfeasible. Typical grid spacings are about 10 -1 a which means that -albeit implicitly accounted for by DNS -lubrication forces are less and less accurately described as particles come to near contact. A first simple approach consists in adding the theoretical lubrication force F th lub (⇠), known for two spheres in near-contact [40], to the computed hydrodynamic force. To put it more precisely, this rather takes the form F th lub (⇠) -F th lub (⇠ cut ) where ⇠ cut is a cut-off separation below which this lubrication correction is activated. This approach is widely employed whatever the numerical technique chosen : fictitious domain [START_REF] Yu | A fictitious domain method for particulate flows with heat transfer[END_REF], dissipative particle dynamics [START_REF] Martys | Study of a dissipative particle dynamics based approach for modeling suspensions[END_REF], lattice Boltzmann method [START_REF] Nguyen | Lubrication corrections for lattice-Boltzmann simulations of particle suspensions[END_REF] or boundary element method [START_REF] Nasseri | Lubrication approximation in completed double layer boundary element method[END_REF]. Yet, such simple approach is not rigorous because -as previously stressed -a part of the lubrication is already partially included in the numerical solution. Just adding the theoretical force consequently results in double-counting the resolved part of the lubrication force. A more rigorous method is proposed in SD [START_REF] Durlofsky | Dynamic simulation of hydrodynamically interacting particles[END_REF][START_REF] Brady | Stokesian dynamics[END_REF]. The analytical two-sphere resistance interactions are added to the long-range hydrodynamic resistance matrix. However, in order not to count twice the far-field part of the exact two-sphere resistance interactions, the authors subtract off this contribution, which is found by inverting a two-sphere mobility matrix to the same level of approximation as in the long-range hydrodynamic problem. This method is however not well suited for DNS because long-range hydrodynamic resistance or mobility matrices are not explicitly computed. In the frame of FCM, Maxey et al. [START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF][START_REF] Dance | Incorporation of lubrication effects into the force-coupling method for particulate two-phase flow[END_REF] developed an approach inspired by the method used in SD.

The main idea is to compute, without any lubrication correction, the interactions between two spheres for different configurations and separation distances ⇠, in order to estimate the numerically resolved part of lubrication. The basic ideas of this method can be retained and readily transposed to a DNS approach, which will be done in the frame of this study. Note that there exist other alternatives such as using X-FEM (Extended Finite Element Method) to enrich the solution with lubrication theory in order to capture the subgrid-scale lubrication layer flow [START_REF] Wagner | Particulate flow simulations using lubrication theory solution enrichment[END_REF]. This is an earlier example of a combination of DNS and lubrication (while avoiding double-counting) but was only demonstrated in 2D and for small numbers of particles (⇠10).

A last essential ingredient for accurate simulations of suspension physics is the modeling of collisions, or contacts, between particles. Contacts inevitably occur in inertial flows but also in non-inertial flows despite lubrication because of many-body interactions or particle roughness. Additionally, from a numerical standpoint, contact may also take place in dense suspensions due to finite time steps. In order to prevent from unphysical particle overlaps, most authors include an additional ad hoc short-range repulsion force. The functional form of this force is usually dictated by numerical convenience rather than physics. Nonetheless, the important role of contact forces in dense suspensions has been highlighted in recent experimental works that unequivocally show that moderately-sheared highly-concentrated suspensions behave like dry granular materials (e.g., see [START_REF] Huang | Flow of wet granular materials[END_REF]). As dry granular materials have long been known to be well modeled by a Discrete Element Method (DEM) (see [START_REF] Radjaï | Discrete element modeling of granular materials[END_REF][START_REF] Pöschel | Computational granular dynamics: models and algorithms[END_REF] for a review), couplings between DEM and DNS solvers have recently emerged [START_REF] Wachs | A DEM-DLM/FD method for direct numerical simulation of particulate flows: Sedimentation of polygonal isometric particles in a Newtonian fluid with collisions[END_REF][START_REF] Mansouri | A 3D DEM-LBM approach for the assessment of the quick condition for sands[END_REF]. Yet, those works do not consider lubrication and a peculiar feature of our present model is that both lubrication and contact using DEM are considered.

The goal of this paper is to present an efficient method suited for the simulation of concentrated suspensions, thereby including a detailed modeling of lubrication and contact. Section 2 describes the fictitious domain method employed. It is close to the works of Yu and Shao [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF] and Apte et al. [START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF] although a significant difference is that our method is fully Eulerian and does not require Lagrangian collocation nodes. In Section 3, a robust and effective lubrication method is presented. It retains the basic ideas of the recent method proposed in FCM [START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF] and is adapted to the present DNS approach. In Section 4, the DEM method is delineated. Section 5 presents different comprehensive validations especially focusing on spheres in linear shear flows. Finally, conclusions are given in Section 6.

Fictitious domain method

Governing equations

Let D be the computational domain including the fluid domain D f and the solid particle domain D p , so that D = D f ⊕D p . The particle domain is assumed to be made up of N particles P i : D p = S i=N i=1 P i . In the frame of the present work, particles are assumed to be rigid and homogeneous with density ⇢ p . The boundary between fluid and particle domains is noted @P. The idea of fictitious domain is to extend the Navier-Stokes equations in the whole domain D including particles. They are thus supposed to be filled with the fluid, which is supposed here to be Newtonian and incompressible with density ⇢ f and viscosity µ. Incompressibility for velocity u in the whole domain

D yields r•u = 0 (1)
whereas the momentum equation in the whole domain D reads:

⇢ f ( @u @t + u •ru) = r•Σ + ⇢ f λ (2)
where the term ⇢ f λ can be seen as a body-force that enforces the rigid body motion inside particles. This force is non-zero in the particle domain D p and zero in the fluid domain D f . The fluid stress tensor Σ reads for a Newtonian fluid :

Σ = -pI + 2µE (3)
where p is the fluid pressure and E the rate-of-strain tensor E = 1 2 (ru + ru T ). Note that in Eq. (2), gravity is not included since the whole domain D is filled with the same fluid of density ⇢ f . Buoyancy forces will be included in a subsequent step. Particles are supposed to be rigid, so that the velocity field inside particles is given by a rigid body motion

u = U + Ω ⇥ r (4)
where U and Ω are the particle translational and rotational velocities, respectively, and r denotes the position vector with respect to particle center of gravity r = xx G .

The forthcoming developments intend to construct some explicit expressions for particle velocities and forcing term λ and are close to Yu and Shao [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF]. Particle motion is given by Newton's equations and reads

M dU dt = F hyd + (1 - ⇢ f ⇢ p )Mg + F int (5) J. dΩ dt + Ω ⇥ (J.Ω) = T hyd + T int ( 6
)
where g is the gravity acceleration and F int and T int are the force and torque due to interactions on particles. In present work, they will describe the contact between particles but they also can represent any other external interactions such as electrostatic forces for instance. F hyd and T hyd are respectively the hydrodynamic force and torque exerted by the fluid on the particle :

F hyd = Z @P n.Σ dS (7)
T hyd = Z @P r ⇥ (n.Σ)dS ( 8
)
where n is the outward-pointing unit normal vector on the boundary @P. In Eq. ( 5) and Eq. ( 6), M and J stand for the particle mass and inertia tensor :

M = Z P ⇢ dx (9) J = Z P ⇢[(r.r)I -r ⌦ r]dx (10)
These integrals are easily handled by introducing a generalized density ⇢(x) given as

⇢(x) = ⇢ f (1 -I P (x)) + ⇢ p I P (x) (11)
where I P (x) is the particle indicator function:

I P (x) = 8 > > > > > < > > > > > : 1 if x 2P 0 if x < P (12)
Integrating momentum equation Eq. ( 2) as well as r ⇥ Eq. (2) on a particle P and using the rigidity constraint Eq. ( 4)

and the definition of hydrodynamic force Eq. ( 7) and torque Eq. ( 8) yields

F hyd = -⇢ f Z P λ dx + ⇢ f ⇢ p M dU dt (13)
T hyd = -⇢ f Z P r ⇥ λ dx + ⇢ f ⇢ p ⇥ J. dΩ dt + Ω ⇥ (J.Ω) ⇤ (14)
Using these expressions in Eq. ( 5) and Eq. ( 6) provides

⇢ f Z P λ dx = (⇢ f -⇢ p ) ⇢ p M[ dU dt -g] + F int (15) ⇢ f Z P r ⇥ λ dx = (⇢ f -⇢ p ) ⇢ p ⇥ J. dΩ dt + Ω ⇥ (J.Ω) ⇤ + T int (16)
Eventually, the fluid-particle problem is governed by the following equations

@u @t + u •ru = - rp ⇢ f + ⌫r 2 u + λ in D (17) r•u = 0 in D (18)
u = U + Ω ⇥ r in D p (19) ⇢ f Z P λ dx = (⇢ f -⇢ p ) ⇢ p M[ dU dt -g] + F int (20) ⇢ f Z P r ⇥ λ dx = (⇢ f -⇢ p ) ⇢ p ⇥ J. dΩ dt + Ω ⇥ (J.Ω) ⇤ + T int ( 21
)
where ⌫ = µ/⇢ f is the fluid kinematic viscosity. At this step, the governing equations are essentially similar to the ones proposed by Yu and Shao [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF] and Apte et al. [START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF] and are a strong-form counterpart of the classical fictitious domain method introduced by Glowinski et al. [START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF][START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF]

Numerical procedure

A fractional-step time scheme is used to decouple the previous system Eq. ( 17)-( 21) into two subproblems. The first subproblem Eq. ( 17)-( 18) is a standard Navier-Stokes problem in D whereas Eq. ( 19)-( 21) form a particle subproblem which mainly consists in enforcing the rigid body motion inside particles. The numerical procedure is presented for the general case of 3D flows with arbitrary-shaped particles although only spheres will be considered in this paper.

Fluid subproblem

The standard incompressible Navier-Stokes equations Eq. ( 17)-( 18) are classically solved using finite differences on half-staggered Cartesian grids. The whole domain D is uniformly meshed with a constant grid spacing ∆. The components of the forcing term λ are half-staggered in the same way as the velocity components. As we primarily focus on low Reynolds number flows, the time integration scheme is implicit (Crank-Nicolson) for the diffusive terms and explicit for the convective terms using a second-order Adams-Bashforth scheme :

u n+1 -u n ∆t + [ 3 2 (u •ru) n - 1 2 (u •ru) n-1 ] + rp n+1 ⇢ f = ⌫ 2 r 2 (u n+1 + u n ) + λ n (22) r•u n+1 = 0 ( 23
)
This formulation is similar to the non-Lagrange multiplier methods already mentioned [START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF][START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF][START_REF] Veeramani | A fictitious domain formulation for flows with rigid particles: A non-Lagrange multiplier version[END_REF] and differs from the original fictitious domain method [START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF][START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF] in the sense that the rigid body motion is here imposed only in an approximate way because the forcing at the previous iteration λ n is used. The system Eq. ( 22)-( 23) is classically solved using a projection method [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF]. In such approach, the momentum equation Eq. ( 22) is first solved while discarding the pressure term (i.e. rp n+1 =0). The obtained velocity field, noted u ⇤ , is consequently not divergence-free (r•u ⇤ , 0). The subsequent projection step is obtained on noting that

u ⇤ = u n+1 + ∆t ⇢ f rp n+1 ( 24
)
which is an Helmoltz-Hodge decomposition of velocity u ⇤ . Taking the divergence of Eq. ( 24) and imposing a divergence-free condition for velocity u n+1 gives a Poisson equation on the pressure

r 2 p n+1 = ⇢ f ∆t r•u ⇤ (25)
The resulting pressure field p n+1 is then used in Eq. ( 24) so as to obtain the desired velocity u n+1 .

The space discretization of gradient or laplacian operators is achieved using classical second-order centered schemes. The four resulting linear systems (three for each velocity component in Eq. ( 22), one for pressure in Eq. ( 25))

are solved by a geometric multigrid with V-cycles and two Gauss-Seidel iterations for pre-and post-smoothing.

Particle subproblem

The particle subproblem is three-fold : first, constrain a rigid body motion inside particles ; second, compute particle velocities ; and third, update the forcing term λ. As previously stressed, the fluid velocity u n+1 does not exactly comply with a rigid body motion inside particles because the used forcing term λ n is taken at the previous time step. At this stage, an explicit direct forcing of the rigid body motion is thus required :

ũn+1 = U n+1 + Ω n+1 ⇥ r in D p ( 26
)
This new velocity ũn+1 now fulfills an exact rigid body motion inside particles. It can be noted that the difference between ũn+1 and u n+1 is a measure of the defect in constraining the exact rigid body condition and can be used to update λ. To this end, the discretized Navier-Stokes equation Eq. ( 22) can be rewritten in the following compact form

u n+1 -u n ∆t = r•F n+1 + λ n ( 27
)
The new velocity ũn+1 describes an exact rigid body motion, so that

ũn+1 -u n ∆t = r• F n+1 + λ n+1 (28)
By difference, we obtain

λ n+1 = λ n + ũn+1 -u n+1 ∆t -r•( F n+1 -F n+1 ) (29)
Rigorously speaking, fluxes F n+1 and F n+1 are not exactly equal because velocities have been slightly modified due to the correction step Eq. ( 26). Indeed, F n+1 is not explicitly known and depends on λ n+1 so that Eq. ( 29) is an implicit equation of limited practical interest. In order to have an explicit equation on λ n+1 , the following approximation is used

F n+1 ⇡ F n+1 (30)
which is a reasonable assumption for moderate time steps. This eventually gives a simple explicit equation for the forcing term :

λ n+1 = λ n + ũn+1 -u n+1 ∆t (31)
The last step consists in computing the particle translational and rotational velocities U n+1 and Ω n+1 that are actually required to enforce the rigid body motion Eq. ( 26). The discretized counterparts of Eq. ( 20) and Eq. ( 21) read

⇢ f Z P λ n+1 dx = (⇢ f -⇢ p ) ⇢ p M[ U n+1 -U n ∆t -g] + F int,n (32)
⇢ f Z P r ⇥ λ n+1 dx = (⇢ f -⇢ p ) ⇢ p [ J.(Ω n+1 -Ω n ) ∆t + Ω n ⇥ (J.Ω n )] + T int,n (33)
Substituting Eq. ( 26) in Eq. [START_REF] Hwang | Direct simulations of particle suspensions in a viscoelastic fluid in sliding bi-periodic frames[END_REF] gives

λ n+1 = λ n + U n+1 + Ω n+1 ⇥ r -u n+1 ∆t (34)
Finally, Eq. ( 34) is integrated on the particle volume and inserted in Eq. ( 32) and, similarly, r ⇥ Eq. ( 34) is integrated and inserted in Eq. ( 33). This leads to the explicit equations

MU n+1 = (⇢ p -⇢ f ) ⇢ p M[U n + g∆t] + ⇢ f Z P (u n+1 -λ n ∆t)dx + F int,n ∆t (35)
J.Ω n+1 = (⇢ p -⇢ f ) ⇢ p [J.Ω n -Ω n ⇥ (J.Ω n )∆t] + ⇢ f Z P r ⇥ (u n+1 -λ n ∆t)dx + T int,n ∆t (36)
Let us remark that in the general case of non-spherical particles, the inertia tensor J is not diagonal which requires an additional 3⇥3 system inversion in Eq. ( 36) so as to compute the rotational velocities. These expressions for velocities are essentially similar to the ones obtained by other non-Lagrange multiplier techniques [START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF][START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF][START_REF] Veeramani | A fictitious domain formulation for flows with rigid particles: A non-Lagrange multiplier version[END_REF].

Particle tracking

Non-boundary-fitted Eulerian methods are attractive since a fixed mesh avoids any remeshing difficulties. However, the interface between fluid and particles must be defined properly. A large amount of work is devoted to this interface tracking problem, such as VOF (Volume-Of-Fluid) methods [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] or level-set methods [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Osher | Level Set Methods: An Overview and Some Recent Results[END_REF] (see also [START_REF] Lakehal | Interface tracking towards the direct simulation of heat and mass transfer in multiphase flows[END_REF] for a review.) Still, this problem is not acute here because only rigid particles are considered and their shape consequently remains unchanged. In present case, the interface tracking is needed so as to delineate the fluid region D f and the particle region D p . The method chosen here rests on the well-known level-set approach [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Osher | Level Set Methods: An Overview and Some Recent Results[END_REF] which considers a level-set function (x) negative inside particles and positive outside and with the interface defined by the hypersurface (x)=0. A widely employed level-set function is the signed distance : for the simple case of a sphere of radius a centered at x G , the level-set reads

(x) = kx -x G k-a (37)
The previously introduced particle indicator function I P (x) can be directly computed from the level-set function :

I P (x) = 1 -H( (x)) ( 38
)
where H is the Heaviside function

H(s) = 8 > > > > > < > > > > > : 0 if s  0 1 if s > 0 (39)
This Heaviside function is smeared over a small region ✏ (proportional to grid spacing) around the boundary giving a smooth variation that will enhance numerical accuracy and robustness, such as

H ✏ (s) = 1 2 (1 + tanh(b s ∆ )) ( 40
)
where ∆ is the grid spacing and b is a free parameter controlling the smearing on a typical size ✏ ⇠O (∆/b). In the present study, a value b=5 is chosen which gives a correct smoothing in typically less than a grid spacing. The corresponding smooth indicator function I ✏ P is then used to compute the particle mass, volume, inertia tensor as well as any integral of a quantity φ on the particle volume P

Z P φ dx = Z D(P) I ✏ P (x).φ dx ( 41
)
where the integration domain D(P) is a subset of D encompassing the particle P (bounding box).

Particles are assumed to be rigid so that their position can be tracked by solving the following transport equation

dX dt = U ( 42
)
where X denotes the vector of particle centers. For non-spherical particles, X also contains orientations (Euler angles or quaternions) and U additionally includes rotational velocities. Note however that the present study will only focus on spheres. Using a second-order Adams-Bashforth scheme, the previous equation reads

X n+1 = X n + 3 2 ∆tU n+1 - 1 2 ∆tU n (43)
A peculiar feature of our approach is to be fully Eulerian : the forcing λ is not defined at some particle collocation nodes but rather at the Eulerian grid points like any other fluid variables. Consequently, an Eulerian transport step is inserted in the method so that the forcing remains attached to the moving particle. Hence, the following advection equation is solved for each particle P @λ P @t + U P .rλ P = 0

where subscript P has here been added to recall that this equation needs to be solved for each particle. Practically speaking, this equation is solved only in the bounding box D(P) to avoid unnecessary computational effort. Equation ( 44) is a classical advection equation for which a large number of methods is available. However, the forcing term λ P has a very sharp variation across the particle boundary (it is zero in the fluid domain) and even high-order schemes, such as a fifth-order WENO (Weighted Essentially Non-Oscillatory) scheme [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF] or a fourth-order WKL (Warming-Kutler-Lomax) scheme [START_REF] Huang | Semi-Lagrangian advection schemes and Eulerian WKL algorithms[END_REF], failed to correctly advect the forcing on moderately coarse grids without any appreciable dissipation and dispersion. The best performance is found to be obtained by semi-Lagrangian schemes [START_REF] Huang | Semi-Lagrangian advection schemes and Eulerian WKL algorithms[END_REF][START_REF] Staniforth | Semi-lagrangian integration schemes for atmospheric models: A review[END_REF][START_REF] Bermejo | The conversion of semi-Lagrangian advection schemes to quasi-monotone schemes[END_REF] that are furthermore ideally suited for Cartesian grids. The advection step Eq. ( 44) is written in Lagrangian form Dλ P /Dt = 0 which means that quantity λ P is conserved on the characteristic path dx/dt = U P . The semi-Lagrangian scheme then simply reads

λ P (x, t +∆t) = λ P (x -U P ∆t, t) (45)
which just requires an interpolation of the forcing at position x-U P ∆t. The overall accuracy of the scheme is governed by the accuracy of this interpolation. Although (tri-)linear interpolation is simple, it induces much dissipation. A second-order Lagrange polynomial interpolation is preferred as a good trade-off between accuracy and computational demand. Note that this semi-Lagrangian scheme is reminiscent of high-order Lagrangian/Eulerian interpolation steps in usual fictitious domain methods as discussed in the introduction. However, present interpolation step is performed only once per time step, just for advecting the momentum forcing, whereas three Lagrangian/Eulerian interpolation steps are generally required per iteration in the aforementioned methods (e.g., see [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF]).

Numerical algorithm

The following steps are performed (where n is the current iteration) :

1. Compute level-set function

For each particle P, define a cubic bounding box D(P) encompassing the particle and compute the level-set function n P (x) Eq. ( 37) and associated smeared indicator function I ✏,n P (x) in the bounding box 2. Advect momentum forcing For each particle P, advect the forcing term from λ n to λ n ⇤ by solving the transport equation Eq. ( 44) in D(P) for one time step ∆t

Solve fluid subproblem

Compute new fluid velocities u n+1 and pressure p n+1 (Eq. ( 22) and Eq. ( 23)) using λ n ⇤ as the momentum forcing.

Update particle velocities

Update particle translational U n+1 and rotational velocities Ω n+1 using Eq. ( 35) and Eq. ( 36) with λ n ⇤ as the momentum forcing and with interaction forces/torques

F int,n =F int (X n , U n ) and T int,n =T int (X n , U n ).

Enforce rigid body motion

Explicitly enforce the rigid body motion ũn+1 = U n+1 + Ω n+1 ⇥ r (Eq. ( 26)) inside each particle

Update momentum forcing

Correct the momentum forcing using Eq. ( 31) :

λ n+1 = λ n ⇤ + ( ũn+1 -u n+1
)/∆t 7. Update particle position Update particle position using Eq. ( 43) :

X n+1 = X n + 3 2 ∆tU n+1 -1 2 ∆tU n 8. Repeat
Set n to n+1 and go to step 1 A last remark concerns the steady Stokes equations that will be extensively used in the forthcoming validations due to the availability of analytic solutions in this regime. Rather than considering extremely large fluid viscosity to obtain vanishing Reynolds numbers -which would inevitably result in a significant time step reduction despite the implicit time-stepping -it is preferred to discard the convective term (u •ru) in the Navier-Stokes equations. Furthermore, we loop over step 3 to 6 in pseudo-time until convergence before moving to the next steps. This actually consists in reaching a steady state Stokes flow before updating particle positions. This steady state is typically achieved within 2-3 subiterations.

Lubrication model

As stressed in the introduction, lubrication forces play a major role in the suspension rheology and need to be modeled properly. Since they are short-range in nature, they can usually not be fully resolved with the typical grids used and consequently require an ad hoc model. The method proposed is partly inspired from Stokesian Dynamics [START_REF] Durlofsky | Dynamic simulation of hydrodynamically interacting particles[END_REF][START_REF] Brady | Stokesian dynamics[END_REF] but has been here adapted to comply with a DNS framework.

Theoretical framework

Let us consider a system of N spherical particles suspended in a linear Stokes flow and let U be the 6N vector of translational/rotational velocities U = (U, Ω) T and F = (F, T) T the 6N vector of hydrodynamic forces/torques exerted by the fluid on the particles. Due to the linearity of the Stokes equations, there are linear relations between the forces/torques and the flow parameters and the velocities of particles can be written in resistance form as [40]

F = R FU • (U 1 -U) + R FE : E 1 ( 46
)
where U 1 and E 1 are the unperturbed velocities and the rate-of-strain tensor due to the prescribed flow field, respectively. The 6N ⇥ 6N resistance matrices R FU et R FE are not known in general except for the case of two spheres. The key idea is to split the hydrodynamic interactions into long-range interactions -explicitly resolved by the numerical model -and a short-range lubrication contribution which can not be resolved since it is subgrid :

R ⇡ R + R sub ( 47
)
where the tilde denotes the resolved part and the superscript sub refers to the subgrid unresolved part of the interaction.

Expanding Eq. ( 46) together with Eq. ( 47) yields

U 1 -U = R-1 FU • [F -RFE : E 1 -R sub FU • (U 1 -U) -R sub FE : E 1 ] (48)
The resolved velocity Ũ which is given by the numerical simulation formally reads

U 1 -Ũ = R-1 FU • [F -RFE : E 1 ] (49)
Combining Eq. ( 48) and Eq. ( 49) gives the exact velocity U in terms of the numerically resolved velocity Ũ

U = Ũ + R-1 FU • [R sub FU • (U 1 -U) + R sub FE : E 1 ] (50)
This basically means that the lubrication-corrected velocity U can be obtained from the resolved velocity Ũ by adding the external force F lub in the numerical procedure with

F lub = R sub FU • (U 1 -U) + R sub FE : E 1 (51)
This force represents the portion of hydrodynamic interactions not resolved by the numerical approach. The subgrid resistance matrices R sub FE and R sub FU will be specified subsequently.

Numerical procedure

Newton's equations for particles read

M • dU dt = F hyd + F lub + F int ( 52
)
where F hyd represents the resolved hydrodynamic interactions (including gravity), F lub the additional lubrication force Eq. ( 51), F int the contact interactions, and M the generalized mass/inertia matrix. Lubrication forces strongly depend on particle positions and velocities and are only weakly coupled with long-range hydrodynamics, so that Eq. ( 52) can be split into

M • Ũ -U n ∆t = F hyd + F int (53) M • U n+1 - Ũ ∆t = F lub ( 54
)
The first step Eq. ( 53) corresponds to the numerical procedure to compute particle velocity and is a formal rewriting of Eqs. ( 35)-( 36) while the second step Eq. ( 54) intends to correct the obtained particle velocity from lubrication. Since lubrication forces are very short-range, Eq. ( 54) is stiff and must be solved implicitly using F lub (U n+1 ). Rearranging

Eq. ( 54) with Eq. ( 51) yields

[ M ∆t + R sub FU ] • U n+1 = M ∆t • Ũ + R sub FU • U 1 + R sub FE : E 1 (55)
From the known uncorrected velocity Ũ, the actual lubrication-corrected particle velocity U n+1 is obtained by solving the linear system Eq. ( 55). Since the resistance matrix R sub FU is symmetric and sparse, an iterative conjugate gradient method is used to solve the system until a prescribed relative tolerance is achieved (typically 10 -6 ). By and large, iterative techniques work well for matrices that are well-conditioned (i.e., with eigenvalues tightly clustered).

However, resistance matrices are generally ill-conditioned due to the singular behavior of the resistance functions near contact, which induces very poor convergence, if any. Therefore, a preconditioned conjugate gradient with an incomplete Cholesky factorization with zero fill-in IC(0) is used. The efficiency of the preconditioning is further improved by a matrix reordering using a Reverse Cuthill-McKee (RCM) algorithm. As already noticed by Sierou and Brady [START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF], the preconditioned iterative procedure works well when the resistance matrix is reasonably well-behaved, in practice when the inter-particle distance is typically no less than 10 -6 a. As a result, a threshold distance of 10 -6 a is prescribed when constructing the resistance matrix with an eye to avoiding poor convergence rates.

In the numerical algorithm already presented (see Sec. 2.4), this lubrication correction step is performed just after the computation of particle velocities (step 4) and is therefore inserted between step 4 and step 5.

Subgrid resistance matrix

The subgrid resistance matrices R sub FU and R sub FE do represent the portion of hydrodynamic interactions not resolved by the numerical procedure. Just like in SD, subgrid resistance matrices are estimated by subtracting the two-body resistance matrices R2B obtained numerically from the exact two-body resistance matrices R 2B,theo known from lubrication theory :

R sub = R 2B,theo - R2B (56)
For two spherical particles (1) and ( 2

u 1 (x 1 ) -U 1 u 1 (x 2 ) -U 2 ! 1 (x 1 ) -Ω 1 ! 1 (x 2 ) -Ω 2 E 1 E 1 3
The various resistance tensors can be written in terms of several scalar functions [40]

A ↵β ij = X A ↵β d i d j + Y A ↵β (δ ij -d i d j ) (59)
B ↵β ij = Y B ↵β ✏ i jk d k (60) C ↵β ij = X C ↵β d i d j + Y C ↵β (δ ij -d i d j ) (61)
G ↵β i jk = X G ↵β (d j d k - 1 3 δ jk )d i + Y G ↵β (d j δ ik + d k δ ij -2d i d j d k ) (62)
H ↵β i jk = Y H ↵β (✏ jil d l d k + ✏ kil d l d j ) (63)
in which d = r/krk and r is the separation vector between particle centers r = x 1x 2 . In the previous expressions, δ and ✏ stand for the Kronecker and Levi-Civita tensor, respectively. The different scalar functions X ↵β and Y ↵β in Eq. ( 59)-( 63) depend on radius a and non-dimensional distance ⇠ = krk/a -2, and their analytic forms are given in [40,[START_REF] Jeffrey | Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow[END_REF].

In order to calculate R2B , it is first assumed that this resistance tensor has the same functional form as the theoretical R 2B,theo hereinbefore given by Eq. ( 58)- [START_REF] Gallier | Approche par domaine fictif pour la simulation des suspensions[END_REF]. The associated resolved scalar functions X↵β and Ỹ↵β are then computed with our DNS solver -without any lubrication correction -for pairs of particles having various configurations (orientation and separation distance). The obtained resolved scalar functions are subsequently tabulated against the non-dimensional separation distance ⇠. This step is done once for all. From a numerical viewpoint however, the solution slightly differs depending where the gap is located on the grid, especially if a coarse grid is used. For example, for two particles separated by one grid length, the computed lubrication force is different depending on whether the gap is spanned by a single grid cell, compared with two half-cells. In order to limit those grid-induced effects, the tabulated scalar functions are averaged over several configurations (for a given gap ⇠) that will span different relative positions on the grid. The resulting average values minimize grid effects.

Finally, subgrid resistance matrices are readily computed using Eq. ( 56), i.e. subtracting the tabulated resolved resistance matrix R2B previously obtained from the exact resistance matrix R 2B,theo given by Eq. ( 58).

Some numerical tests performed show that when the separation distance a⇠ between particles is larger than the grid spacing ∆, the numerical model is by itself able to resolve all the hydrodynamic interactions (i.e., R2B = R 2B,theo ) and therefore does not require any lubrication model. In contrast, when particles are separated by less than one grid spacing, the numerical solver can no longer resolve the near-field interactions accurately and the present lubrication correction technique must be activated. With the typical grid resolution used in this study (∆ ⇡ a/5), it means that the cut-off distance ⇠ lub for activating the lubrication correction (the so-called lubrication barrier) is about 0.2, a value also noticed in other methods like lattice-Boltzmann techniques [START_REF] Nguyen | Lubrication corrections for lattice-Boltzmann simulations of particle suspensions[END_REF].

For a many-particle system, the resistance matrices R sub are computed in a pairwise additive fashion using the classical assumption of pairwise additivity of forces. This basically corresponds to constructing the grand resistance matrix by summing up all the pair resistance matrices. Note that at low volume fractions, there are many particles which do not have particles in their neighborhood. Only pairs of particles separated by a distance smaller than the lubrication barrier ⇠ lub are actually considered in the lubrication correction.

As a last remark, it is worthwhile to note that the resolved subgrid scalar functions X↵β and Ỹ↵β have been tabulated

for wall-particle interactions as well using a similar procedure. The theoretical wall-particle scalar functions are given in Yeo and Maxey [START_REF] Yeo | Dynamics of concentrated suspensions of non-colloidal particles in couette flow[END_REF] and some validations of our approach for wall-particle interactions may be found in [START_REF] Gallier | Approche par domaine fictif pour la simulation des suspensions[END_REF].

Stresslet correction

The effective stress tensor of a suspension is linked to the first moment of fluid stress on particle surface [64]

Σ p = Z P (Σ • n) ⌦ x dx (64)
The symmetric part of this tensor Σ p is referred to as the hydrodynamic stresslet S and represents a portion of the total particle stress. It plays a major role in the rheology of suspensions and must be corrected from lubrication forces as well. The method employed for the stresslet is identical to the one previously described for velocities. The deviatoric stresslet can be written in resistance form as

S = R SU • (U 1 -U) + R SE : E 1 (65)
and can be similarly decomposed into a resolved and a subgrid part

S = S + R sub SU • (U 1 -U) + R sub SE : E 1 ( 66
)
where S corresponds to the resolved stresslet computed by the code. This stresslet can be obtained directly from the momentum forcing λ and, neglecting inertia, is given as [START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF] S = -⇢ f

Z P (r ⌦ λ)dx (67)
At this step, the exact velocity U is known and can be used in Eq. ( 66). The subgrid resistance matrices R sub SU and R sub SE are obtained as described in the previous section. Their theoretical expression for a pair of particles can be found in [40]. Finally, a similar correction procedure is also applied to the trace of Σ p using the theoretical resistance functions from Jeffrey et al. [START_REF] Jeffrey | The pressure moments for two rigid spheres in low-reynolds-number flow[END_REF].

By and large, the global lubrication model presented in the frame of this study is general although it will be here applied only for equally-sized spheres. The method can be extended to spheres with different size since lubrication theory is known in this case [40]. The only difference is that the scalar functions also depend on the size ratio. The numerically resolved scalar functions, needed to construct subgrid matrices, can be tabulated in a similar way although computations are a bit more tedious as different size ratios must be considered. Yet, an extension to arbitrary-shaped particles is limited because of the current lack of a complete theoretical framework for lubrication in that case. Some first elements have been however proposed by Claeys and Brady [START_REF] Claeys | Lubrication singularities of the grand resistance tensor for two arbitrary particles[END_REF]. Similarly, it seems that no general lubrication theory is available for non-Newtonian fluids. Thus, the limitations of our approach concerning lubrication arise from the availability of a general lubrication theory.

Collision model

Although lubrication forces prevent collisions in the case of two smooth non-inertial spherical particles, multibody interactions as well as inertia or surface roughness are liable to involve contacts between particles. In the field of suspension physics, contact forces are well known to induce irreversibilities which markedly alter the microstructure of suspensions [START_REF] Stickel | Fluid mechanics and rheology of dense suspensions[END_REF][START_REF] Morris | A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow[END_REF]. Moreover, contact is generally promoted by particle asperities which means that a general collision model must also account for surface roughness. Instead of using ad hoc short-range repulsion force as is usually done, this work intends to consider a specific model describing the physics of collisions. A popular model used in granular physics is the so-called DEM (Discrete Element Method) which considers particles individually as discrete entities. This method actually encompasses different models (see [START_REF] Radjaï | Discrete element modeling of granular materials[END_REF][START_REF] Pöschel | Computational granular dynamics: models and algorithms[END_REF] for a review) and the most common is the molecular dynamics variant developed in the 70s [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. The consideration of DEM as a collision model in DNS simulations has recently emerged [START_REF] Wachs | A DEM-DLM/FD method for direct numerical simulation of particulate flows: Sedimentation of polygonal isometric particles in a Newtonian fluid with collisions[END_REF][START_REF] Mansouri | A 3D DEM-LBM approach for the assessment of the quick condition for sands[END_REF][START_REF] Feng | Combined three-dimensional lattice boltzmann method and discrete element method for modelling fluid-particle interactions with experimental assessment[END_REF]. As the collisional time scale is much smaller than the fluid time scale, both phenomena can be decoupled so that DEM can be easily implemented as a separate submodel in the fluid solver.

Still, this will not be the case here because of lubrication. Since lubrication and contact are governed by similar time scales, they need to remain coupled.

Review of the Discrete Element Method (DEM)

The DEM approach considers granular media as a collection of particles that will interact through collisions described by forces accounting for elastic deformation and friction. Let us consider a pair of spherical particles P i and P j (of radius a i and a j ) undergoing contact. The contact interaction force F int exerted by particle P j on P i is classically decomposed into its normal F int n and tangential F int t components :

F int = F int n + F int t ( 68
)
The contact is modeled by a Kelvin-Voigt behavior : the normal force is assumed to be the sum of an elastic restoring force, proportional to the overlap distance δ ij = kr ij k-a ia j where r ij = x ix j , and a dissipative component proportional to the relative normal velocity

F int n = [-k n δ ij -γ n dδ ij dt ] • n ij ( 69
)
where k n and γ n are the normal stiffness and damping coefficients. The normal unit vector is defined as

n ij = r ij /kr ij k.
Parameter k n can be expressed in terms of mechanical properties (Young modulus and Poisson coefficient) or, alternatively, chosen sufficiently large to describe a rigid solid. The damping coefficient γ n is generally estimated in order to match a given restitution coefficient e n (ratio between post-collisional and pre-collisional relative normal velocities),

which is an easier quantity to measure. Defining the damping parameter  = γ n /γ cr n where γ cr n = 2 p k n M is the critical damping (with M the particle mass), the normal restitution coefficient e n is

e n = exp -⇡ p 1 - 2 (70)
Similarly, the tangential force is given by

F int t = -k t Υ ij -γ t dΥ ij dt ( 71
)
where Υ ij is a tangential spring defined by integrating the slip velocity U s ij during the contact

Υ ij = Z t 0 U s ij dt ( 72
)
where the slip velocity is

U s ij = U i -U j -[(U i -U j ) • n ij ] • n ij -(a i Ω i + a j Ω j ) ⇥ n ij ( 73
)
Because the tangential plane may also vary with time, the obtained tangential force must be projected back on the current tangential plane after each time step as

F int t = F int t -(F int t • n ij ) • n ij ( 74
)
Using the classical Amontons-Coulomb law of friction, the actual tangential force is modified if it exceeds the friction limit µ d |F int n | and is then given by

F int t = µ d |F int n | F int t |F int t | (75)
where µ d is the dynamic friction coefficient. The tangential stiffness k t is k t = 2k n /7 [69, 70] and γ t is such that e t =e n . Note that low velocity impacts have generally a restitution coefficient e ⇡ 1 (e.g., see [START_REF] Labous | Measurements of collisional properties of spheres using high-speed video analysis[END_REF]). Moreover for low Stokes numbers, substantial energy is dissipated due to lubrication forces [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF] so that mechanical dissipation can be neglected. This legitimates the use of e n =e t =1 for low Reynolds suspensions. The reader is referred to Ref. [START_REF] Radjaï | Discrete element modeling of granular materials[END_REF] for a general discussion on the DEM parameters. Finally, the corresponding torque is

T int = a i n ij ⇥ F int ( 76
)
Contact forces also induce an additional particle stress tensor given, for each pair of spheres, by

Σ int p = F int ⌦ r ij (77)

Numerical implementation

For each time step, forces F int and torques T int are computed using Eq. ( 68) and Eq. ( 76) for all pairs of particle in contact and are integrated explicitly through Eqs. ( 35)- [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF]. One of the main computational effort actually lies in finding the pairs of particles in contact. An efficient collision detection technique is generally preferred ; in present case, a Verlet list method is used [START_REF] Verlet | Computer experiments on classical fluids. i. thermodynamical properties of lennard-jones molecules[END_REF].

Roughness is an important feature needed and can be readily accounted for in the model. Assuming sparse asperities of size h rug , contact now occurs between asperities and surface whenever kr ij k 6 a i + a j + h rug . An easy way to implement roughness thus consists in defining the overlap distance as

δ ij = kr ij k-a i -a j -h rug ( 78
)
In classical DEM, used for dry contacts (i.e., without fluid), the numerical time step ∆t must be a fraction of the collisional time scale ⌧ coll = (M/k n ) 1/2 . Since contact force is modeled by a mass-spring system, this quantity appears naturally as the time constant of the differential equation M δ + k n δ = f where f represents forces other than contact. (Note that mechanical dissipation during contact is generally neglected for stability issues.) For real material properties, k n can be quite large which eventually results in a drastic time step reduction. This is however different for a lubricated collision which is mainly controlled by viscous effects [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF]. The contact dynamics is now rather described by an equation of the form M δ + q δ/(δ + h rug ) + k n δ = f where q is linked to the lubrication force and depends on viscosity and particle size. Assuming δ ⌧ h rug , this equation becomes linear and has generally an overdamped behavior because of the predominance of lubrication damping. It consequently has two time scales : a very short time scale ⌧ 1 = Mh rug /q only due to lubrication and a long time scale

⌧ 2 = ⌧ 2 coll /⌧ 1 .
As lubrication is treated implicitly, the short time scale poses no stability problems. The second time scale ⌧ 2 is, for reasonable values of k n , much larger than the numerical time step and involves no time step limitations.

Validations

Single sphere in a creeping shear flow

The simple case of a single force-free torque-free sphere in a linear creeping shear flow is of interest due to the availability of analytic solutions. A neutrally-buoyant sphere of radius a is suspended at the center of a cubic Couette cell. The shear rate is γ and the cell size is L=20a with a grid spacing ∆=a/4.9, thereby leading to 97 3 grid points.

Periodic conditions are used in the spanwise and streamwise directions. A velocity of γL/2 and -γL/2 in the flow direction is prescribed at the top and bottom boundaries, respectively. The time step is ∆t=3.10 -3 γ-1 , corresponding to a CFL number based on the diffusional time scale CFL d =⌫∆t/∆ 2 of 50. Steady Stokes equations are solved ; as already highlighted in Sec. 2.4, this means that unsteady Stokes equations are solved to a steady state.

Velocity field

In a flow with prescribed rate-of-strain tensor E 1 and vorticity ! 1 , the theoretical fluid velocity induced by a rigid force-free torque-free spherical particle reads [START_REF] Batchelor | The hydrodynamic interaction of two small freely-moving spheres in a linear flow field[END_REF]

u i = u 1 i + E 1 ij x j + ✏ i jk ! 1 j x k -a 5 E 1 ik x k r 5 - 5a 3 2 (1 - a 2 r 2 )E 1 jk x i x j x k r 5 ( 79
)
given that the sphere is at the origin and r = kxk. Classically, subscript (1) refers to the direction of the flow while subscripts ( 2) and ( 3) denote the direction of the velocity gradient and vorticity, respectively. parameters chosen, the particle translational velocity U 1 is computed to be of the order of 10 -4 -the expected value is 0 -and the rotational velocity is found to be Ω 3 /γ=-0.501 quite close to the theoretical value (-0.5).

Computations are also conducted while discarding the forcing term λ (filled symbols in Fig. 1) so as to highlight its effect on the results. It can be clearly noted that the role of the momentum forcing is significant in order to obtain the expected velocity field. Discarding this term therefore results in poorly predicted hydrodynamic interactions although the rotational and translational particle velocities are nevertheless correct because the rigid body motion is still enforced by Eq. ( 26). This result can seem surprising since some published studies present satisfactory results without any momentum forcing [START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF][START_REF] Blasco | A Fictitious Domain, parallel numerical method for rigid particulate flows[END_REF][START_REF] Carlson | Rigid fluid: animating the interplay between rigid bodies and fluid[END_REF]. The difference is actually linked to the flow regime studied. The aforementioned studies consider inertial flows with Reynolds number Re ⇠O (10 1 -10 3 ) for which an explicit timestepping is used. The time step ∆t consequently remains small compared to the fluid time scale meaning that within one iteration, the velocity field inside particles has almost kept its previous rigid body motion : ũn+1 ⇡ u n+1 . From

Eq. ( 31), it is expected that λ remains small and superfluous. It is thus possible to predict correctly inertial flows without momentum forcing because the explicit direct forcing Eq. ( 26) is sufficient to impose the rigid body motion.

However, things are different for very low Reynolds number flows. Due to the necessary use of implicit time-stepping -to avoid vanishingly small time steps -the time step ∆t becomes much larger than the numerical diffusional time scale ⌧ d ⇠ ∆ 2 /⌫. This means that within one time step, the flow field is significantly modified by diffusion. It becomes thus necessary to introduce a body-force λ when solving the momentum equation so that the fluid "feels" the presence of the particle at this step.

Convergence study

A convergence study is conducted on the previous sheared single-sphere configuration. The quantity investigated here is the particle stresslet due to its importance in rheology. This stresslet can be computed directly from the momentum forcing λ as given by Eq. ( 67). Note that an accurate computation of the stresslet is inherently difficult because it only depends on λ which has very sharp variations across particle boundary. For the simple shear flow considered, only the S 12 (= S 21 ) component is non-zero. The theoretical value for a single sphere in a simple shear is S 12,theo = 10/3⇡µγa 3 . Time convergence is studied using a fixed grid spacing ∆=a/4.9. We note in Figure 2-a) that the steady solution on the stresslet depends on the time step. This is actually due to the assumption used in Eq. ( 30), which is inherently first-order in time. When the time step grows, this gives a poorer approximation of the momentum forcing λ through Eq. ( 29). Particle velocity remains correct -as already discussed in the previous section -but not the stresslet since it is computed using the momentum forcing λ, see Eq. [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. Although the solver eventually reaches a steady state, there is still an error in the calculation of the stresslet that explains the noted first-order time accuracy. Numerical validations -some of which will be presented hereafter -therefore show that moderate time steps are required for accurate simulations, about 10 -3 γ-1 , which are values also typically found in FCM or SD methods. This basically corresponds to a diffusional CFL number lower than roughly 50. About five grid points per particle radius (∆ ⇡ a/5) are found to be a good trade-off between accuracy and computational cost. This point is essential to keep the ability to simulate large many-particle systems while maintaining reasonable simulation time and resources.

Two smooth spheres in creeping shear flow

Let us now consider the case of two neutrally-buoyant spheres in a linear shear flow. This case is relevant to suspensions where hydrodynamic interactions are of primary importance. Two equally-sized spherical particles [START_REF] Coussot | Rheophysical classification of concentrated suspensions and granular pastes[END_REF] and ( 2) are freely suspended in a linear shear flow (with shear rate γ) and let r be the separation vector connecting the two sphere centers. Numerical parameters used are ∆=a/4.9 and ∆t ⇡ 10 -3 γ-1 (CFL d =20). The computational domain is still 20a in each direction, which has been checked to be sufficient to fulfill an infinite domain assumption.

Hydrodynamic interactions arise between spheres and alter their velocity. The theoretical expressions for the particle relative translational velocity Û = U (1) -U (2) as well as the particle rotational velocity are given by Batchelor and

Green [START_REF] Batchelor | The hydrodynamic interaction of two small freely-moving spheres in a linear flow field[END_REF] :

Ûi = ✏ i jk ! 1 j r k + r j E 1 ij -r k E 1 jk [A(r) r i r j r 2 + B(r)(δ ij - r i r j r 2 )] (80)
Ω (1) i =Ω (2) i = ! 1 i + C(r)✏ i jk E 1 kl r j r l r 2 (81)
where E 1 and ! 1 are the unperturbed rate-of-strain tensor and vorticity, respectively. The mobility functions A(r), B(r) and C(r) are known functions of the distance r between particles [START_REF] Batchelor | The hydrodynamic interaction of two small freely-moving spheres in a linear flow field[END_REF][START_REF] Dacunha | Shear-induced dispersion in a dilute suspension of rough spheres[END_REF].

Numerical simulations are performed for various configurations with different separation distances r. The computed translational and rotational velocities are then recast in terms of the mobility functions A(r), B(r) and C(r) using Eq. ( 80)-( 81) and compared with their theoretical counterparts. Figure 3 shows the obtained results without lubrication correction (open symbols) and with lubrication correction (filled symbols). Far-field hydrodynamics is well predicted, which confirms that our model correctly captures hydrodynamic interactions when particles are remote. Without lubrication correction, discrepancies between predictions and theory arise when particles almost touch, around ⇠ = r/a -2 ⇠ 0.2, which roughly corresponds to one grid spacing in this case. In this near-contact subgrid region ⇠  0.2 (delineated by the dotted line in Fig. 3), the flow field can no longer be accurately resolved and mobility functions are consequently underestimated. With the lubrication correction technique proposed, predictions now coincide with theory even in near-contact configurations. Note that computed results with or without lubrication are identical in the far-field region as lubrication correction is only activated below the lubrication barrier ⇠ lub = 0.2. The motion of particles (translation as well as rotation) is well predicted by the model even in the case of almost-touching configurations where lubrication forces become dominant. The ability to resolve very small separation distances while keeping rather coarse grids is an essential feature to simulate dense suspensions with many particles. Finally, note that particles are numerically allowed to touch (or even slightly overlap), even though there are no longer grid points between particles. This induces no singularity in a fictitious domain approach because a fluid problem is basically solved in the whole domain. When particles are allowed to move, the same configuration can also be studied in terms of relative trajectoryi.e. relative position r between particles -as a function of time. Theoretical trajectories can be obtained by a time integration of theoretical translational velocities Eq. ( 80). Both spheres are initially separated by vector r = r 0 ; they lie in the shear plane (r 0 3 = 0) and are separated by r 0 1 = -6a in the streamwise direction. Three different vertical 22 separations r 0 2 are investigated here : 2a, a and 0.5a and are presented in Fig. 4 along with the theoretical trajectories. Note that in this figure, the vertical coordinate has been stretched for the sake of clarity. The reference particle is depicted in black while the steric exclusion limit (non-overlapping region) by a dotted line. Computed trajectories are in good with theory. An important point is that they remain symmetric by virtue of the reversibility of the Stokes equations. For the case r 0 2 =2a, particles remain well separated : hydrodynamic interactions are weak and trajectories are only slightly affected. The minimal distance is well above one grid spacing so that the lubrication correction is not activated. In contrast, for r 0 2 =a or 0.5a, the minimum distance between particles can become quite small (about 10 -4 a for r 0 2 =0.5a). For these cases, the lubrication correction allows us to obtain the correct trajectories. Lubrication forces play a major role and effectively impede contact between particles. Clearly, an incorrect treatment of lubrication would result in a spurious particle overlap. Note that the authors are not aware of any similar DNS simulations since only FCM and SD methods have been used so far to compute such flows accurately. This work clearly shows that a DNS approach -like fictitious domain -is suited for low Reynolds suspension flows.

Two rough spheres in creeping shear flow

Similar computations are also conducted making allowance for particle surface roughness. It is recalled that particles come into contact when the apparent distance between their surface is smaller than the roughness height h rug , also defined as h rug =a✏ rug where ✏ rug is the non-dimensional roughness. This can be expressed as δ ij < 0 with the modified distance δ ij given by Eq. [START_REF] Rampall | The influence of surface roughness on the particle-pair distribution function of dilute suspensions of non-colloidal spheres in simple shear flow[END_REF]. When in contact, the normal contact force (which only resists compressive forces and not tension) acts to keep the inter-particle distance constant. Since roughness promotes non-hydrodynamic contact forces, it can profoundly modify suspension rheology (see for instance [START_REF] Blanc | Experimental signature of the pair trajectories of rough spheres in the shear-induced microstructure in noncolloidal suspensions[END_REF][START_REF] Rampall | The influence of surface roughness on the particle-pair distribution function of dilute suspensions of non-colloidal spheres in simple shear flow[END_REF]). Because contact forces are compressive but not tensile, they eventually result in a break of the fore-aft symmetry and the development of anisotropic microstructures. This asymmetry is clearly visible in the relative trajectories in Fig. 5 as a net displacement in the vertical direction which is seen to increase with roughness. This means that particles separate on streamlines further apart than on their approach. The DEM parameters used here are e n =e t =1 (no dissipation during contact), µ d =0.1 and k n is chosen such that the largest deformation of the rugosity is 10%. A balance between hydrodynamic and contact forces thus gives the estimation 6⇡µa 2 γ = 0.1h rug k n , which is checked a posteriori to be correct.

Note that a very similar case has been theoretically computed in [START_REF] Dacunha | Shear-induced dispersion in a dilute suspension of rough spheres[END_REF][START_REF] Blanc | Experimental signature of the pair trajectories of rough spheres in the shear-induced microstructure in noncolloidal suspensions[END_REF] and -although not reported here -the results are extremely close to our predictions which means that the proposed approach is able to address correctly the physics of contact. A triplet of spheres includes many-body interactions that are typically found in actual many-particle systems and can therefore be considered as a valuable validation for a lubrication correction method. Despite its simplicity and relevance, the case of three spheres in a shear flow has -as far as the authors know -no theoretical solutions nor reported simulation results. Two simple configurations are studied here and are sketched in Fig. 6. In each configuration, the triplet consists of equally-sized spheres suspended in a linear shear flow with shear rate γ. The cubic domain size is 30a, with a grid resolution ∆=a/4.9 and a diffusional CFL number of 20. The non-dimensional separation distance ⇠ between spheres is set to 0.01. Sphere B is located at the center of the domain. In configuration (a), the orientation angle is ✓=30 o (see Fig. 6).

As no theoretical solutions are available, reference numerical solutions for those two triplet configurations are first obtained using an alternative numerical method partly based on the commercial simulation software COMSOL r . The rigid motion is enforced by penalizing the strain tensor on the particle domain which leads to a generalized Stokes variational formulation. This basically corresponds to consider particle as a fluid with a very high viscosity. The method is implemented in the COMSOL r finite element Stokes solver and more details on the penalty method used here may be found in [START_REF] Janela | A penalty method for the simulation of fluid-rigid body interaction[END_REF][START_REF] Lefebvre | Fluid-particle simulations with freefem++[END_REF]. The COMSOL r mesh is fitted to the particles and refined on their surface. Grid convergence studies have been conducted so as to obtain a grid-independent solution. It is found that a grid spacing about 0.1a on particle surface and at least five grid points in the small 0.01a gap between particles are needed to resolve the lubrication flow properly. Preliminary computations with this penalty technique have been undertaken on two-sphere configurations and compared favorably with available theoretical solutions.

Comparisons between the COMSOL r solution (reference) and our computations are presented in Tab. 1 and Tab. 2 for the aligned and chair configurations. Translational and rotational velocities are non-dimensional using aγ and γ, respectively. As usual, subscripts (1, 2, 3) respectively denote the direction of velocity, velocity gradient and vorticity and superscripts (A, B, C) refer to the particles as depicted in Fig. 6. The average error is also specified in the bottom row of each table. Because particles are in the shear plane, velocities U 3 , ⌦ 1 and ⌦ 2 are equal to zero and are thus not given in the tables.

In both configurations, results are very encouraging especially in the aligned configuration where the average discrepancy is about 1 %. Results are not as excellent in the chair configuration but still remain acceptable. In this particular case, differences might possibly arise due to the pairwise approximation of lubrication forces. This approximation is well justified for short-range forces which is the case in the aligned configuration since interactions mostly come from a squeezing flow with a strong ⇠ -1 lubrication force. By contrast, the chair configuration rather induces a shearing flow between particles with a weaker log ⇠ singularity.

Sphere rolling down an inclined plane

The motion of a sphere rolling down an inclined planar surface due to gravity provides a convenient configuration for validating the coupling between lubrication and contact forces. A sphere of radius a, having microscopic asperities of uniform height h rug , is placed on a wall in a gravity field g = (g sin ✓, -g cos ✓, 0) to mimic an inclined plane with inclination angle ✓. Creeping flow conditions are assumed. Figure 7 particle velocity along with the experimental results from Davis et al. [START_REF] Davis | Solid-solid contacts due to surface roughness and their effects on suspension behaviour[END_REF]. Translational and rotational velocities are non-dimensionalized using U s sin ✓ and U s sin ✓/a, respectively, where U s = 2a 2 ∆⇢g/9µ is the Stokes settling velocity. The roughness height h rug = 0.0015a and the friction coefficient µ d = 0.12 are taken from [START_REF] Davis | Solid-solid contacts due to surface roughness and their effects on suspension behaviour[END_REF]. Simulations and experiments are in very good accordance showing that the model is capable of describing accurately the expected physics. For small angles of inclination, the sphere rolls without slipping so that U = a⌦. As inclination is increased, the critical angle (✓ cr ⇠ 15 o ) is reached and the sphere begins to slip with partial rolling, i.e. U > a⌦. The limit values U=0.231 and ⌦=0.036, for ✓ = 90 o , correspond to a motion completely dictated by lubrication and can be obtained using theoretical resistance functions for wall-sphere interactions. This test shows that a suitable contact model is thus needed, in addition to lubrication, to resolve the correct physics, particularly the rolling/slipping behavior. Accounting for roughness is important as well since it controls the surface-to-surface separation and the magnitude of lubrication forces accordingly. A perfectly smooth sphere would -at least theoretically -stick to the plane and not move due to the singularity of the hydrodynamic resistance functions at contact.

High-frequency shear viscosity in random arrangements

This simulation considers a many-particle system with N neutrally buoyant spherical particles (radius a) in a cubic domain of volume V with volume fraction φ = 4⇡a 3 N/3V. Particles are placed randomly using a Monte-Carlo packing method as described by Torquato [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]. A shear flow of magnitude γ is generated by two opposite velocities at top and bottom boundaries ; other boundaries being periodic. The size of the computational domain is 30a and the grid spacing is ∆=a/4.9 so that the number of grid points is 147 3 . Creeping flow conditions are again assumed. The relative high-frequency shear viscosity µ r can be defined by [START_REF] Guazzelli | A Physical Introduction to Suspension Dynamics[END_REF] µ where the stresslet tensor S (i) for particle i can be obtained from Eq. ( 66). The high-frequency limit here means that the microstructure is frozen and viscosity is consequently computed for a fixed (i.e. non time-evolving) particle configuration. For each prescribed volume fraction φ, ranging from 0.05 to 0.6, the shear viscosity is obtained by averaging over ten different random particle configurations. For φ=0.6, the configuration contains roughly 3,900 particles. Results are shown in Fig. 8 along with similar computations from Sierou and Brady [START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF] (using SD) and Yeo and Maxey [START_REF] Yeo | Dynamics of concentrated suspensions of non-colloidal particles in couette flow[END_REF] (using FCM). Our predictions appear to be quite close to simulations from the literature. Some slight differences are noticed with Sierou and Brady [START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF] for high volume fractions. This can be anticipated since present computations account for upper and lower walls whereas simulations from Sierou and Brady [START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF] are performed for an unbounded domain. A divergence of viscosity as φ ! φ max is noted. Our results were fitted to a Krieger-Dougherty law which is widely employed in the field of rheology and reads [START_REF] Krieger | A mechanism for non-newtonian flow in suspensions of rigid spheres[END_REF]

r = 1 + 1 γµV N X i=1 S (i)
µ r = (1 - φ φ max ) -[µ]φ max ( 83
)
The obtained best-fit parameters are [µ]=2.32 and φ max =0.68 which is not far from the usual φ max ⇡ 0.64 for a monodisperse system of hard spheres. In this dense regime, the distance between particles can be quite small (as low as 10 -6 a), so that the lubrication system Eq. ( 55) is ill-conditioned. Nonetheless, the preconditioned conjugate gradient solver behaves well and requires at most 20-25 iterations to solve the system with a relative tolerance level of 10 -6 . This simulation shows that actual concentrated suspensions can be handled with our method and that quantitative predictions of rheological quantities, such as viscosity, can be obtained. Once again, this suggests that DNS can be a valuable alternative to SD or FCM to address the rheology of dense suspensions. As the code is not yet parallelized, the computation was run on a single core of a Linux IBM Idataplex cluster, based on Intel Xeon 5600 processors at 3,1 GHz, and needed 25 GB of memory. It took about 72 h of CPU time for 24153 iterations (0.74 µs/iteration/grid point). This is a reasonable restitution time for such a simulation and it is expected that computations including much more particles will become feasible thanks to the forthcoming parallelization.

Figure 9 shows some snapshots of particle configuration for different non-dimensional times t/t s . The core of the suspension first swiftly settles followed by the rest of the particles. However, particles near the walls fall very slowly due to the strong wall/particle interactions. Indeed, it has been seen previously (see section 5.5) that a particle very close to a wall (separation distance of about the roughness h rug =0.0015a) had a settling velocity of only 0.23 times its Stokes velocity. This is even more pronounced for particles trapped in a corner of the domain : as can be seen in For this simulation, Tab. 3 presents the CPU performance obtained by profiling. It can be noted that the most important computational effort lies in solving the large linear system from the implicit momentum equation. Similarly, the linear system for the Poisson equation on pressure also needs much CPU usage. They both represent more than 42 % of the total time and are independent of the number of particles. Body-force advection as well as computation of particle velocities, level-sets or rigid body motions are another 36 % of the computation time and, interestingly, scale linearly with the number N of particles. Only the lubrication system is not linear in N since it rests on a conjugate gradient solver but it represents less than 10 % of the total load. Note that it also depends on the configuration, more specifically on the size of clusters of particles separated by less than the lubrication barrier.

Conclusions

In this paper, we have proposed a Lagrange multiplier-free fictitious domain approach adapted for the simulation of dense suspensions. The particle-fluid solver is close to some recent works [START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF][START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF]. The main difference however is that our approach is fully Eulerian and does not need any Lagrangian nodes attached to particles.

Lubrication forces as well as contact forces play a major role in low-Reynolds suspensions and have here been addressed in much detail. Because they are very short-range, lubrication forces can not be fully resolved with the typical grids used and consequently require an additional model. The method proposed is inspired from Stokesian Dynamics but modified here so as to comply with our fictitious domain solver. The non-resolved (subgrid) part of lubrication is first identified by subtracting tabulated numerically-resolved lubrication interactions from their theoretical counterparts. This subgrid part is subsequently considered in a resistance problem to correct particle velocities.

This leads to an ill-conditioned resistance linear system solved successfully using an IC(0)-preconditioned conjugate gradient method. This lubrication method is quite general and can be used in any flow solver. It has also been adapted to account for solid walls in order to address bounded flows. Note that although the numerical method can handle non-spherical particles or non-Newtonian fluids, the lubrication method is strongly linked with available lubrication theory and is therefore restricted to spheres in a Newtonian fluid.

Unlike most reported suspension simulations, contact forces are here computed using a physical model based on the discrete element method (DEM) which is widely used in granular matter physics. The DEM is implemented so that it remains strongly coupled with lubrication forces. To this end, contact forces are simply added as a background force in the lubrication problem. Roughness is also included in the DEM by altering the definition of the distance between particles. Some detailed validations are provided for one-particle, two-particle and three-particle configurations in a linear shear flow. The good correlations obtained indicate that the method is well adapted for the simulation of particulate flows ; it particularly gives good results even for flows governed by lubrication. Notably, simulations of a single particle on an inclined wall as well as rough particles interacting in a shear flow suggest that lubrication and contact are both accurately described. Simulation of many-particle systems up to O(10 4 ) particles can be conducted in a reasonable CPU time, even though the code is sequential so far . It is expected that parallelization may open ways to simulations of O(10 6 ) particles that will allow for a better knowledge of suspension physics.

This study shows that actual concentrated simulations can be handled by the present method and that quantitative predictions of rheological behavior, such as viscosity, are possible. This suggests that direct numerical simulation is a valuable alternative to address the rheology of dense suspensions.

Figure 1

 1 Figure1presents a comparison between the theoretical solution Eq. (79) and the computed non-dimensional velocity u 2 /γa on the centerline (x 2 =0) and in the shear plane (x 3 =0). The particle is located at x=0. The numerical results (open symbols in Fig.1) clearly show a very good correlation with the analytic solution outside the particle (|x 1 |/a > 1) but also inside the particle (|x 1 |/a < 1). In the latter region, the linear profile is imposed by the rigid body motion due to a rotation with rate -γ/2.e 3 , from which a velocity |u 2 |/γa=0.5 is found at |x 1 |/a=1. With the numerical

Figure 1 :

 1 Figure 1: Centerline fluid velocity u 2 /γa against position x 1 /a for a single sphere at the origin (x=0). Analytical (solid line) and computed solution with momentum forcing (⌥) and without momentum forcing (⌅).

Figure 2 -

 2 scheme is about one, which is expected according to theory. Figure2-b) presents the error in space and time which is evaluated by reducing the grid size and time steps simultaneously. A diffusional CFL number of 50 is held fixed. The overall order is about 2.

Figure 2 :

 2 Figure 2: Error in the computed stresslet against time step (a) and grid resolution (b).

Figure 3 :

 3 Figure 3: Simulated (symbols) and theoretical (solid line) mobility functions as a function of non-dimensional separation distance ⇠ = r/a -2: A (triangles), B (squares) and C (circles). Open symbols (4,⇤,⌥) are without lubrication correction ; filled symbols (N,⌅, ) include lubrication correction. The dotted line delineates the lubrication barrier.

Figure 4 :

 4 Figure 4: Relative trajectories of two spheres in a shear flow for initial vertical separations r 0 2 =2a, a and 0.5a : theory (solid line) and computations (symbols). The black circle depicts the reference particle and the dotted line delineates the non-overlapping region.

Figure 5 :

 5 Figure 5: Relative trajectories of two rough spheres in a shear flow for different non-dimensional roughness ✏ rug

Figure 6 :

 6 Figure 6: Considered triplet configurations : (a) aligned configuration with orientation angle ✓=30 o ; (b) chair configuration. The non-dimensional distance is ⇠=0.01

Figure 7 :

 7 Figure 7: Non-dimensional translational U (diamonds) and rotational Ω (circles) velocities of a sphere moving down an inclined plane : computations (open symbols) and experimental results [81] (filled symbols)

Figure 8 :

 8 Figure 8: Predicted relative high-frequency viscosity against volume fraction for a random arrangement. Present results are fitted to a Krieger-Dougherty law Eq. (83) with [µ]=2.32 and φ max =0.68

Fig 9 f

 9 Fig 9 f), those particles have a very low settling velocity. The simulation was therefore stopped at t/t s =1.61 before all particles settled.

Figure 9 :

 9 Figure 9: Particle configuration for different non-dimensional time : a) t/t s =0.00 ; b) t/t s =0.06 ; c) t/t s =0.14 ; d) t/t s =0.32 ; e) t/t s =0.61 ; f)

  ), the theoretical resistance relation for the force and torque is[40]

				2
	2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 T 1 F 1 F 2 T 2	3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5	= µR 2B,theo	6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4

  A 11 A 12 (B 11 ) t (B 21 ) t G 11 G 21 A 12 A 22 (B 12 ) t (B 22 ) t G 12 G 22

						7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7	(57)
	with
		2					3
	R 2B,theo =	6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4	B 11 B 12 B 21 B 22	C 11 C 12	C 12 C 22	H 11 H 21 H 12 H 22	7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5

Table 1 :

 1 presents the computed translational and rotational Reference and present simulations for the aligned configuration

		Reference Computations		Reference Computations
	U (A) 1	0.322	0.321	U (A) 1	1.765	1.765
	U (A) 2	0.463	0.467	U (A) 2	0.132	0.125
	Ω (A) 3	-0.388	-0.389	Ω (A) 3	-0.669	-0.669
	U (B) 1	0.000	<0.001	U (B) 1	0.135	0.136
	U (B) 2	0.000	<0.001	U (B) 2	0.132	0.125
	Ω (B) 3	-0.341	-0.339	Ω (B) 3	-0.500	-0.507
	U (C) 1	-0.322	-0.323	U (C) 1	0.134	0.135
	U (C) 2	-0.463	-0.463	U (C) 2	-0.246	-0.250
	Ω (C) 3	-0.388	-0.380	Ω (C) 3	-0.329	-0.324
	Err (%)		0.6	Err (%)		1.9

Table 2 :

 2 Reference and present simulations for the chair configuration

Table 3 :

 3 CPU performance for the sedimentation of 10,648 particles

		% CPU time
	Momentum solve	28.6
	Body-force advection	20.4
	Pressure solve	13.7
	Lubrication solve	9.8
	Particle velocities	6.9
	Level-set computation	5.7
	Distance computation	3.3
	Rigid motion forcing	2.6
	Collision (DEM)	0.1
	Others	8.9

[START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]

Acknowledgments

This work was funded by the French Procurement Agency (DGA)