Stany Gallier
email: stany.gallier@herakles.com

Francois Peters

Laurent Lobry

Rheology of sheared suspensions of rough frictional particles

This paper presents three-dimensional numerical simulations of non-Brownian concentrated suspensions in a Couette flow at zero Reynolds number using a fictitious domain method. Contacts between particles are modelled using a DEM-like approach (Discrete Element Method), which allows for a more physical description, including roughness and friction. This work emphasizes the effect of friction between particles and its role on rheological properties, especially on normal stress differences. Friction is shown to notably increase viscosity and second normal stress difference |N 2 | and decrease |N 1 |,i n better agreement with experiments. The hydrodynamic and contact contributions to the overall particle stress are particularly investigated. This shows that the effect of friction is mostly due to the additional contact stress since the hydrodynamic stress remains unaffected by friction. Simulation results are also compared with experiments, such as normal stresses or effective friction coefficient µ(I v ), and the agreement is improved when friction is accounted for. This suggests that friction is operative in actual suspensions.

Introduction

Concentrated suspensions of small particles in low Reynolds number flows are ubiquitous in industry (food transport, cosmetic products, civil engineering, to mention but a few) as well as in biological or natural flows (blood, mud or lava flows). Fresh concrete or uncured solid rocket fuel are two examples of industrial concentrated suspensions in which a very high particle volume fraction is desired while keeping correct rheological properties and flowing behaviour. Such dense suspensions do exhibit an intricate physics which is hitherto far from being understood completely. This complexity partly arises from the wide variety of interactions between particles : Brownian, colloidal, hydrodynamic, frictional or collisional (see [START_REF] Coussot | Rheophysical classification of concentrated suspensions and granular pastes[END_REF] for a review on those different regimes) as well as from the physical properties of particles (roughness, shape, size distribution,...). Even the case of non-Brownian non-colloidal single-sized spherical particles embedded in a Newtonian fluid -which will be the system investigated in this study -is likely to show complex non-Newtonian behaviours [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]).

The rheological behaviour of non-Brownian suspensions was initially described by an effective suspension viscosity ⌘ s , and considerable effort focused on determining this viscosity as a function of particle volume fraction φ. Experimental studies have confirmed that suspension viscosity ⌘ s increases with φ and diverges as φ approaches a maximal fraction φ m about 0.6 [START_REF] Ovarlez | Local determination of the constitutive law of a dense suspension of noncolloidal particles through magnetic resonance imaging[END_REF]. Despite some differences among the various experimental results, the relation between viscosity and particle volume fraction is now reasonably well documented [START_REF] Stickel | Fluid mechanics and rheology of dense suspensions[END_REF]. By contrast, the normal stresses Σ xx ,Σ yy , and Σ zz (where x, y, and z refer to the direction of velocity, velocity gradient, and vorticity, respectively) are much less studied and available results are still controversial. From experiments in parallel-plate and Couette geometries, [START_REF] Zarraga | The characterization of the total stress of concentrated suspensions of noncolloidal spheres in newtonian fluids[END_REF] measured negative normal stress differences N 1 =Σ xx -Σ yy and N 2 =Σ yy -Σ zz with similar dependence on the volume fraction. The second normal stress difference |N 2 | is found to be larger than |N 1 | with N 2 /N 1 ⇡ 3.6. Later measurements by [START_REF] Singh | Experimental measurements of the normal stresses in sheared stokesian suspensions[END_REF] using similar devices have confirmed that N 1 and N 2 are both negative but with a ratio N 2 /N 1 which depends on φ. [START_REF] Dai | Viscometric functions for noncolloidal sphere suspensions with newtonian matrices[END_REF] used a combination of a parallelplate rheometer and the open semi-circular trough method to measure N 1 and N 2 .T h e y found a good agreement with previous works with negative normal stress differences and |N 1 |⌧| N 2 |. Those three studies measured negative normal stress differences with |N 1 | smaller than |N 2 |. On the contrary, some recent experimental works have obtained very different results, especially on the magnitude of |N 1 |. By studying the deformation of the free surface of a suspension flowing in a narrow inclined channel, [START_REF] Couturier | Suspensions in a tilted trough: second normal stress difference[END_REF] concluded that N 1 was very small compared to N 2 and within the experimental uncertainty, they were unable to determine the sign of N 1 . Likewise, a detailed experimental study by [START_REF] Dbouk | Normal stresses in concentrated non-brownian suspensions[END_REF] using a parallel-plate geometry has found this time that N 1 is positive over the entire range of volume fraction investigated 0.2  φ  0.47. The value of N 2 however remains as in the other studies. Numerical simulations of concentrated suspensions using Stokesian Dynamics [START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF] or Force-Coupling Method (Yeo & Maxey 2010a) show that N 1 and N 2 are both negative with N 1 ⇡ N 2 , which actually does not match any experiments.

Numerical simulation is instrumental in shedding light on the complex physics of suspensions. Among the most widely employed method for low Reynolds number suspensions is the 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] Brady | Stokesian dynamics[END_REF]. SD has provided among the most relevant results in the field of suspensions [START_REF] Drazer | Deterministic and stochastic behaviour of non-brownian spheres in sheared suspensions[END_REF][START_REF] Drazer | Microstructure and velocity fluctuations in sheared suspensions[END_REF][START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF]. This technique resembles a molecular dynamics method in which the driving forces are the hydrodynamic interactions, possibly complemented with Brownian or collisional forces. Hydrodynamic interactions are written as a truncated multipole expansion and further split into a far-field multi-body contribution and near-field lubrication contribution. Another popular method relevant to low Reynolds number suspension flows is the Force-Coupling Method (FCM) [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF][START_REF] Lomholt | Force-coupling method for particulate two-phase flow: Stokes flow[END_REF]Yeo & Maxey 2010c,a,b). Similarly, hydrodynamic interactions are written using a multipole expansion in which the standard Dirac delta function is replaced by a localized force envelope. The latter is incorporated as a body-force in the Stokes equations that are subsequently solved using classical flow solvers. More recently, Direct Numerical Simulations (DNS) have emerged as a valuable alternative to the aforementioned SD and FCM methods, such as Lattice-Boltzmann methods (Ladd 1994a,b;[START_REF] Ladd | Lattice-Boltzmann simulations of particle-fluid suspensions[END_REF] or fictitious domain methods [START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[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] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF]Wachs 2009;[START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF]. Because they directly solve the governing equations (Navier-Stokes or Stokes equations) without any further assumptions other than numerical approximations, they are more general and appropriate to deal with inertial flows, arbitrary-shaped particles or non-Newtonian fluid. For any methods however, dealing with concentrated suspensions requires a detailed modelling of lubrication forces since the average separation distance becomes extremely small. Lubrication forces arise between particles in near-contact because of the draining of interstitial fluid in the gap, and are known to be singular in the limit of touching particles. Consequently, lubricationinduced dissipation increases and modifies the rheology of suspensions. A last essential physical ingredient for accurate simulations of suspension is the modelling of collisions, or contacts, between particles. Contacts inevitably occur in concentrated stokesian suspensions despite lubrication because of particle roughness. It is now well-accepted that contacts are responsible for fore-aft asymmetry and non-zero normal stress differences. They are also the source of irreversibility and chaos in sheared suspensions [START_REF] Metzger | Irreversibility and chaos: Role of long-range hydrodynamic interactions in sheared suspensions[END_REF][START_REF] Metzger | Irreversibility and chaos: Role of lubrication interactions in sheared suspensions[END_REF].

As stressed previously, there is a large discrepancy between experiments on N 1 and N 2 , but also between experiments and simulations, and this has not received much attention so far. On the other hand, normal stresses are crucial since they determine the migration of particles in suspensions subjected to an inhomogeneous shear field [START_REF] Nott | Pressure-driven flow of suspensions: simulation and theory[END_REF][START_REF] Mills | Rheology and structure of concentrated suspensions of hard spheres. shear induced particle migration[END_REF][START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF][START_REF] Nott | The suspension balance model revisited[END_REF][START_REF] Lhuillier | Migration of rigid particles in non-brownian viscous suspensions[END_REF]. This paper therefore intends to address this issue by considering numerical simulations of suspensions accounting for friction and roughness. Friction implies an additional tangential force during contact that is likely to alter normal stresses. Using Stokesian Dynamics simulations, [START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF] investigated the effect of sliding friction on the rheological behaviour of a concentrated suspension at φ=0.4 and showed that friction slightly changes the microstructure which, in turn, modifies N 1 and N 2 . For a friction coefficient µ d =0.5, they noticed an overall 50 % reduction in |N 1 | and a similar increase in |N 2 | as well as a moderate rise in the suspension viscosity ⌘ s (about 10 %). The study was however limited to a single volume fraction and the model used was deliberately simple (for instance, particles are supposed to slip irrespective of the magnitude of the normal force). This effect of friction was also confirmed by [START_REF] Wilson | Shear stress of a monolayer of rough spheres[END_REF]; [START_REF] Davis | Solid-solid contacts due to surface roughness and their effects on suspension behaviour[END_REF] who reported a moderate effect of friction with a slight increase in the viscosity and decrease in |N 1 |. The friction coefficient is known to depend strongly on normal force, surface roughness, etc. and this could possibly explain the discrepancy between experiments. Note that friction is also believed to play a major role in jamming and discontinuous shear-thickening [START_REF] Seto | Discontinuous shear thickening of frictional hard-sphere suspensions[END_REF][START_REF] Wyart | Discontinuous shear thickening without inertia in dense nonbrownian suspensions[END_REF]. A second potential cause for the noted scatter in normal stresses is the presence of particle roughness. Particles are never perfectly smooth and a dimensionless roughness height of about 10 -3 ⇠ 10 -2 is generally reported [START_REF] Smart | Measurement of the hydrodynamic surface roughness of noncolloidal spheres[END_REF]. Such surface asperities promote an early contact and impose the surface-to-surface separation between particles and the magnitude of lubrication stress consequently. Surface roughness can therefore significantly modify the rheological properties and the microstructure of suspensions, as confirmed by many numerical studies [START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF][START_REF] Drazer | Deterministic and stochastic behaviour of non-brownian spheres in sheared suspensions[END_REF][START_REF] Drazer | Microstructure and velocity fluctuations in sheared suspensions[END_REF][START_REF] Dacunha | Shear-induced dispersion in a dilute suspension of rough spheres[END_REF]. Theoretical studies in dilute regimes show that large roughness can decrease viscosity as well as increase |N 1 | and |N 2 | [START_REF] Wilson | An analytic form for the pair distribution function and rheology of a dilute suspension of rough spheres in plane strain flow[END_REF][START_REF] Wilson | Shear stress of a monolayer of rough spheres[END_REF][START_REF] Davis | Solid-solid contacts due to surface roughness and their effects on suspension behaviour[END_REF][START_REF] Zarraga | Normal stress and diffusion in a dilute suspension of hard spheres undergoing simple shear[END_REF]. [START_REF] Zarraga | Normal stress and diffusion in a dilute suspension of hard spheres undergoing simple shear[END_REF] found for small roughness that |N 1 | > |N 2 | due to a depletion of particles in the extension quadrant. The effect is reversed for large roughness due to the additional stress induced by contact forces and N 2 becomes the prevailing normal stress difference. In dense regimes however, Stokesian Dynamics simulations from [START_REF] Wilson | Shear stress of a monolayer of rough spheres[END_REF] attested that |N 1 | decreased with increasing roughness height, which suggests that the role of roughness is not completely understood so far.

This paper intends to contribute to the physics of sheared non-Brownian suspensions with emphasis on rheology -especially normal stress differences -using numerical simulations. We aim at providing insights that could help understand the discrepancy between computations and experiments as well as between experiments themselves. The study specifically addresses the role of friction forces between particles. In present work, we make use of a fictitious domain approach as detailed in [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF]. Our method explicitly solves long-range hydrodynamics and incorporates a modelling of lubrication interactions as well as DEM (Discrete Element Method)-based contact forces. This allows for a straightforward implementation of tangential friction forces or particle roughness. Section 2 presents a brief description of the fictitious domain approach used. In section 3, we consider the simulation of a pair of particles in a shear flow with an eye to demonstrating that friction and roughness are correctly modelled. Section 4 presents suspension simulation results in a Couette geometry and highlights the important role of friction forces during contact.

DNS simulation : fictitious domain method

This section briefly describes the numerical method used ; more details can be found in a previous paper [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF]. In a fictitious domain method, solid particles are supposed to be filled with a fluid having the same properties as the actual fluid. From a computational viewpoint, this means that a classical fluid problem is solved in the whole domain. Particles are thus considered as some regions of the fluid constrained to have a rigid body motion.

Governing equations

Let D be the whole computational domain including a fluid domain D f and a solid particle domain D p with D = D f ⊕D p . The solid domain is made up of particles P that are supposed to be rigid and homogeneous. In this study, the fluid is assumed incompressible and Newtonian and is governed by the Stokes equations (although Navier-Stokes can also be readily considered). Mass and momentum conservation equations are written in the whole domain D and read

r•u = 0 (2.1) r•Σ + ⇢λ = 0 (2.2)
where ⇢ and u are the fluid density and velocity, respectively, while λ is a momentum forcing term used to enforce the rigid body motion inside particles. Consequently, λ is non-zero in the particle domain D p and zero elsewhere. For a Newtonian fluid, the stress tensor Σ reads

Σ = -pI +2⌘E (2.3)
where p is the pressure, ⌘ the fluid viscosity and E the rate-of-strain tensor E = 1 2 (ru + ru T ). The fluid velocity inside each particle must comply with a rigid body motion, so that

u = U + Ω ⇥ (x -x g )i n D p (2.4)
where U and Ω stand for the particle translational and rotational velocities and x g is the position of the centre of gravity of the particle. Particle motion is given by Newton's equations. Assuming spherical particles and neglecting inertia, they read

F h + F c + F e = 0
(2.5)

T h + T c + T e = 0 (2.6)
where forces F and torques T are decomposed into their hydrodynamic part (h), contact part (c) and external part (e) that includes any external forces, such as gravity. In absence of inertia, hydrodynamic interactions can be defined as [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF][START_REF] Yu | A direct-forcing fictitious domain method for particulate flows[END_REF])

F h = -⇢ Z P λ dx (2.7) T h = -⇢ Z P (x -x g ) ⇥ λ dx (2.8)
Contact forces and torques are modelled as usually in Discrete Element Method (DEM), as addressed in a subsequent section.

Numerical procedure

The system Eq. ( 2.1)-( 2.8) forms the set of governing equations to be solved. A fractionalstep time scheme is used to decouple this system into two sub-problems. The first subproblem Eq. ( 2.1)-( 2.3) is a standard fluid problem. It is solved using finite differences on a staggered Cartesian grid and incompressible equations are then resolved by a standard projection method [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF]) with implicit time-stepping for the diffusive terms (Crank-Nicolson scheme). The resulting linear systems are inverted using a geometric multigrid technique. At each iteration n, solving this fluid sub-problem yields the new fluid velocity u n+1 and pressure p n+1 . The particle sub-problem Eq. (2.4)-(2.8) consists in enforcing the rigid body motion inside particles, as well as updating the new forcing term λ n+1 and eventually computing particle velocities. Details are skipped in the frame of this overview and may be found in [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF]. Once particle velocities are known, their position can be tracked by solving dX dt = U (2.9)

where X and U are the vectors of particle positions and velocities, respectively. This equation is solved using a second-order Adams-Bashforth scheme.

Lubrication correction

Lubrication forces play a major role in the rheology of concentrated suspensions. Since they are very short-range in nature, they can usually not be fully resolved with the typical grids used in the numerical model and they consequently require an ad hoc model. Our lubrication model rests on the ideas already used in SD [START_REF] Brady | Stokesian dynamics[END_REF][START_REF] Durlofsky | Dynamic simulation of hydrodynamically interacting particles[END_REF]) and FCM (Yeo & Maxey 2010a).

Let us consider a system of N p spherical particles suspended in a linear Stokes flow and let U be the 6N p vector of translational/rotational velocities U =( U , Ω) T and F =(F , T ) T the 6N p 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 [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF])

F = R FU • (U ∞ -U )+R FE : E ∞ (2.10)
The key idea is to split the resistance matrix describing the hydrodynamic interactions into long-range interactions -explicitly resolved by the numerical model -and a shortrange lubrication contribution that can not be resolved numerically since it is subgrid 2.11) where the tilde (.) denotes the explicitly resolved part and the superscript sub refers to the subgrid unresolved part of the interaction. It can be shown that the actual lubricationcorrected particle velocity U can be obtained from the numerically resolved velocity Ũ by adding in the numerical procedure an external force F lub given as [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF])

R ⇡ R + R sub (
F lub = R sub FU • (U ∞ -U )+R sub FE : E ∞ (2.12)
This force/torque represents the portion of hydrodynamic interactions that can not be resolved by the numerical approach and is directly included in Eq. (2.5)-(2.6) as an external force and torque. The subgrid resistance matrices R sub FE and R sub FU are computed by subtracting the two-body resistance matrix R2B -obtained numerically on two-sphere configurations -from the exact theoretical two-sphere resistance matrices R 2B theo known from lubrication theory [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF]. For a many-particle system, matrices R sub are constructed assuming a pairwise additivity of forces.

The hydrodynamic stresslet S h is corrected from lubrication as well using a similar procedure. The deviatoric stresslet is written in resistance form as

S h = R SU • (U ∞ -U )+R SE : E ∞ (2.13)
and is similarly decomposed into a resolved and a subgrid part

S h = S + R sub SU • (U ∞ -U )+R sub SE : E ∞ (2.14)
where S corresponds to the resolved stresslet computed by the numerical method. It 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 = - ⇢ 2 Z P [(x -x g ) ⌦ λ + λ ⌦ (x -x g )] dx (2.15)
Subgrid resistance matrices R sub SU and R sub SE are obtained as described previously while theoretical expressions are found in [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF]. Finally, a similar correction procedure is also applied to the trace of S h -which represents the hydrodynamic contribution to particle pressure Π -using the theoretical resistance functions from [START_REF] Jeffrey | The pressure moments for two rigid spheres in low-reynolds-number flow[END_REF]. Particle-wall lubrication interactions are also included using a similar strategy.

Despite its complexity, an accurate treatment of lubrication is crucial for quantitative simulations of dense suspensions. Some computations without lubrication correction were attempted but generally failed unless extremely small time steps were considered. In the latter case, very noisy signals were observed. This is likely to result from large particle relative velocities implying very intense collisions. Likewise, a too crude lubrication correction model (e.g. neglecting logarithmic singularities) also involves discrepancies, such as an underestimation of viscosity.

Contact modelling

In order to implement a more physical modelling of contact between particles, contact forces are modelled as in a Discrete Element Method (DEM), which is a very popular method in the field of granular physics [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 pair of spherical particles P i and P j (of radius a) undergoing contact, the contact force F c exerted by particle P j on particle P i is classically decomposed into its normal F c n and tangential F c t components : The normal contact force is modelled using a Hertz law

F c = F c n + F c t (2.16)
F c n = -k n |δ| 3/2 n (2.17)
in which δ = krk-2a is the overlap distance with r = x jx i , and n is the normal vector n = r/krk. Surface roughness can be readily accounted for in the model, although in a simple way. Assuming sparse asperities of size h r , contact is supposed to occur whenever krk 6 2a + h r . An easy way to implement roughness thus consists in defining a modified overlap distance as δ 0 = δ -h r . Hence, contact occurs if δ 0  0. Figure 1 presents a sketch depicting this roughness modelling between two particles. Note that lubrication forces are however still evaluated with actual distance δ = δ 0 + h r . The normal stiffness k n in Eq. (2.17) can be expressed theoretically using particle mechanical properties (Young modulus and Poisson coefficient). For real materials however, the resulting stiffness is extremely high, which is prejudicial to numerical stability unless the time step is reduced drastically. For numerical reasons, it is therefore preferred to prescribe a lower stiffness, which is here given in terms of an average roughness deformation ✏ = | δ0 |/h r . Balancing contact and hydrodynamic forces (in dilute regimes) yields 6⇡⌘ γa 2 = k n (✏h r ) 3/2 from which k n can be estimated. In the forthcoming computations, we set ✏=0.05, i.e. the stiffness k n is such that the roughness deformation is 5 % under the hydrodynamic force 6⇡⌘ γa 2 . Parametric simulations on ✏ (and so, k n )h a v es h o w n that there are no effects on rheology as long as it is sufficiently small, typically below 0.2. However, the actual force exerted on particles is larger than 6⇡⌘ γa 2 for dense systems. This means that the roughness deformation can exceed the prescribed value of 0.05. Its value was checked systematically and was found to be always lower than 0.1. This limited variation in ✏ is related to the non-linear nature of the Hertzian force which scales as δ 3/2 . The relative importance of contacts with respect to hydrodynamics can be estimated using parameter Γ, which is analogous to a Péclet number :

Γ= 6⇡⌘ γa 2 F 0 (2.18)
where F 0 is a typical scale for the contact force. Taking F 0 as the force needed to deform completely a roughness

F 0 = k n h 3/2 r
gives Γ=✏ 3/2 ⇡ 10 -2 .S i n c e✏ is fixed, this means that the relative contribution of contact remains similar regardless of roughness height.

The physical relevance of keeping a constant Γ is questionable. However present work primarily focuses on friction and the effect of Γ will not be investigated in detail here. Its role on shear-thickening for instance was studied in other works, e.g. see [START_REF] Seto | Discontinuous shear thickening of frictional hard-sphere suspensions[END_REF] for very dense systems. We have already mentioned that parametric simulations (for φ bulk =0.4 and frictionless particles) show no effects of ✏ if it is small enough. This is tantamount to varying Γs i n c e Γ=✏ 3/2 . Results are found independent of Γi nt h e range 10 -3 ⇠ 10 -1 . Values below 10 -3 have not been attempted but are expected to give similar results because we stay in the asymptotic regime of nearly hard spheres. For Γ above 10 -1 , a shear-thickening behaviour is however noticed : since roughness is softened, the interparticle distance is reduced and the viscosity increases because of lubrication-induced dissipation. All the results presented in this paper are obtained with Γ ⇡ 10 -2 but are valid at least in the range 10 -3 ⇠ 10 -1 .K e e p i n g Γ constant implies

that k n h 3/2 r
is also held constant for the same shear rate. Obviously, this is irrelevant for zero roughness (h r =0), causing an infinite stiffness k n . In that case, an arbitrary high value of k n should be prescribed instead. However, dense suspensions of perfectly smooth particles were not investigated in this work.

The tangential force is given by

F c t = -k t Υ (2.19)
in which Υ is defined by integrating the slip velocity U s during the contact

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

U s = U i -U j -[(U i -U j ) • n] • n +(aΩ i + aΩ j ) ⇥ n (2.21)
Using the classical Amontons-Coulomb law of friction, the actual tangential force is modified if it exceeds the friction limit µ d |F c n | and is then given by

F c t = µ d |F c n | F c t |F c t | (2.22)
where µ d is the dynamic friction coefficient. This work does not consider the static friction coefficient µ s and we implicitly set µ s =µ d . The tangential stiffness k t is linked with the normal stiffness k n and a classical choice is . 2001;[START_REF] Shäfer | Force schemes in simulations of granular materials[END_REF]. Finally, the corresponding contact torque is

k t /k n =2 δ 1/2 /7( S i l b e r tet al
T c = an ⇥ F c (2.23)
Contact forces also induce an additional contact stresslet which is given for a particle as

S c = 1 2 (F c ⌦ r 2 + r 2 ⌦ F c ) (2.24)
This stresslet is not traceless and its trace determines the contribution of contact to particle pressure. Also note that since F c ⌦ r may not be symmetric, there exists an antisymmetric contact contribution to the particle stress due to contact torques. Contact torques are balanced by hydrodynamic torques by virtue of Eq. (2.6).

A pair of rough frictional spheres in a shear flow

The numerical model presented previously has been validated thoroughly on different configurations in [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF]. However, before turning to the suspension simulations in §4, we here consider the case of a sheared pair of rough frictional spheres. This preliminary computational study can be viewed as an additional relevant validation since results will be compared to a reference theoretical solution. It also gives a foretaste of the forthcoming suspension simulations and allows the reader to become familiar with rough frictional particles interacting in a shear flow.

Two equally-sized spherical particles (of radius a) are freely suspended in a simple shear with separation vector r connecting the two sphere centres. A shear flow is generated by moving the upper and lower boundaries with opposite velocities while the other boundary conditions (in x and z directions) are periodic. Figure 2 presents a sketch of this two-particle configuration. A cubic computational domain of size 20a in the three directions is used. The spheres are initially separated by r init =(-6a, 0.5a, 0). Due to the imposed shear flow, the spheres will hydrodynamically interact and possibly touch in the case of rough spheres. A theoretical reference solution is computed using the approach already presented in other works [START_REF] Dacunha | Shear-induced dispersion in a dilute suspension of rough spheres[END_REF][START_REF] Zarraga | Normal stress and diffusion in a dilute suspension of hard spheres undergoing simple shear[END_REF][START_REF] Metzger | Irreversibility and chaos: Role of lubrication interactions in sheared suspensions[END_REF]. The relative trajectories are integrated in time using the theoretical resistance functions R 2B theo . We have implemented this approach using Matlab R .I n addition, this model is supplemented with a DEM-like contact model such as described in §2.4 to allow us to consider theoretical solutions accounting for friction.

Simulations are three-dimensional and performed using a grid spacing ∆ = a/5 and a time step ∆t = 10 -3 γ-1 where γ is the imposed shear rate.

Figure 3 displays the computed and theoretical relative trajectories in the plane of shear (z=0) for particles having different non-dimensional roughness ✏ r =h r /a and friction coefficient µ d . The case ✏ r =0 corresponds to smooth particles. Note that in this figure, the vertical coordinate y 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. A very good accordance between simulation and theory is noted, which validates lubrication and contact modelling. As expected, the smooth particle trajectory remains symmetric by virtue of the reversibility of the Stokes equations. On the contrary, rough particles undergo contact since the distance between their surface is smaller than the roughness height at some point. Since contact forces are compressive but not tensile, they eventually result in a fore-aft symmetry breaking and the development of anisotropic microstructures. This asymmetry is clearly visible in the relative trajectories 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. Surprisingly, friction does not modify trajectories, i.e. frictional or frictionless particles result in hardly distinguishable trajectories. Because this paper primarily focuses on rheology -especially normal stress differences -, we plot in Fig. 4 the non-dimensional normal stress differences N * 1 and N * 2 as a function of the pair orientation angle ✓ measured from the flow direction (positive x). The non-dimensional normal stress differences are here defined by

N * 1 =2 (S h xx + S c xx ) -(S h yy + S c yy ) S ∞ (3.1) N * 2 =2 (S h yy + S c yy ) -(S h zz + S c zz ) S ∞ (3.2)
They are made dimensionless using the stresslet S ∞ of a sheared pair of spheres without hydrodynamic interactions S ∞ = 20/3⇡a 3 γ. The results presented are obtained for the case of rough frictionless particles. Here again, the agreement between theory and simulations is satisfactory. The main effect of roughness lies in the extension quadrant (⇡ -✓>⇡/ 2) where particles are further away, thereby reducing hydrodynamic stresses.

In the compression quadrant (⇡ -✓<⇡/ 2), the contact force itself is expected to increase the normal stress due to the additional contact dipole F c ⌦r. As this contact dipole rises, the hydrodynamic dipole decreases, which leads to a very limited effect on the total stress as can be seen in Fig. 4. In [START_REF] Zarraga | Normal stress and diffusion in a dilute suspension of hard spheres undergoing simple shear[END_REF], similar results are obtained but with a stronger contact contribution in the compression region because non-dimensional roughness ✏ r could be as high as 0.5. Surface roughness ✏ r in the range 10 -3 ⇠ 10 -2 is typical of particles used in suspension studies [START_REF] Smart | Measurement of the hydrodynamic surface roughness of noncolloidal spheres[END_REF] [START_REF] Zarraga | Normal stress and diffusion in a dilute suspension of hard spheres undergoing simple shear[END_REF].

The role of friction is investigated in the case of the largest roughness ✏ r =10 -2 and the obtained results are presented in Fig. 5 for µ d =0 (frictionless case) and µ d =0.3. Unlike roughness, friction alters normal stresses only in the compression quadrant and not in the extension region. This is coherent with the fact that friction does not change the pair relative trajectories as pointed out previously. Friction only acts during contact since it implies an additional tangential contribution F c t ⌦ r to the overall stress. The kind of result presented in Fig. 5 can be used to estimate the role of friction on the

Suspension flow simulations : results and discussion

The objective of this study is to investigate the role of roughness and friction on suspension rheology. Numerical simulations of dense suspensions are performed for different volume fractions 0.1  φ  0.45. The number N p of particles is in the range 600⇠1000, depending on the volume fraction. The computational domain is a Couette cell of size L y =20a, L z =20a, and L x ≥ 20a. The length L x is actually modified so as to adjust the number of particles in the desired range and can reach L x =60a for the lowest volume fraction φ=0.1. A shear flow of magnitude γ is imposed by moving upper and lower walls with opposite velocities. Periodic boundary conditions are used in x (velocity direction) and z (vorticity direction). The numerical parameters used are the same as in the previous section, i.e. a grid spacing ∆ = a/5 and a time step ∆t = 10 -3 γ-1 . For simulations including friction, the time step is however decreased to 5 ⇥ 10 -4 γ-1 to allow for a better time-resolution of the rapidly varying tangential forces. Note that unlike inertial collisional regimes (where the collision duration is very small), dense suspensions are rather in a frictional regime and a contact lasts longer, about γ-1 . During most of the contact duration, the tangential force generally changes mildly, except at the beginning of the contact (initial tangential spring stretching) where rapid variations are captured in only af e w( ⇠5) time steps. This explains the time step reduction for frictional simulations. A further decrease does not change results significantly.

All runs are started using random hard-sphere equilibrium configurations obtained from a Monte-Carlo procedure, and the first 50 strains γt are discarded so that the system is allowed to reach a steady state configuration before computing averages. After this transient regime, the simulation is continued for another 100-150 strains. The run is then divided into five statistically independent time-segments of 20-30 strains each. The standard deviation is computed over the means of those five subintervals and then divided by p 5. This represents the statistical variation of the properties and will be given as error bars on the numerical plots. The resulting statistical error is usually less than 1 % for ⌘ s , typically 2⇠4 % for N 2 , and can reach 5⇠10 % for N 1 or Π. Despite the large number of particles, normal stress differences -most notably N 1 -d oe x p e r i e n c e large intermittent fluctuations. As postulated by [START_REF] Singh | Normal stresses and microstructure in bounded sheared suspensions via stokesian dynamics simulations[END_REF], this might be in connection with the formation of large gap-spanning network of particles, although this was not investigated further in this work.

To determine the rheological properties, we recall that the bulk stress in a suspension of rigid particles subjected to a strain rate E ∞ ij is defined as [START_REF] Batchelor | The hydrodynamic interaction of two small freelymoving spheres in a linear flow field[END_REF])

Σ ij = -pδ ij +2⌘E ∞ ij +Σ p ij (4.1)
where Σ p ij is the particle contribution to the bulk stress. For a statistically significant number of particles in the suspension, it is given by

Σ p ij =Σ h ij +Σ c ij (4.2) with Σ h ij = nhS h ij i (4.3) Σ c ij = nhS c ij i (4.4)
where S h ij and S c ij are the hydrodynamic and contact stresslets, respectively, n is the number density of particles and the brackets h•i imply an ensemble average. The splitting of Σ p ij into its hydrodynamic part Σ h ij and contact part Σ c ij will be often used throughout this study.

For a linear shear flow, the relative viscosity ⌘ r = ⌘ s /⌘ of the suspension is

⌘ r =1+ Σ p xy ⌘ γ (4.5)
The normal stress differences are given by

N 1 =Σ p xx -Σ p yy (4.6) N 2 =Σ p yy -Σ p zz (4.7)
and the particle pressure by

Π=- 1 3 (Σ p xx +Σ p yy +Σ p zz ) (4.8)
Since we compute a wall-bounded flow with periodicity imposed only in x and z directions, average quantities depend on vertical position y. As proposed by Yeo & Maxey (2010a), an average volume fraction hφ(y)i can be defined as

hφ(y)i = 1 L x L z h ZZ χ(x)dxdzi (4.9)
where χ(x) is the particle indicator function which is 1 in the particle and 0 elsewhere. Note that hφ(y)i is rather an areal fraction but it is known from stereology theory to be equal to the volume fraction (Delesse principle). Figure 6 presents the average volume fraction hφ(y)i for three bulk fractions φ bulk in the case L y = 20a. Some simulation results by Yeo & Maxey (2010a) at φ bulk =0.4 are also plotted. Local peaks in the wall region indicate the presence of a stable particle layering, which has been attested in other computations (Yeo & Maxey 2010a;[START_REF] Kromkamp | Lattice boltzmann simulation of 2d and 3d non-brownian suspensions in couette flow[END_REF] and in experiments [START_REF] Blanc | Rhéologie et microstructure des suspensions concentrées non-browniennes[END_REF]. Results presented in Fig. 6 are obtained for frictionless particles : friction actually has a weak influence but the effect of friction on wall-induced ordering is not detailed here and will be part of a future work. In this case L y =20a, and for bulk volume fraction φ bulk below 0.5, there is still a flat profile in the core flow which indicates that the suspension is devoid of wall effects there and is therefore expected to behave like an unbounded suspension. On the opposite, for φ bulk =0.5, wall effects are dominant across the whole channel. For low volume fractions, this suggests -as proposed by Yeo & Maxey (2010a) -to split the whole suspension domain D into a core region D C and a wall region D W . For L y =20a, this core region D C exists only if φ bulk < 0.5 and is chosen as

D C =[0,L x ] ⇥ [L y /4, 3L y /4] ⇥ [0,L z ] (4.10)
This choice follows the works of Yeo & Maxey (2010a) and here corresponds to 5 < y/a < 15. The present study intends to highlight the role of friction and roughness for homogeneous suspensions because we believe that those effects must be well understood in that case before adding wall and confinement effects. The reader may wonder why the authors did not consider triperiodic Lees-Edwards boundary conditions to remove wall effects. The reason is that we are also currently studying the effects of confinement, although this will be deferred to a subsequent paper. Thus we have kept the walls in our simulations ; this is moreover expected to be closer to experimental reality. Consequently, all the results presented in this paper are for the core region D C . As shown by Yeo & Maxey (2010a), rheological properties and microstructure in this core region agree well with results obtained in an unbounded domain. For computational time reasons, most results are presented for a channel width L y =20a, which is suitable for bulk volume fractions φ bulk up to 0.4. Some simulations at φ bulk =0.4 and L y =40a did not result in significant changes. Due to the need to keep a homogeneous core, results presented for φ bulk =0.45 are obtained for a wider channel L y =40a. Here again, the results remain unchanged when considering a larger channel width. Therefore, all the results presented in this papereven for the highest fractions -are checked to be independent of wall effects. Such wallinduced layering effects on suspension rheology are complex and out of the scope of the present study. Note that Lees-Edwards boundary conditions are not so difficult to implement but for higher volume fractions (above 0.5), ordering in monodisperse suspensions is nevertheless expected even in fully periodic conditions [START_REF] Kulkarni | Ordering transition and structural evolution under shear in brownian suspensions[END_REF][START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF]). An improvement would rather be to consider polydispersity.

A last remark is that due to confinement effects, the average volume fraction φ core in D C might slightly differ from φ bulk , typically a few percent. We do consider the actual value φ core when presented on the numerical plots. Similarly, the particle-phase velocity profile in the core region is linear but with a slope γcore that can vary by a few percent from the imposed shear rate γ. This is an apparent wall slip, as also observed experimentally in concentrated suspensions. The actual value γcore is considered when computing rheological properties, such as Eq. (4.5).

Table 1 compiles the range of the various parameters investigated in this study.

Effect of roughness

Although this study focuses on friction, we have briefly investigated the role of roughness for frictionless particles with an eye to assessing to which extent roughness could affect suspension rheology. Figure 7 presents the relative viscosity (a) and the normal stress

Volume fraction φ bulk 0.1 ⇠ 0.45 Friction coefficient µ d 0 ⇠ 1 Reduced roughness ✏r 10 -4 ⇠ 10 -2 Contact Péclet Γ1 0 -3 ⇠ 10 -1
Table 1. Range of parameters investigated differences (b) as a function of the reduced roughness ✏ r = h r /a for a suspension of frictionless particles at φ bulk =0.4. Rheological properties are evaluated in the core region D C as discussed above. The total relative viscosity is decomposed into its hydrodynamic and contact parts as suggested by Eq. ( 4.2). Note that the solvent contribution (the 1 in Eq. (4.5)) is lumped into the hydrodynamic contribution.

The relative viscosity decreases monotonically as the roughness increases as seen in Fig. 7(a). Hydrodynamic and contact parts vary in a similar way so that their relative contribution to suspension viscosity remains constant and independent of roughness : ⌘ h ⇡ 0.68⌘ r and ⌘ c ⇡ 0.32⌘ r . Contact forces are balanced by hydrodynamic forces and as the hydrodynamic contribution to viscosity decreases, so does the contact contribution. The hydrodynamic contribution decreases because as asperities grow, the separation distance between particles becomes larger which, in turn, lowers the lubrication dissipation. The decrease in viscosity is however modest and is less than 20 % for a hundred-fold increase in roughness. This weak decrease in viscosity with increasing roughness is also found in SD [START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF][START_REF] Singh | Normal stresses and microstructure in bounded sheared suspensions via stokesian dynamics simulations[END_REF] and FCM simulations (Yeo & Maxey 2010a). Note that in SD and FCM simulations, there is no roughness properly speaking but a repulsive inter-particle force whose range is analogous to an interaction distance that can be viewed as a roughness.

Figure 7(b) shows the normal stress differences normalized by the shear stress ⌧ = ⌘ r ⌘ γ. The first and second normal stress differences are negative and have very similar values N 1 /⌧ ⇡ N 2 /⌧ ⇡-0.15. This is in good agreement with SD simulations in which N 1 /⌧ ⇡ N 2 /⌧ ⇡-0.135 [START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF]). An interesting feature is that N 1 /⌧ and N 2 /⌧ are independent of roughness within statistical uncertainty. This means that, as roughness becomes larger, |N 1 | and |N 2 | slightly decreases. This is in accord with the monolayer SD simulations of [START_REF] Wilson | Shear stress of a monolayer of rough spheres[END_REF] that have showed that |N 1 | increases with increasing roughness for dilute systems (this is also expected from our two-particle results in §3) but that the trend is reversed in concentrated regimes. This effect of roughness on the normal stress differences can not explain the experimental results since |N 1 |/|N 2 |⇡1, in contrast with experiments. Roughness could, at most, be involved in some discrepancies between experiments.

Additional simulations were also performed for another volume fraction (φ bulk =0.2) : the ratio normal stress differences over shear stress has a different value (N 1 /⌧ ⇡ N 2 /⌧ ⇡ -0.03) but remains independent of roughness. The fact that the differences N 1 /⌧ and N 2 /⌧ are independent of roughness could suggest that normal stresses may persist in the limit of vanishing roughness. Nevertheless, we have not investigated this hardly realistic case of perfectly smooth spheres. As addressed by [START_REF] Melrose | The pathological behaviour of sheared hard spheres with hydrodynamic interactions[END_REF], simulations of sheared suspensions of perfectly smooth spheres inevitably result in overlaps (hence, contacts) even for vanishing time steps and high-order integration schemes. As a consequence, we believe that -at least numerically -normal stress differences may exist in dense sheared suspensions with zero roughness. Furthermore, we recall that, in our contact model, the contact Péclet Γ (see Eq. (2.18)) is constant, irrespective of roughness. We have then conducted simulations with constant normal stiffness k n : in this case, Γ now depends on roughness. The same results were obtained (at least, in the tested range 10 -3 . Γ . 10 -1 ). This means that the constant N i /⌧ with respect to roughness is not related to our peculiar choice of roughness-independent Γ.

Full three-dimensional simulations reporting the effect of roughness on normal stress differences are scarce. [START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF] investigated the role of the range of the inter-particle force (as pointed out, this range might be viewed as a kind of roughness) but only on the particle pressure Π. They noticed that Π was rather constant, implying a moderate 25 % increase in Π/⌧ when the force range was increased from 10 -4 to 10 -2 . Actually, our results are quite similar with Π being relatively constant and Π/⌧ increasing by 35 % for a roughness in the range 10 -4 ⇠ 10 -2 , as seen in Fig. 8(a). A similar trend is noted for the three normal stresses Σ p ii with a limited rise in |Σ p ii |/⌧ as roughness height grows (Fig. 8(b)). Interestingly, this variation is similar for the three normal stresses, which results in differences N 1 /⌧ and N 2 /⌧ relatively unchanged by roughness. As already pointed out by [START_REF] Singh | Normal stresses and microstructure in bounded sheared suspensions via stokesian dynamics simulations[END_REF], the fact that simulations predict rather constant normal stresses (or particle pressure) contrasts with the theoretical predictions of [START_REF] Brady | Microstructure of strongly sheared suspensions and its impact on rheology and diffusion[END_REF] that -in dilute suspensions -the normal stresses vanish as ✏ 0.22 r .

Effect of friction on viscosity

We now turn to simulations with frictional particles having a roughness maintained to ✏ r =5.10 -3 throughout the forthcoming simulations. This value is typical for particles used in suspension studies [START_REF] Smart | Measurement of the hydrodynamic surface roughness of noncolloidal 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]. Figure 9 presents the φ-dependence of the suspensions relative viscosity ⌘ r for different friction coefficients µ d . A number of experimental results by [START_REF] Zarraga | The characterization of the total stress of concentrated suspensions of noncolloidal spheres in newtonian fluids[END_REF] and [START_REF] Dbouk | Normal stresses in concentrated non-brownian suspensions[END_REF] are also shown. For frictionless particles (µ d =0), the agreement is moderate with predicted viscosities systematically lower than experimental results. This was already noticed in previous simulations (Yeo & Maxey 2010a;[START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF]. When friction is taken into account, a substantial increase in viscosity is noted and the agreement with experiments is significantly improved. The effect of friction is quite limited at low volume fractions φ . 0.3 because there are less particles in contact. For dense suspensions however friction gives rise to a notable increase in viscosity. This result suggests that the use of frictionless contact model could explain the systematic underprediction of viscosity reported in previous simulations. Indeed, it has been stressed in the previous section that roughness can only result in limited variations in the viscosity and can therefore not explain experimental results per se. [START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF] studied a suspension at φ=0.4 and noted a limited 10 % rise in ⌘ r for µ d =0.5. For present simulations at φ=0.4, the effect is stronger with almost a two-fold increase. This quantitative difference is likely due to the friction model. The model used in [START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF] was deliberately kept simple in order to explore the effect of friction. It was assumed that particles were sliding (no roll-slip behaviour) and that the frictional force was always opposing the bulk shearing motion.

The role of friction on viscosity can be better understood by studying the hydrodynamic and contact contribution to viscosity, as already done in §4.1. Figure 10 presents the hydrodynamic and contact contribution to viscosity as a function of volume fraction for the frictionless case (µ d =0) and a frictional case µ d =0.5. The high-frequency dynamic viscosity ⌘ ∞ is also plotted. It was obtained by simulations on frozen equilibrium configurations of hard smooth spheres. Note that this viscosity is purely hydrodynamic in origin and corresponds to the viscous contribution to the stress for an equilibrium microstructure, i.e. not affected by the flow. A first significant result is that the hydrodynamic contribution to viscosity is not sensitive to friction : it can be seen in Fig. 10 that frictionless and frictional hydrodynamic contributions collapse onto the same curve. The only difference between the frictionless and frictional cases lies in the contact contribution which significantly increases with friction. Another interesting result is that the hydrodynamic contribution to viscosity ⌘ h remains very close to the high-frequency vis- cosity, at least for moderate volume fractions. This agreement might be fortuitous since it is clear from Fig. 7 that ⌘ h depends on roughness. However this dependence is very limited for usual roughness with a mere 5 % decrease between ✏ r =10 -3 and 10 -2 . As a reasonable assumption, it can be argued that ⌘ h ⇡ ⌘ ∞ . This result hints at a very limited impact of microstructure anisotropy on hydrodynamic viscosity ⌘ h . As a consequence, this means that the well-known excess viscosity ∆⌘ = ⌘ r -⌘ ∞ can be assimilated to the contact viscosity, i.e. ∆⌘ ⇡ ⌘ c . (2002) but the qualitative behaviour is similar. In particular, the effect of friction is more pronounced on N 2 . Because contacts predominantly take place in the shear plane (x,y), friction gives rise to larger normal stresses Σ xx and Σ yy , resulting in a marked increase in N 2 =Σ yy -Σ zz .

The trends observed are an interesting result since experimental studies show that the second normal stress difference is always significantly larger than the first normal difference, i.e. N 2 /N 1 > 1. Figure 12 plots our results in terms of N 2 /N 1 as a function of volume fraction and friction coefficient. In accordance with simulations from the literature, frictionless particles (µ d =0) are reported to give N 2 /N 1 ⇡ 1. This ratio is relatively independent of the volume fraction. For low φ however, N 2 /N 1 is slightly smaller than unity in qualitative agreement with the theoretical results of [START_REF] Zarraga | Normal stress and diffusion in a dilute suspension of hard spheres undergoing simple shear[END_REF] Figure 13 displays the hydrodynamic ("H") and contact ("C") contribution to the normal stress differences N 1 /⌧ (a) and N 2 /⌧ (b) for frictionless (µ d =0) and frictional case (µ d =0.5). For the first normal stress difference in Fig. 13(a), it is found that the hydrodynamic contribution is independent of friction while the contact contribution is less negative in the frictional case. The noted decrease in |N 1 | is therefore solely due to a reduction in the contact contribution. This predominant role of contact contribution in the case of friction has also been noticed previously for viscosity. For moderate volume fractions, the hydrodynamic stress represents the major contribution to N 1 while for denser regimes, contact and hydrodynamic stresses contribute in a similar way. For the second normal stress difference in Fig. 13(b), similar conclusions hold : friction weakly affects the hydrodynamic contribution while it considerably increases the contact contribution. Figure 13(b) clearly shows that the hydrodynamic contribution to N 2 is always small. Therefore normal stress differences do not have the same origin : for N 1 ,i ti s mostly hydrodynamic -at least for low and intermediate volume fractions -while for N 2 it comes entirely from contacts. This is anticipated since contacts preferentially occur in the shear plane and in the compression quadrant, so that |Σ c zz |⌧|Σ c yy |⇡|Σ c xx |.

Effect of friction on normal stresses

The previous section was entirely dedicated to normal stress differences because of their importance in rheology. From a more fundamental viewpoint however, the study of normal stresses is essential. Furthermore, there exist only very scarce experimental data on normal stresses in suspensions [START_REF] Zarraga | The characterization of the total stress of concentrated suspensions of noncolloidal spheres in newtonian fluids[END_REF][START_REF] Dbouk | Normal stresses in concentrated non-brownian suspensions[END_REF][START_REF] Couturier | Suspensions in a tilted trough: second normal stress difference[END_REF]Boyer et al. 2011b;[START_REF] Garland | Normal stress measurements in sheared non-brownian suspensions[END_REF]. This section considers the role of friction on normal stresses and comparisons with experiments will be provided in the forthcoming sections.

Figure 14 presents the obtained results on normal stresses Σ p ii and particle pressure Π for frictionless particles (µ d =0) and frictional particles with µ d =0.5. All those quantities are made non-dimensional using the shear stress ⌧ . The effect of friction is similar for all the normal stresses, i.e. friction results in a substantial increase in the stresses for φ & 0.3. Friction implies a doubling of Π/⌧ , which basically corresponds to a three-fold increase in Π. It will be shown hereinafter ( §4.7) that the contact particle pressure Π c also rises significantly because of friction. Since Π c /h F c .ri is a measure of normal forces, this shows that friction also involves a much higher level of normal forces.

Microstructure

In Fig. 15, the pair-distribution function g(r) is plotted as a function of the angle ✓ for µ d =0 and µ d =0.5. Given a particle at the origin, the pair-distribution function g(r) describes the probability density of finding a second particle at distance r, normalized by the particle number density in the suspension. It is here determined in the usual manner by discretizing the area around each particle and counting the number of neighbouring particles within the area. In spherical coordinates, r =( r, ✓, )w h e r er is the radial distance, ✓ the azimuthal angle measured from the flow direction (positive x) and the polar angle measured from the vorticity direction (positive z). Figure 15(a) shows the pair-distribution computed in the plane of shear and for r<2.1a,i . e .g(✓)=g(r/a < 2.1,✓,⇡/2) while Fig. 15(b) presents the actual contact pair-distribution, i.e. g c (✓)= g(r/a < 2+✏ r ,✓,⇡/2) with ✏ r =5.10 -3 . For both cases, this pair-distribution function is calculated in the core region. Considering the mirror symmetry, only 0 <✓<⇡is shown in Fig. 15.

Both pair-distribution functions have an asymmetric shape with a high probability to encounter another particle in the compression region (⇡ -✓<⇡ / 2) and a marked depletion zone around the extension axis ⇡ -✓ ⇡ 3⇡/4. This depletion zone is much more apparent for particles at contact (Fig. 15(b)). This fore-aft asymmetry is well-known to be responsible for the non-Newtonian behaviour of dense suspensions and has been observed in experiments [START_REF] Blanc | Experimental signature of the pair trajectories of rough spheres in the shear-induced microstructure in noncolloidal suspensions[END_REF][START_REF] Parsi | Fore-and-aft asymmetry in a concentrated suspension of solid spheres[END_REF]. Our results in Fig. 15(a) show an excellent agreement with FCM simulations by Yeo & Maxey (2010a) for frictionless particles (this paper does not report the contact pair-distribution, so that the comparison is not shown in Fig. 15(b).)

The effect of friction on the microstructure is clear in the compression region, especially for particles at contact. Friction involves less particles in the direction of velocity ⇡-✓ ⇡ 0 but more particles in the region ⇡/4 <⇡-✓<⇡/ 2. This corresponds to a smaller density of normal force that would give rise to a negative N 1 and a larger density of normal force that would give rise to a positive N 1 . As a result, the values of N 1 become less negative with friction. The opposite effect occurs for N 2 . Overall, this angular redistribution of particles leads to a rather flat g(✓) profile in the compression region, which is consistent with small N c 1 as noted in Fig. 13(a). This effect of friction on microstructure is in qualitative agreement with the simulations from [START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF].

Dependence on the friction coefficient

We have investigated the dependence of viscosity and normal stress differences on the friction coefficient for a suspension at φ bulk =0.4, as shown in Fig. 16. The relative viscosity ⌘ r increases monotonically with friction coefficient, although a slope break is noted for µ d ⇡ 0.5. On the contrary, the non-dimensional normal stress differences N 1 /⌧ and N 2 /⌧ vary notably at low friction coefficients but saturate at some point, for µ d ⇡ 0.5. This means that the ratio N 2 /N 1 rapidly changes for low friction coefficients but then becomes independent of µ d for µ d & 0.5. Note that even in that regime, N 1 and N 2 keep changing -albeit moderately -since ⌧ = ⌘ r ⌘ γ increases with µ d .

It is apparent from Fig. 17

N 1 /⌧ =N h 1 /⌧ +N c 1 /⌧ ⇡ N h 1 /⌧
, which is independent of µ d as already seen in Fig. 13(a). Figure 17(b) shows that the contact normal stresses Σ c ii /⌧ are rather constant for sufficiently high friction coefficients. This, in turn, implies that N c 2 /⌧ =(Σ c yy -Σ c zz )/⌧ is independent of friction coefficient and so does N 2 /⌧ =N h 2 /⌧ +N c 2 /⌧ . This splitting between contact and hydrodynamic contributions reveals more clearly the important role of contacts in dense suspensions. Results in Fig. 17 can be used to provide the relative contribution of contacts at φ=0.4 as a function of µ d . For µ d =0.5 for instance, we calculate Σ c xx /Σ p xx ⇡ 0.83, Σ c yy /Σ p yy ⇡ 0.91, and Σ c zz /Σ p zz ⇡ 0.82. It means that in this case, contacts are responsible for 80% ⇠ 90 % of the total particle stress. This predominant role of contact in dense frictional suspensions will be confirmed in the following sections as well.

Our results show that the normal stress differences are very sensitive to friction, especially for µ d . 0.5, which is the typical range of friction coefficient for usual materials. This means that friction could be at the origin of some discrepancies between experiments since particles used in suspension studies are made of different material and hence have different friction coefficients. It is also known that the dynamic friction coefficient µ d may depend upon normal contact force, slip velocity, roughness, temperature, etc. which may also contribute to the overall discrepancy.

Continuum models

Particle-resolved simulations offer an inestimable insight into the detailed physics of suspensions. However they are too computationally expensive to apply for large-scale industrial systems with complex geometries. Therefore there have been attempts to develop continuum models where the suspension is assimilated to an homogeneous fluid with some constitutive laws for the particle stress Σ p . A successful approach for shear flows was proposed by [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF] -referred to as the suspension balance model (SBM) -, in which where η p n is the normal stress viscosity and λ p 2 and λ p 3 are anisotropy parameters λ p 2 = Σ p yy /Σ p xx and λ p 3 =Σ p zz /Σ p xx .S u p e r s c r i p tp is here added to recall that those quantities are computed using the particle stress Σ p .

Σ p = -⌘ γQ +2⌘⌘ r E (4.11)
Due to the lack of experimental data, [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF] proposed constant anisotropy parameters λ p 2 ⇡ 0.8 and λ p 3 ⇡ 0.5. The SBM rests on the idea that the overall behaviour of the suspension is driven by the particle stress Σ p . However the exact nature of the particle stress to be included in the model is still controversial. In particular, [START_REF] Lhuillier | Migration of rigid particles in non-brownian viscous suspensions[END_REF] and [START_REF] Nott | The suspension balance model revisited[END_REF] demonstrated that the contact contribution to the particle stress, namely Σ c , should be considered in the SBM model instead of Σ p . This possible flaw in the SBM might not impair the good prediction capabilities of this model, for two reasons. First, the anisotropy parameters have been chosen by [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF] to match experimental results on migration, so that their exact nature is unimportant ; and second, we have shown that in dense frictional suspensions, contact stress is the prevailing stress.

In this section, we present simulation results on the SBM parameters, i.e. normal viscosity and anisotropy parameters. The main purpose is to study the role of friction but also to analyse whether the particle stress or the contact stress should be considered in the model. Direct measurements of SBM parameters are very scarce [START_REF] Dbouk | Normal stresses in concentrated non-brownian suspensions[END_REF] because the three normal stresses must be measured simultaneously. Furthermore, strong theoretical arguments are still missing to determine whether experiments actually measure the particle stress Σ p (as defined by Batchelor, see Eq. (4.1)), or the contact contribution Σ c to particle stress, or another stress.

Our results for anisotropy parameters λ p 2 =Σ p yy /Σ p xx and λ p 3 =Σ p zz /Σ p xx are presented in Fig. 18 for frictionless (µ d =0) as well as frictional particles µ d =0.5. For frictionless particles, our simulations are close to Yeo & Maxey (2010a) and a moderate linear increase in the anisotropy parameters is found. Friction has a modest effect and induces a slight rise in λ p 2 and λ p 3 . The obtained values remain close to the ones proposed in the SBM model. However, computations suggest that λ p i may have a φ-dependence. Anisotropy parameter λ p 2 is experimentally noted to increase linearly with φ but with values slightly above the predictions. Conversely, experimental λ p 3 is found to be relatively constant and close to 0.5, which is the value introduced by [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF] in their model to represent the lack of migration in torsional flows.

Figure 19 presents similar plots for the contact anisotropy parameters, defined as λ c

2 =Σ c yy /Σ c xx and λ c 3 =Σ c zz /Σ c xx . Those parameters are expected to replace the usual λ p i if the contact stress were to be considered instead of the particle stress. On the figure are also shown the experimental results from [START_REF] Dbouk | Normal stresses in concentrated non-brownian suspensions[END_REF]. They are the same as on Fig. 18 since it is not clear whether they are representative of Σ p or Σ c . Interestingly, the agreement with experiments is clearly improved for both parameters. In particular, there is an excellent match for λ c 3 which is relatively constant and close to 0.5 for all volume fractions. λ c 2 is closer to experiments as well, although the φ-dependence seems to differ since λ c 2 decreases with φ. Just as for λ p i , the effect of friction on λ c i is weak and the values remain close to the ones chosen by [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF]. The fact that λ c i parameters match experiments better than λ p i is puzzling. We however believe it is too premature to conclude whether experiments actually measure the contact particle stress.

The computed normal viscosity η p n and contact normal viscosity η c n are shown in Fig. 20 together with experiments from [START_REF] Dbouk | Normal stresses in concentrated non-brownian suspensions[END_REF]. Experimental data are here again duplicated in both Fig. 20(a) and Fig. 20(b) due to the uncertainty on their physical meaning. Simulation results from Yeo & Maxey (2010a) are also displayed (available only for η p n and frictionless particles) and a good agreement is noted. Unlike anisotropy parameters, normal viscosities η p n and η c n strongly increase with friction -this was already the case for the shear viscosity η r . Accounting for friction (µ d =0.5) allows us to match closely the experimental data. This good agreement with experiments is obtained for both usual normal viscosity η p n and contact normal viscosity η c n . This is an additional evidence that friction must be considered for relevant simulations of dense suspensions.The ratio η c n /η p n can be calculated from those results and is found to be relatively constant for φ & 0.3w i t hη c n /η p n ⇡ 0.8. This again confirms that contact stress prevails in dense suspensions. [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF] proposed an empirical law for the normal viscosity given as In a recent work, Boyer et al. (2011a) intended to unify suspension and granular rheology. They showed that dense confined suspensions have a viscoplastic behaviour similar to dry granular media that can be described using a single viscous parameter I v . This viscous number is defined as I v = η γ/P p where P p is the confining pressure. This parameter is the ratio between a flow strain time γ-1 and a viscous time η/P p . It is analogous to the inertial number I used in granular rheology I = a γp ρ p /P p where ρ p is the particle density. Using a pressure-imposed cell, they found that the effective friction coefficient µ and the volume fraction φ were uniquely defined by I v ,i . e .µ = µ(I v ) and φ = φ(I v ) whatever the confining pressure and for two different types of particles. However, only dense suspensions (φ>0.4) and consequently very small viscous numbers I v were investigated.

η p n = K n ✓ φ/φ m 1 -φ/φ m ◆ 2 (4.13)
Here, we intend to compare our simulations -for frictional as well as frictionless particles -to their results. Our goal is to gauge the effect of friction as well as to check whether their experimental observations still hold for large I v regimes. In their measurements, the confining pressure was prescribed by applying a normal force to the suspension using a porous plate. In order to parallel experiments, this suggests to consider the normal stress Σ p yy as the confining pressure P p . To be more specific, we actually have investigated Σ c yy as well as Σ p yy and have found that results are closer to Boyer's experiments if we consider Σ c yy . This could again suggest that actual measurements determine Σ c rather than Σ p . (2011a). For low volume fractions, a deviation between Boyer's results and our simulations is clearly noted and grows with increasing I v . Yet, let us recall that their experimental results are obtained for dense suspensions with I v in the range 10 -6 ⇠10 -1 while our results are restricted to much higher viscous numbers 10 -1 ⇠10 2 . It is then questionable whether their experimental law can be extended to very large I v . Simulations show that φ m and n should be slightly adjusted in the semi-dilute case with φ m ⇡ 0.64 and n ⇡ 0.4. Note that Eq. (4.14) is obtained considering the divergence of viscosity near maximum packing fraction and there is little chance that very dilute suspensions can be modelled by laws obtained through the behaviour of suspensions near jamming. The fact that our fitted value φ m ⇡ 0.64 is close to the random-close packing fraction is probably fortuitous.

The second constitutive law addresses the effective friction coefficient µ (not to be mistaken with material dynamic friction coefficient µ d ) defined as µ = τ /P p (4.15) where τ is the shear stress τ = η r η γ. Here again, P p is assimilated to Σ c yy in our simulations. Boyer et al. (2011a) proposed to model the friction law as the sum of the hydrodynamic µ h and contact µ c contribution by

µ(I v )=µ 1 + µ 2 -µ 1 1+I 0 /I v | {z } µ c + I v + 5 2 φ m I 1/2 v | {z } µ h (4.16)
Using granular rheology results, they prescribed µ 1 =0.32, µ 2 =0.7 and I 0 =0.005. Figure 22 presents our simulations for the friction coefficient µ. Similarly, there is a reasonable collapse of all our results onto a single curve. This is an interesting result since computations span very different roughness or friction coefficient. Like the φ(I v )c u r v e , the agreement with Boyer's results is good for small I v but a deviation is noted for dilute suspensions (large I v ). Again, this might be expected since their fit was obtained for I v smaller than 10 -1 . The dashed line represents an improved model for dilute systems which is described hereafter.

In order to improve the friction law in dilute regimes, we first study separately the hydrodynamic and contact contribution to µ, defined as

µ h = η h η γ Σ c yy (4.17) µ c = η c η γ Σ c yy (4.18)
This is straightforward in computations but complex from an experimental viewpoint. Those two contributions can then be compared to the expressions of µ h and µ c proposed in Eq. (4.16). Results are presented in Fig. 23(a) and Fig. 23(b) for the hydrodynamic and contact parts, respectively. As seen in Fig. 23(b), the contact contribution to µ shows a reasonable agreement with the model Eq. (4.16). Note that because I v is rather large (I v I 0 ), Eq. ( 4.16) boils down to µ c ⇡ µ 2 =0 .7. This good agreement is a striking result because the value of µ 2 is taken from experiments reported in dry granular media. This means that the contact contribution in dry granular media and dilute suspensions is similar. Note that the contact contribution is moderately dependent on friction -frictional particles would give a mere µ 2 ⇡ 0.8 -, at least for this range of I v . Simulation results are however in a too dilute regime to provide a refined modelling of µ c for a wide range of I v , and only the high-I v limit µ 2 can be inferred here.

For moderately dense suspensions, the most important contribution comes from hydrodynamics as seen in Fig. 23(a) (note the difference in scale with Fig. 23(b)). Actually, Eq. ( 4.16) predicts µ c ' µ h for I v ⇡ 10 -1 . Therefore for dilute systems it is of prime importance to first improve the hydrodynamic contribution to µ. We have previously mentioned that η h ⇡ η ∞ seems a reasonable assumption in the semi-dilute regime. The relation η ∞ (φ) can be readily inferred from simulations unequivocally. Since it is obtained from a hard-sphere equilibrium non-touching configuration, it does not involve contacts and is neither related to a distorted microstructure. Past computations [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF] show that the high-frequency viscosity can be correctly modelled using a Krieger-Dougherty law

η ∞ =(1- φ φ 0 m ) -[η]φ 0 m (4.19)
with [η] ⇡ 2.4 and φ 0 m ⇡ 0.68. The approximation η h ⇡ η ∞ is used, with η ∞ given by Eq. (4.19), and further combined with Eq. (4.14). This finally gives characteristics (e.g. roughness or friction coefficient). Although the overall collapse is reasonable, a closer look reveals some slight dispersions, meaning that an effect of particle characteristics -although limited -does exist. Among the investigated parameters (stiffness, roughness, friction), friction is found to play the most important role. Figure 24 reconsiders the φ(I v ) results already presented in Fig. 21, but only retains the frictionless and frictional (µ d =0.5) particles. The aforementioned deviations are here more distinct and this suggests that -rigorously speaking -a unique master curve is not likely to exist due to friction. This effect was not noted in Boyer's measurements : first, because of the experimental difficulty to vary the friction coefficient significantly ; and second, because they considered very small I v numbers where the role of friction might be different. Interestingly, Fig. 24 shows that frictional simulations are closer to the experimental correlation, suggesting that friction is operative in actual suspensions. This is especially true for high volume fractions (small I v ). A separate fit on our µ d =0 and µ d =0.5 results gives φ m ⇡ 0.69 and φ m ⇡ 0.62, respectively, i.e. the frictional φ m is in better agreement with the experimental φ m ⇡ 0.585.

µ h = I v ✓ 1 - φ m φ 0 m 1 1+I n v ◆ -[η]φ 0 m
Let us summarize this section. The overlapping range (in volume fraction) between Boyer's experiments and our simulations is very small which hinders any further validations in dense regimes. Yet, one of our primary purpose was to check the existence of the master curves µ = µ(I v ) and φ = φ(I v ) and our results show that this reasonably works. The experimental laws seem to be less accurate in dilute regimes and we have therefore proposed a modified Boyer-like correlation to fit our dilute results. Also, it seems that friction is responsible for some slight deviations from a unique φ(I v ) curve ; as a consequence, master curves might not exist, rigorously speaking. Frictional simulations are much closer to experiments, substantiating a role of friction in actual suspensions.

Conclusions

In this paper, we have presented some three-dimensional numerical simulations of concentrated suspensions in a Couette flow using a fictitious domain method. Our numerical method includes a lubrication model as well as a contact model between particles similar to DEM (Discrete Element Method). This contact model assumes a Hertz law and includes roughness and friction. Some validations on a pair of particles are presented in order to demonstrate that our numerical model is adequate to deal with rough frictional particles.

This work particularly focuses on the effect of friction and its role on rheological properties, especially the relative viscosity η r and the normal stress differences N 1 and N 2 . A relatively new and significant result is that friction has a profound impact on suspension rheology. Friction increases η r and |N 2 | and decreases |N 1 |, in better agreement with available experiments. In particular, the N 2 /N 1 ratio is typically about 3⇠4, which corroborates experiments from [START_REF] Zarraga | The characterization of the total stress of concentrated suspensions of noncolloidal spheres in newtonian fluids[END_REF]. Yet, we have not found quasi-zero or positive N 1 as attested by some other measurements [START_REF] Couturier | Suspensions in a tilted trough: second normal stress difference[END_REF][START_REF] Dbouk | Normal stresses in concentrated non-brownian suspensions[END_REF]. Interestingly, friction is shown to act mostly through the contact stress since the hydrodynamic stress is found to remain unaffected by friction. Slight modifications of the microstructure in the compression region are likely to explain the noted effects on N 1 and N 2 . The systematic splitting between hydrodynamic and contact contributions to the stress reveals that N 1 and N 2 have not the same origin since N 1 is mostly hydrodynamic while N 2 arises from contact. This also shows that the contact stress is the prevailing stress in dense suspensions ; this is particularly true for the normal stresses. The dependence on the friction coefficient µ d is investigated as well and is found to be significant in the range 0⇠0.4.

The link with continuum models is also studied in the frame of the suspension balance model. Some important parameters of the model, such as the normal viscosity or the anisotropy parameters, are computed and are found to be in better agreement with experiments when friction is accounted for. Similarly, our results are compared with recent µ(I v ) global rheology and appear to collapse onto a single curve, as expected from the works of Boyer et al. (2011a). Here again, the agreement with experiments is improved for frictional particles, meaning that friction is operative in actual suspensions.

Figure 1 .

 1 Figure 1. Sketch of roughness modelling : δ is the actual separation distance, hr the roughness height, and contact occurs when δ 0 = δ -hr  0.

Figure 2 .

 2 Figure 2. Sketch of a pair of particles interacting in a shear flow.

Figure 3 .

 3 Figure 3. Relative trajectories for a pair of particles having non-dimensional roughness ✏r=0 ;1 0 -3 ;1 0 -2 : theory (solid lines) and computations for µ d =0 (•);µ d =0.3 (⇤).

Figure 4 .

 4 Figure 4. Non-dimensional normal stress differences N * 1 and N * 2 as a function of pair orientation angle ✓ for three roughness ✏r=0 ; 10 -3 ;10 -2 : computations (open symbols) and theory (lines). Particles are frictionless.

  rheology of dilute suspensions by integrating over all possible pair trajectories. This was done by Wilson & Davis (2000) who found that friction induces a slight rise in |N 1 | and |N 2 | with negative values of N 1 and N 2 , and |N 1 | > |N 2 |. However, this may hold only for dilute suspensions since our simulations for dense systems will show a reverse trend, i.e. friction leads to a decrease in |N 1 | and increase in |N 2 | and |N 1 | < |N 2 |.T h i s was already pointed out by Wilson & Davis (2002) using monolayer SD simulations. Let us conclude by stressing again that the previous single-trajectory simulations should be only viewed as a validation. As pointed out, statistical consideration is required to infer relevant information on the rheology of dilute suspensions.

Figure 5 .

 5 Figure 5. Non-dimensional normal stress differences N * 1 and N * 2 as a function of pair orientation angle ✓ for frictionless (µ d =0) and frictional particles (µ d =0.3) : computations (open symbols) and theory (lines).

Figure 6 .

 6 Figure 6. Local volume fraction hφ(y)i in the channel width for Ly =2 0 a and three bulk fractions φ bulk =( 0 .3, 0.4, 0.5). Open symbols are computations by Yeo & Maxey (2010a). Particles are frictionless.

Figure 7 .

 7 Figure 7. Relative viscosity (a) and normal stress differences (b) as a function of roughness for φ bulk =0.4. Total ("T") relative viscosity is decomposed into hydrodynamics ("H") and contact ("C") parts.

Figure 8 .

 8 Figure 8. Non-dimensional particle pressure Π/⌧ (a) and normal stresses -Σ p ii /⌧ (b) as a function of roughness for φ bulk =0.4.

Figure 9 .

 9 Figure 9. Relative viscosity ⌘r as a function of volume fraction φ for different friction coefficients µ d . Also shown are experiments from Zarraga et al. (2000) and Dbouk et al. (2013). Non-dimensional roughness is ✏r=5.10 -3 .

Figure 10 .

 10 Figure 10. Hydrodynamic ("H") and contact ("C") contribution to relative viscosity ⌘r as a function of volume fraction φ for friction coefficient µ d =0 (white) and µ d =0.5 (grey). Non-dimensional roughness is ✏r=5.10 -3 . The high-frequency viscosity ⌘∞ is also shown.

Figure 11 . 4 . 3 .

 1143 Figure 11. Non-dimensional N1/⌧ (a) and N2/⌧ (b) as a function of volume fraction φ for different friction coefficients µ d . Non-dimensional roughness is ✏r=5.10 -3 .

  Figure11presents the results obtained on the normal stress differences (non-dimensional using the shear stress ⌧ = ⌘ r ⌘ γ) for different friction coefficients µ d . The effect of friction on normal stress differences is significant with a reduction in |N 1 | and a rise in |N 2 |.This effect is in accordance with the scarce previous studies[START_REF] Wilson | Shear stress of a monolayer of rough spheres[END_REF][START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF]. As already discussed, the only three-dimensional simulations with friction are the works from[START_REF] Sierou | Rheology and microstructure in concentrated noncolloidal suspensions[END_REF] who reported for φ =0 .4 and µ d =0.5, a 40 % reduction in |N 1 | and a 75 % increase in |N 2 |. Our simulations show a slightly different effect with a 30 % reduction in |N 1 | as well as a three-fold increase in |N 2 |.T h e differences are likewise possibly due to a simpler friction model used by Sierou & Brady(2002) but the qualitative behaviour is similar. In particular, the effect of friction is more pronounced on N 2 . Because contacts predominantly take place in the shear plane (x,y), friction gives rise to larger normal stresses Σ xx and Σ yy , resulting in a marked increase in N 2 =Σ yy -Σ zz .The trends observed are an interesting result since experimental studies show that the second normal stress difference is always significantly larger than the first normal difference, i.e. N 2 /N 1 > 1. Figure12plots our results in terms of N 2 /N 1 as a function of volume fraction and friction coefficient. In accordance with simulations from the literature, frictionless particles (µ d =0) are reported to give N 2 /N 1 ⇡ 1. This ratio is relatively independent of the volume fraction. For low φ however, N 2 /N 1 is slightly smaller than unity in qualitative agreement with the theoretical results of[START_REF] Zarraga | Normal stress and diffusion in a dilute suspension of hard spheres undergoing simple shear[END_REF] in dilute regimes. Friction results in a larger N 2 /N 1 ratio that grows with friction coefficient and volume fraction. For dense regimes φ & 0.4 and large friction coefficients µ d & 0.3, N 2 /N 1 is in the range 3⇠4, which matches quite well the experiments from Zarraga et al. (2000) who found N 2 /N 1 ⇡ 3.6. Unlike their results however, our simulations suggest that N 2 /N 1 has a φ-dependence which is well-marked in the range 0.2⇠0.4. Below φ ⇡ 0.2, the effect of friction is very limited. Even though our simulations with frictional particles are closer to some experimental results from the literature, they are not likely to explain all of them. Although the values of |N 1 | in the core flow are

Figure 12 .

 12 Figure 12. Ratio N2/N1 as a function of volume fraction φ for different friction coefficients µ d . Non-dimensional roughness is ✏r=5.10 -3 .

Figure 13 .

 13 Figure 13. Hydrodynamic ("H") and contact ("C") contribution to N1/⌧ and N2/⌧ as a function of volume fraction φ for friction coefficient µ d =0 (white) and µ d =0.5 (grey). Non-dimensional roughness is ✏r=5.10 -3 .

Figure 14 .

 14 Figure 14. Normal stress -Σ p xx (a) ; -Σ p yy (b) ; -Σ p zz (c) and particle pressure Π as a function of volume fraction φ for friction coefficient µ d =0 (•)a n dµ d =0.5 (⇤). Non-dimensional roughness is ✏r=5.10 -3 .

Figure 15 .

 15 Figure 15. Angular dependence of the pair-distribution function g(✓) for a system at φ bulk =0.4 and µ d =0 (•);µ d =0.5 (⇤). g(✓)i sc o m p u t e df o rr/a < 2.1 (a) and for contact r/a < 2+✏r (b). Non-dimensional roughness is ✏r=5.10 -3 . Also shown are results from Yeo & Maxey (2010a).

  (a) that the normal stresses Σ p ii are actually more affected by friction than the shear stress ⌧ because the ratio |Σ p ii |/⌧ keeps increasing with µ d . It typically doubles for µ d between 0 and 0.8. Figure 17(b) presents the contact contribution Σ c ii to the normal stress Σ p ii . The values of contact normal stresses Σ c xx and Σ c yy come closer as friction coefficient grows and for µ d & 0.5, we have Σ c xx ⇡ Σ c yy . As a consequence, N c 1 ! 0 as friction increases and this results in

Figure 16 .

 16 Figure 16. Relative viscosity ⌘r (a) and non-dimensional normal stress differences N1/⌧ and N2/⌧ (b) as a function of µ d . Simulations are performed at φ bulk =0.4. Non-dimensional roughness is ✏r=5.10 -3 .

Figure 17 .

 17 Figure 17. Normal stresses -Σ p ii (a) and contact normal stresses -Σ c ii (b) as a function of µ d . Simulations are performed at φ bulk =0.4. Non-dimensional roughness is ✏r=5.10 -3 .

Figure 18 .

 18 Figure 18. Anisotropy parameters λ p 2 (a) and λ p 3 (b) as a function of volume fraction φ for friction coefficient µ d =0 and µ d =0.5. Non-dimensional roughness is ✏r=5.10 -3 . Also shown are simulations from Yeo & Maxey (2010a)( µ d =0) and experiments from Dbouk et al. (2013). The dotted lines are the values taken in the SBM model (Morris & Boulay 1999).

Figure 19 .

 19 Figure 19. Contact anisotropy parameters λ c 2 (a) and λ c 3 (b) as a function of volume fraction φ for friction coefficient µ d =0 and µ d =0.5. Non-dimensional roughness is ✏r=5.10 -3 . Experiments are from Dbouk et al. (2013) (same as on Fig. 18). The dotted lines are the values taken in the SBM model (Morris & Boulay 1999).

Figure 20 .

 20 Figure 20. Normal viscosity ⌘ p n (a) and contact normal viscosity ⌘ c n (b) as a function of volume fraction φ for friction coefficient µ d =0 and µ d =0.5. Non-dimensional roughness is ✏r=5.10 -3 . Also shown in (a) are simulations from Yeo & Maxey (2010a)( µ d =0). Experiments are from Dbouk et al. (2013) and are the same in (a) and (b).

Figure 21 .

 21 Figure 21. Volume fraction φ as a function of Iv for frictionless (•) and frictional (⇤)s i mulations. Solid line is Boyer's results as represented by Eq. (4.14) with φm=0.585 and n=0.5. Dashed line is for φm=0.64 and n=0.4.

Figure 22 .

 22 Figure 22. Friction coefficient µ = ⌧/Σ c yy as a function of Iv for frictionless (•) and frictional suspension (⇤) simulations. Solid line is Boyer's results as represented by Eq. (4.16) with φm=0.585. Dashed line is Eq. (4.20) with φm=0.64 and n=0.4.

Figure 23 .

 23 Figure 23. Hydrodynamic friction coefficient µ h (a) and contact friction coefficient µ c (b) as a function of Iv for frictionless (•) and frictional (⇤)s i m u l a t i o n s .S o l i dl i n ei sE q .( 4 . 1 6 )w i t h φm=0.585. Dashed line is Eq. (4.20) with φm=0.64 and n=0.4.

0 m

 0 I v is large, and noting that φ ⇡ φ m , this expression becomesµ h = I v +[η]φ m I 1-n v (4.21) which essentially gives back Boyer's model Eq. (4.16) with a slight change in the exponent, i.e. 1-n=0.6 instead of 0.5 in the original model. The previous results Fig. 21 to Fig. 23 display all of our available simulations in order to illustrate that results do collapse onto a single master curve regardless of particle

Figure 24 .

 24 Figure 24. Volume fraction φ as a function of Iv for frictionless (•) and frictional µ d =0.5 (⇤) particles. Non-dimensional roughness is ✏r=5.10 -3 . Solid line is Boyer's experimental fit Eq. (4.14) with φm=0.585 and n=0.5.

  . Since it decreases the positive values of N * 1 (✓) and N * 2 (✓) in the extension quadrant, roughness contributes to overall negative N * 1 and N * 2 . Likewise, |N * 1 | and |N * 2 | are expected to increase with surface roughness in agreement with Wilson & Davis (2002); Davis et al. (2003); Zarraga

Acknowledgements

This work has been funded by the French Defence Procurement Agency (DGA) and has been granted access to the HPC and visualization resources of "Centre de Calcul Interactif" hosted by the University of Nice-Sophia Antipolis.