Donia Marzougui
email: donia.marzougui@3sr-grenoble.fr

Bruno Chareyre
email: bruno.chareyre@3sr-grenoble.fr

Julien Chauchat
email: julien.chauchat@grenoble-inp.fr

Microscopic origins of shear stress in dense fluid-grain mixtures

Keywords: granular suspension, rheology, lubrication, shear flow, discrete element method, hydromechanical coupling

A numerical model is used to simulate rheometer experiments at constant normal stress on dense suspensions of spheres. The complete model includes spheresphere contacts using a soft contact approach, short range hydrodynamic interactions defined by frame-invariant expressions of forces and torques in the lubrication approximation, and drag forces resulting from the poromechanical coupling computed with the DEM-PFV technique. Series of simulations in which some of the coupling terms are neglected highlight the role of the poromechanical coupling in the transient regimes. They also reveal that the shear component of the lubrication forces, though frequently neglected in the literature, has a dominant effect in the volume changes. On the other hand, the effects of lubrication torques are much less significant.

The bulk shear stress is decomposed into contact stress and hydrodynamic stress terms whose dependency on a dimensionless shear rate -the so called viscous number I v -are examined. Both contributions are increasing functions of I v , contacts contribution dominates at low viscous number (I v <0.15) whereas lubrication contributions are dominant for I v > 0.15, consistently with a phenomenological law infered by other authors. Statistics of microstructural variables highlight a complex interplay between solid contacts and hydrodynamic interactions. In contrast with a popular idea, the results suggest that lubrication may not necessarily reduce the contribution of contact forces to the bulk shear stress. The proposed model is general and applies directly to sheared immersed granular media in which pore pressure feedback plays a key role (triggering of avalanches, liquefaction).

Introduction

Dense suspensions of particles immersed in a viscous fluid are ubiquitous in natural phenomenon, such as sediment transport or debris flows, and in numerous industrial applications such as form filling with fluid, concrete in civil engineering or slurry transport in petroleum industries. The understanding of dense suspension rheology has lead to an important research effort over the past decades [START_REF] Bagnold | Experiments on a gravity-free dispersion of large solid spheres in a newtonian fluid under shear[END_REF][START_REF] Frankel | On the viscosity of a concentrated suspension of solid spheres[END_REF][START_REF] Brady | Stokesian dynamics[END_REF]. The complexity of this problem arises from its two-phase nature involving a fluid phase (continuous) and a particulate phase (discrete) for which particle-particle interactions and fluid-particle interactions contribute to the behavior of the system in the dense limit.

Classical rheometer experiments impose simple shear of the suspensions at constant volume (type I). In such case, the interpretation of the shear stress in terms of effective viscosity η e suggests that η e diverges at high solid fraction (φ ≈ 0.6 for spheres) [START_REF] Stickel | Fluid mechanics and rheology of dense suspensions[END_REF]. Rheometer experiments at constant normal stress P have been performed only recently [START_REF] Boyer | Unifying suspension and granular rheology[END_REF] (type II). Under such conditions the volume of the suspension is free to change as a response to the imposed shearing. On this basis, a description of dense suspensions has been proposed, which unifies classical suspension rheology, described in terms of a shear and normal effective viscosity which depend on the solid fraction, and the dense granular flow rheology, described in terms of shear to normal stress ratio (µ) and solid fraction (φ). For this purpose a dimensionless shear rate was introduced, the so called viscous number (I v ), controlling both frictional and viscous contributions to the shear stress. An advantage of such a visco-plastic vision is that the behaviour of the suspension in the very dense limit is not associated with a divergence of the viscosity. Instead, the shear to normal stress ratio reaches a constant value corresponding to a Coulomb-type bulk friction.

From an analytical point of view, the rheology of suspensions has been studied since the beginning of the 20st century. Einstein in 1905 [START_REF] Einstein | Über die von der molekularkinetischen theorie der wärme geforderte bewegung von in ruhenden flüssigkeiten suspendierten teilchen[END_REF] derived the effective viscosity of a dilute suspension based on long range hydrodynamic interactions. Frankel and Acrivos [START_REF] Frankel | On the viscosity of a concentrated suspension of solid spheres[END_REF] proposed another derivation in the limit of dense packings. The later was based on the so-called lubrication approximation of the hydrodynamic interactions between nearly touching spheres. The lubrication terms are singular, they appear as pair interactions between particles and they diverge at the approach of contact.

The first discrete numerical simulations of particle suspension has been proposed in the framework of 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] using resistance and mobility matrices [START_REF] Jeffrey | The forces and couples acting on two nearly touching spheres in low-reynolds-number flow[END_REF]. This technique is able to quantitatively reproduce the divergence of effective viscosity for solid fraction approaching close packing. In the general framework of SD, the hydrodynamic forces include both long range and short range interactions which are defined independently, the later are the diverging terms as found in the lubrication approximations. In addition the forces depend on fluid velocities and particle veloci-ties as independent variables. Practical implementation of this general framework in computer codes for many particles is still a great challenge, however. Commonly, the fluid velocity unknowns are eliminated by assuming that the fluid comoves with the particles at large scales [START_REF] Bossis | Dynamic simulation of sheared suspensions. i. general method[END_REF][START_REF] Brady | The rheology of concentrated suspensions of spheres in simple shear flow by numerical simulation[END_REF][START_REF] Durlofsky | Dynamic simulations of hydrodynamically interacting particles[END_REF][START_REF] Ladd | Hydrodynamic interactions and the viscosity of suspensions of freely moving spheres[END_REF][START_REF] Ladd | Dynamic simulations of sedimenting spheres[END_REF][START_REF] Cichocki | Stokes drag on conglomerates of spheres[END_REF][START_REF] Ball | A simulation technique for many spheres in quasi-static motion under frameinvariant pair drag and brownian forces[END_REF]. This assumption also entails that the pair drag forces must be frame invariant [START_REF] Ball | A simulation technique for many spheres in quasi-static motion under frameinvariant pair drag and brownian forces[END_REF], which exclude many components of the full resistance matrices. In fact, all the long range interactions must be removed, and only the lubrication terms are left.

The assumption of co-movement is acceptable for shear flow at constant volume (type I), but it is otherwise a severe limitation of the numerical technique. One consequence of this assumption is that the solid fraction must be constant over space and time. It is not uncommon in practical applications to violate this condition. This is the case namely in sedimentation, or when the flow of a suspension has a free surface (e.g. debris flow or sheet flow), or when the shear occurs at imposed confining stress (type II experiments). In such case, the divergence of the large scale velocity fields of both phases balance each other, giving rise to long range hydromechanical coupling, also known as poromechanical couplings in porous media theory and soil mechanics [START_REF] Coussy | Poromechanics[END_REF]. This strong coupling governs a range of phenomena such as liquefaction of loose materials or, conversely, delays in the solid-fluid transition in dense materials [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]. Note also that, even if both phases comove at large scales, the implicit fluid velocity field corresponding to the lubrication terms is not divergence free, and therefore not fully consistent. Beside Stokesian Dynamics, numerical methods solving the Stokes or Navier-Stokes equations directly for the interstitial fluid have been introduced [START_REF] Ladd | Lattice-boltzmann simulations of particle-fluid suspensions[END_REF][START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF]. In all cases, the so-called lubrication terms are singular and can not be solved explicitly by computational fluid dynamics (CFD). Capturing the divergence of these terms when particles approach contact would need to shrink the mesh to unrealistically small element sizes around the contact region. A practical approach is to let CFD compute the non-singular terms and to add the lubrication terms directly using analytical expressions [START_REF] Nguyen | Lubrication corrections for lattice-boltzmann simulations of particle suspensions[END_REF][START_REF] Ding | Extension of the latticeboltzmann method for direct simulation of suspended particles near contact[END_REF].

In the present work, dense suspensions are simulated, using a particle-based method and accounting for three effects: contact interactions, drag forces resulting from the poromechanical coupling, and lubrication forces. The contact forces and the motion of the particles are computed using a Discrete Element Method (DEM). The poromechanical coupling is accounted for using the pore-scale coupling DEM-PFV (Pore Finite Volume) developed recently [START_REF] Chareyre | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF][START_REF] Catalano | Porescale modeling of fluid-particles interaction and emerging poromechanical effects[END_REF]. The DEM-PVF code has been extended to periodic boundary problems for the purpose of the present study. And finally the lubrication forces are introduced using frame invariant expressions.

One question that we will examine in this paper is whether accounting for the divergence free nature of fluid velocity field at the local scale can significantly affect the rheology of dense suspensions at constant volume. The poromechanical coupling will be exhibited as a transient effect during volume changes. Another question of interest concerns the definition of the lubrication forces and torques. As particles move one relative to one another, normal, tangential, rolling and twisting motions generate different effects. These effects are sometimes introduced selectively in numerical models, considering that some of them must have negligible effects, a priori. Typically, only normal and shear forces are computed [START_REF] Rognon | Flowing resistance and dilatancy of dense suspensions: lubrication and repulsion[END_REF][START_REF] Nguyen | Lubrication corrections for lattice-boltzmann simulations of particle suspensions[END_REF][START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF], and sometimes even the shear force is neglected [START_REF] Lemaître | What do dry granular flows tell us about dense non-brownian suspension rheology?[END_REF][START_REF] Ancey | A theoretical framework for granular suspensions in a steady simple shear flow[END_REF][START_REF] Seto | Discontinuous shear thickening of frictional hard-sphere suspensions[END_REF]. Hereafter, we introduce all possible terms in order to evaluate, a posteriori, which ones can be neglected. Finally, the contributions ofrespectively -the solid contacts, the lubrication forces, and the poromechanical coupling will be investigated.

Numerical model

Discrete Element Model

An explicit finite difference scheme is employed for updating the position of each particle in a time-marching algorithm. The particles move according to the Newton's second law. In the absence of gravitational effects, the motion is driven by elastic-frictional contact forces defined using a soft contact approach classical in the DEM [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. The contact parameters are the normal and shear stiffnesses k n and k s , and the angle of contact friction φ. The contact forces are supplemented hereafter with forces coming from the interstitial fluid. A three-dimensional implementation of the DEM as found in the open source software YADE is used herein. For more details about the implementation, please refer to [START_REF] Šmilauer | Yade Dem Formulation[END_REF].

DEM-PFV coupled model

The DEM-PFV method is used to solve a pore-scale version of the mass balance equation which appears in the continuous theory of porous media and leads to the so-called poromechanical coupling. Only the main steps of the method are outlined hereafter since the details can be found in previous papers. Here, we assume incompressible phases as in [START_REF] Chareyre | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF][START_REF] Catalano | Porescale modeling of fluid-particles interaction and emerging poromechanical effects[END_REF] (for compressible phases see [START_REF] Scholtes | Modeling wave-induced pore pressure and effective stress in a granular seabed[END_REF]). A tetrahedral decomposition of the pore space is introduced based on regular triangulation (figure 1), where that part of a tetrahedron occupied by the fluid is called a pore. From now on V i denotes the volume of pore i. It is uniquely defined by the positions x i and sizes of the solid particles, while the rate of change Vi also depends on their velocities ẋi .

An exchange of fluid between adjacent pores i and j is represented by the interface flux q ij . The volume balance equation for one one pore leads to

Vi = j=4 j=1 q ij . ( 1
)
Assuming a Stokes regime entails a linear relation between q ij and the local pressure gradient (p i -p j )/l ij between two pores, where l ij is a reference length (see [START_REF] Chareyre | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF]). It leads to the following relationship between pressure and rate of volume change:

Vi = j=4 j=1 k ij (p j -p i )/l ij = j=4 j=1 K ij (p j -p i ). ( 2
)
In this equation K ij is the local hydraulic conductivity.

It must reflects the small scale geometry of the packing. In details, the proposed expression of K ij depends on a local hydraulic radius R h ij (area of the fluid-solid interface divided by the fluid volume -again see [START_REF] Chareyre | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF]) as

K ij = α S f ij R h ij 2 2η l ij (3)
where S f ij is the cross-sectional area of the pore-throat, η is the viscosity of the fluid, and α can be interpreted as a calibration parameter. α = 1 is known to give good estimates of the actual permeability of glass beads [START_REF] Tong | Pore-Scale Flow Simulations: Model Predictions Compared with Experiments on Bi-Dispersed Granular Assemblies[END_REF] but we used α < 1 in this study. This is further discussed in section 3.

Substituting Vi by its expression in terms of particles velocity and writing equation 2 for every element gives a system of linear equations. In a matrix form and including boundary conditions, the sytem reads

KP = E ẋ + Q q + Q p , ( 4
)
where K is the conductivity matrix containing the terms K ij , P the column vector of pressure unknowns, and E is the matrix defining the rates of volume change of the elements such that Vi = (E ẋ) i . Q q and Q p are flux vectors reflecting the boundary conditions, respectively source terms (imposed fluxes) in Q q and imposed pressures in Q p . Solved at each time step, equation 4 gives the discrete field of fluid pressure P as function of the particles velocity.

The drag forces can be deduced from the pressure field. They are the integrals of the pressure p and the viscous stress τ on the surface of the particle (for a non-inertial fluid, the second integral can be evaluated based on P using momentum balance):

F f k = ∂Γ k pnds + ∂Γ k τ nds (5)
The forces F f k are linearly dependant of P, hence a matrix form giving the drag forces for all particles

F f = I f P (6)
The F f k are introduced in Newton's second law together with the forces coming from solid contacts (F c ) and lubrication effects (F L defined in the next section). I.e.

Mẍ = F c + F L + I f P, (7)
The strong two-way coupling defined by equations 4 and 7 is the poromechanical coupling. It is integrated with an explicit scheme whose accuracy has been verified in [START_REF] Catalano | Porescale modeling of fluid-particles interaction and emerging poromechanical effects[END_REF].

Lubrication

We assume that lubrication reflects the dominant dissipation process locally in the sheared pore fluid. The lubrication effects are computed between particle pairs as soon as they share an edge in the triangulation (figure 1) (the triangulation is updated dynamically during the simulation). They are defined for all the elementary motions described in figure 2. We consider particles k and k ′ with radii a k and a k ′ , linear velocities v k and v k ′ and angular velocities ω k and ω k ′ , respectively. Their average radius is defined as a = (a k + a k ′ )/2 and h denotes the inter-particle distance (surface to surface). An arbitrary relative motion between two particles can be decomposed in four elementary motions corresponding to normal displacement (subscript n), shear displacement (s), rolling (r) and twisting (t). This decomposition is illustrated in figure 2. In addition, we introduce the angular velocity of the local frame attached to the interacting pair:

ω n = (v k ′ -v k ) × n/(a k + a k ′ + h).
Lubrication forces and torques induced by these elementary motions are defined as follow:

F L n = 3 2 π η a 2 h v n (8)
F L s = πη 2 -2a + (2a + h) ln 2a + h h v t (9)
C L r = π η a 3 3 2 ln a h + 63 500 h a ln a h [(ω k -ω k ′ ) × n] (10)
C L t = π η a 2 h a ln a h [(ω k -ω k ′ ) • n] n (11)
where

v n = ((v k ′ -v k ) • n) n is the normal relative velocity and v t = (a k (ω k -ω n )+a k ′ (ω k ′ -ω n ))
×n is an objective expression of the tangential relative velocity. In this set of equations, the normal and shear forces, F n and F s , are based on Frankel & Acrivos [START_REF] Frankel | On the viscosity of a concentrated suspension of solid spheres[END_REF][START_REF] Van Den Brule | Modeling of concentrated suspensions[END_REF] whereas C r and C t are based on Jeffrey & Onishi [START_REF] Jeffrey | The forces and couples acting on two nearly touching spheres in low-reynolds-number flow[END_REF][START_REF] Jeffrey | Calculation of the resistance and mobility functions for two unequal rigid spheres in low-reynolds-number flow[END_REF] -the reason of this choice will be discussed later. The total lubrication force F L k (resp. F L k ′ ) applied by particle k ′ on particle k (resp. by particle k on particle k ′ ) and the total torque C L k (resp. C L k ′ ) applied by particle k ′ on particle k (resp. by particle k on particle k ′ ) relative to the particle center read:

F L k = -F L k ′ = F n + F s , (12)
C L k = (a k + h 2 ) F s + C r + C t , (13)
C L k ′ = (a k ′ + h 2 ) F s -C r -C t . (14)
Note that F s contributes to the total torques. In order to check the validity of the different lubrication approximations for different interparticle distance h 3D Finite Element Method (FEM) simulations of Stokes flow have been carried out. Periodic boundary conditions were used to represent an infinite array of identical spheres, fixed in space but all rotating at the same velocity (inset of figure 3). Figure 3 presents the comparison of the FEM results with that of equation [START_REF] Chareyre | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF], where F L s is determined alternatively using the expression from Jeffrey & Onishi [START_REF] Jeffrey | The forces and couples acting on two nearly touching spheres in low-reynolds-number flow[END_REF]:

F s = π η a ln a h v t , = + + + h Normal motion

Shear motion

Rolling motion

Twist motion

Fig. 2 Relative motion between particles.

or from Frankel & Acrivos [START_REF] Frankel | On the viscosity of a concentrated suspension of solid spheres[END_REF]:

F s = π η 2 -2 a + (2a + h) ln 2a + h h v t .
Both expressions are asymptotically equivalent for h → 0. However, our results show that the second one is in much better agreement with the FEM result for small but finite distances (h/2a < 0.1). Both expressions underestimate the FEM result for h/2a > 0.1, which is not surprising keeping in mind that they have been obtained from asymptotic expansions. However, an additional defect of the former is that it leads to negative torques (i.e. torques with the same sign as the angular velocity) for large h. This feature is unphysical as it leads to a net creation of energy and can severely alter the stability of the numerical scheme. It was thus concluded that the expressions of Frankel & Acrivos were more suitable for implementation.

We account for the deformability of the particles near the contact region by combining the above normal lubrication model with a linear elastic model via a Maxwell-type visco-elastic scheme ( [START_REF] Marzougui | Numerical simulation of dense suspension rheology using a dem-fluid coupled model[END_REF], partly inspired by [START_REF] Rognon | Internal relaxation time in immersed particulate materials[END_REF]) (figure 4). The parameters are k n the contact stiffness and ν n (h) is the instantaneous viscosity of the interaction as defined in eq. ( 3), such that

F L n = ν n (h)v n .
The real velocity of approach between the two surfaces is ḣ = v n -ue n , where u e n is the elastic deformation given by u e n = F L n /k n . Hence, the evolution of the normal lubrication force obeys the differential equation

F L n = ν n (h) v n - Ḟ L n k n , (15)
that we have to integrate over time-steps using the form Fig. 3 Comparaison of viscous shear torques for the case of rotating sphere in a regular assembly of identical particles. h is the surface-to-surface distance and a is the particle's radius.

Ḟ L n = k n v n - F L n ν n (h) . (16)
The same stiffness k n is used for both the DEM contact model and the visco-elastic lubrication model.

Lastly, when the particles approach each other, the pressure in the gap tends to infinity and the normal lubrication force would theoretically prevent contact. Like many others (see e.g. [START_REF] Rognon | Flowing resistance and dilatancy of dense suspensions: lubrication and repulsion[END_REF]), we assume that contact will actually occur when h is of the order of the particle roughness ε (figure 5). Thus the solid contact model will generate a repulsive interaction even before h = 0. As a result of this combination between a visco-elastic model and roughness at contact, h will never reach 0 practically in numerical simulations. . There is no gravity. The numerical sample is H = 18a high, L = 12a long and l = 12a wide (figure 6). It is first confined between two parallel plates then sheared by moving the top and the bottom plates at constant velocity ±V /2 = 1.5 m/s. The boundary conditions for the top plate are the velocities v x = V /2, v z = 0, the total vertical stress T y = 750 Pa and the fluid pressure p = 0. At the bottom plate, v x = -V /2, v z = 0 and the fluid velocity along the y axis v f y = 0 (impermeable boundary). Periodic boundary conditions are defined along the horizontal axis for both the particles and the fluid. For the latest we impose a null pressure gradient at the macro-scale, i.e. ▽p x = ▽p z = 0 (see [START_REF] Marzougui | Numerical simulation of dense suspension rheology using a dem-fluid coupled model[END_REF]). In order to avoid preferential slip zones near the plates, the first layer of spheres in contact with a plate is fixed to the plates by highly cohesive contacts. We introduce the boundary stress vector T = F /S where F is the total force on the top plate and S is the horizontal cross sectional area. T y is constant during the deformation, while T x is a result of the imposed shear.

Transient vs. steady state

Figure 7 shows the evolution of the shear stress T x , the pressure p and the solid fraction φ function of the deformation γ(T ) = T 0 γ(t)dt where γ(t) = V /H(t) is the shear rate. The numerical results are presented for two cases: a first case where the poromechanical coupling is considered and a second one without it (i.e. ignoring the last term in equation 7). When the shear velocity is applied, a transient regime is observed, characterized by an increase followed by a decrease of the shear stress, a decrease of the solid fraction and a negative pressure of the fluid in the coupled case. This later effect entails a higher effective stress in the coupled problem, explaining why the shear stress reaches higher values. The system evolves toward a steady state for large deformations, in which the shear stress and the solid fraction are approximately constant and the pressure is nearly zero. The poromechanical coupling has no visible effect at steady state: the shear stress and the solid fraction reach similar values for both cases.

It is to be noted that the poromechanical coupling entails long range effects in the system and, ultimately, a dependency on the problem size (H in our case). It is known since Terzaghi that the characteristic time of such process scales with ηH 2 /κ where κ is the intrinsic permeability. Since κ scales with a 2 , the relaxation time of the transient regime is proportional to η(H/a) 2 . A consequence is that the peak pore pressure in figure 7 scales with (H/a) 2 , pressure gradients scale with H/a 2 , and finally drag forces scale with H. Important consequences of this feature are that the drag forces and the lubrication forces are not comensurable, and that poromechanical effects cannot be reflected as rheological properties of the bulk material -it is always necessary to solve a coupled problem.

Since our main focus in this study was the bulk viscosity of suspensions (a rheological property), we did not seek a fully realistic combination of drag forces and lubrication forces. Practically, it let us reduce the computation times tremendeously by setting α = 100 in equation 3, while it should be close to 1 for more realistic simulations. The duration of the transient regime would have been multiplied by 100 approximately, and the timesteps of the time marching algorithm would have been reduced by 100, leading to an increase of the total simulation time by a factor 10 4 . Though qualitatively correct, the trends seen in fig. 7 are thus quantitatively wrong (and remember that in any case they are only relevant for a specific value H).

Stress decomposition

Besides T , a tensor representing the average stress in the suspension can be calculated as

σ = σ C + σ L + pI + σ I ( 17
)
σ C = 1 V ij F C ij ⊗ l ij
is the contact stress applied on particles in contact where l ij denotes the branch vector between the centers of the particles i and j.

Similarly, σ L = 1 V ij F L ij ⊗ l ij
is the contribution from lubrication forces [START_REF] Ancey | A theoretical framework for granular suspensions in a steady simple shear flow[END_REF], which will be further decomposed hereafter by considering separately the normal and shear components of the lubrication force. p is the pressure associated to the poromechanical coupling [START_REF] Catalano | Porescale modeling of fluid-particles interaction and emerging poromechanical effects[END_REF]. σ I = k m k v k ⊗ v k reflect the inertial effects as defined in [START_REF] Savage | The stress tensor in a granular flow at high shear rates[END_REF] where m k is the mass of particle k and v k is its velocity. It can be verified in figure 7 -where both T x and σ xy are plotted -that the two expressions compare consistently. Hereafter, eq. ( 17) will be used to assess the microscale origins of the shear stress and their rate dependency. will be further investigated for different values of the viscous number I v in the following.

Results and Discussion

Viscous number

I v is a dimensionless form of the shear rate [START_REF] Boyer | Unifying suspension and granular rheology[END_REF], reflecting the magnitude of viscous effects, and is defined as:

I v = η γ T y , (18)
The key idea here is that, in the viscous regime and at steady state the stress ratio µ = T x /T y and the solid fraction φ are entirely controlled by this unique parameter. In other words, all possible combinations of confining pressure, fluid viscosity, and shear rate corresponding to a given value of I v should give the same result. In order to confirm this property, three series of simulations were conducted in which the control parameters were changed independently in each series to produce different values of I v . The results are plotted in fig. 9 versus the corresponding values of I v . µ and Φ are nearly the same whatever the method to change I v . This result offers a numerical confirmation of the conclusion of Boyer et al. [START_REF] Boyer | Unifying suspension and granular rheology[END_REF]. This property holds only in the viscous regime, i.e. as long as the inertial effects can be neglected. As suggested in [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF], this condition may be characterized by the ratio I/I v , where I = ρ γ2 /T y is the so called inertial number. Here the ratio is at most I/I v = 5, and it is much smaller in most cases, while [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF] suggest I/I v ≃ 10 for the transition from the viscous to the inertial regime. For the rest of this study it was decided to keep I constant and equal to 0.14, which leaves η as the only free parameter. An exception is the dry case (I v = 0) where I = 0.005.

Dropping terms of the hydrodynamic model

The consequences of neglecting some of the terms defined by eqs. ( 5)-( 11) will be examined, in terms of stress ratio µ and of solid fraction φ. This question is of interest for the developpers of numerical models, since many models found in the literature are not including all terms. As we have seen before, neglecting the poromechanical coupling (eq. 5) has a detrimental effect on the transient state but has no strong effect at steady state. We now focus on the different lubrication terms and how they affect the result at steady state. Figure 10 presents the comparison of the numerical results with the phenomenological laws proposed by Boyer et al [START_REF] Boyer | Unifying suspension and granular rheology[END_REF] The result for I v = 0 corresponds to vanishing hydrodynamic forces. It is obviously independent of the hydrodynamic assumptions, and practically it is the result of a simulation for a dry material. The friction coefficient and the solid fraction obtained in this case closely match the values measured experimentally by Boyer et al. The results obtained with normal lubrication forces only (blue squares) match at least qualitatively the empirical evolution of shear stress with I v . It is to be noted however that the solid fraction obtained with this model is almost constant for I v > 0.02, while the experiments suggest a monotonic decrease. As soon as the shear lubrication forces are included, the results (red circles) get closer to the phenomenological law for both µ(I v ) and φ(I v ). It can be concluded that shear lubrication forces play a key role in the rate dependent dilatancy, and they contribute significantly to the shear stress. The normal lubrication alone lead to a sat-isfactory stress ratio but significantly overestimate the solid fraction. Considering this result, one may expect that simulations at imposed volume (called type I in the introduction) and including only normal lubrication forces would underestimate the shear stress even more than in our case.

Further sophisticating the hydrodynamic model does not yield other significant changes. The rolling torques (yellow stars), the twist torques (green diamond) and the poromechanical coupling (cyan triangles) have only marginal effects. We recall that this conclusion holds at the steady state only. As discussed previously, the poromechanical coupling can significantly interfere in the transient regimes.

Contact stress versus hydrodynamic stress

From now on, the results that will be analyzed are obtained with the full model including all possible hydrodynamic interactions. Figure 11 shows the contribution to the total stress of the different terms in eq. ( 17). The contact forces play a significant role for all the values of I v investigated. The contact stress slightly increases for 0 < I v < 0.1 and saturates to an almost constant value for larger values of I v . Lubrication stresses, both normal and shear components, increase almost linearly, and the shear stress due to the shear components is approximately twice smaller than the stress coming from the normal components. For values of I v ≥ 0.15 the sum of the two lubrication stresses exceeds the contact stress. This result highlights the fact that depending on the value of I v two regimes are observed. At low I v the contact interactions are dominant whereas for I v ≥ 0.15 the lubrication interactions dominate. This result confirms a constitutive property inferred by Boyer et al [START_REF] Boyer | Unifying suspension and granular rheology[END_REF]. Herein, the shear stress in dense suspensions is split into two contributions, one coming from the contacts and represented by the same phenomenological law as in dry granular flow, the other one coming from hydrodynamic interactions similar to a Krieger-Dougherty viscosity. However Boyer et al inferred this stress partition from macroscopic measurements. The present results allow to further assess the respective contributions of contacts and hydrodynamic interactions. Also, the choice of a frictional rheology for the contact stress which saturates for high values of I v is predicted by our discrete numerical simulations. We believe that this property is not trivial. Based on solid fraction at the largest I v indeed (φ ≃ 0.38, far below any values that can be reached in dry quasi-static granular systems), one could expect that no solid contacts persist. Nevertheless, the numerical simulations confirm the proposal of Boyer et al. on this aspect.

Roughness

The particles roughness ǫ appears as a key parameter in the lubrication model since the surface-to-surface distance h may vanish as ǫ → 0, a situation in wich the lubrication terms diverge. This situation is peculiar from a physical point of view but it also causes major troubles from a numerical point of view. Arguably perfectly 11 The stress ratio µ and the decomposition in contact stress, the normal lubrication stress and the shear lubrication stress. The solid line is the phenomenological law of [START_REF] Boyer | Unifying suspension and granular rheology[END_REF].

smooth surfaces are rare and this parameter can be justified on a physical ground. It is not always clear how this should be accounted for in models however, and in our case we don't have a precise knowledge of what value of ǫ would be relevant for the spheres used by Boyer et al. In order to evaluate the role of this parameter simulations were reproduced for three different values of ǫ. They are reported in figure 12, which includes the total stress and the contributions from contacts and lubrication forces.

The lesser the roughness, the larger the contribution of lubrication, as one could expect. A less expected result is that the contribution of contacts is not strongly modified. No clear trend can be distinguished, as the points corresponding to the smallest ǫ may be below or above the others depending on the value I v . This suggests that the different values are simply due to imprecisions in the evaluation of the steady state, and that the contribution of contacts could be independent of roughness. Overall, the difference on the total stress is of the order of 10% or less, showing a relatively moderate effect. An effect on dilatancy is visible, the smoother particles dilate more, but again moderate. All three values seem to reflect the main features of the behaviour in spite of quantitative differences. Roughness does not appear as a key parameter here. The limit ǫ 0 remains as an open question, which can't be studied easily with our method due to numerical difficulties.

Microstructure

In order to link averaged quantities to micro-scale variables, we examine how various quantities depends on the orientation of the particle pairs (figure 13). Based on the orientation of the unit normal, every interaction corresponds to a position on the unit sphere. For an arbitrary point M on the unit sphere it is possible to compute averages of interaction variables. For instance, the average distance between the spheres is defined as: h(M ) = n k ∈dS h(k)/N M where dS is a small angular sector centered on M and N M is the number of interactions associated to dS. In figure 13 this value is normalized by the particle diameter (2a). We define similarly the average normal velocity and average shear velocity normalized by 2a γ. The other plots of figure 13 are density functions. The density of contacts is obtained by counting the number N c of interactions with h < ǫ in dS, and dividing by the total number of spheres N s , so that P (M ) = N c (M )/(dS N s ). The densities of the lubrication stress term σ L (M ) and of the contact stress term σ C (M ) are obtained by restricting the sums defined for eq. ( 17) to the subset of interactions associated to M . For instance σ

L (M ) = 1 V dS nij ∈dS F L ij ⊗ l ij .
On figure 13 are the component xy of both stress tensors, normalized by the confining pressure. The lubrication stress is further decomposed into one part due to normal forces and another part due to shear forces.

Despite the fact that all the functions introduced above define surfaces in the 3D space, in figure 13 Figure 13.a shows that there is a minimum in the density of contacts near θ = π/4 for all cases. Conversely, the higher density is observed for orientations between 3π/4 and π. A noticeable effect of increasing I v is that the number of contacts near θ = π/4 vanish. This peculiar effect makes a clear difference between the viscous regime and the dry regime as P (θ = π/4) is always strictly positive in the latest.

The average distance (figure 13.b) is increasing with I v on overall. Since increasing I v corresponds to a decrease of solid fraction, this trend is not surprising. The average distance is anisotropic and takes larger values near θ = π/4, consistently with the lower density of contacts in this region.

The normal component of normalized relative velocity (figure 13.c) is positive on [0, π/2] and negative on [π/2, π]. The extrema are of the same order in both cases, although it can be noted that the velocity of approaching particles ([π/2, π]) is slightly lower. The shear component (figure 13.d) is positive on -approximately -[π/4, 3π/4] and negative on [0, π/4] ∪ [3π/4, π]. In this case the extrema are clearly different, with the largest relative velocities near θ = π/2. We note that the graph of the shear velocity is not exactly symmetric as the maximum values are in fact a bit before θ π/2. It can be explained by smaller h and more solid contacts preventing sliding when θ > π/2.

The magnitude of the normalized relative velocity is increasing slightly with I v . This is observed for both the normal and the shear components. This trend can be explained by considering the growing average distance between particles. If all particles were simply following the mean field velocity, then the relative velocity would obviously grow with h. Since the local fluctuations of velocity with respect to the mean field are not modifying the relative velocities in average, this correlation holds.

The lubrication forces for a given relative velocity are decreasing functions of h. Thus, it could be anticipated from figure 13.b that the contribution of viscous interactions to the bulk stress is dominated by interactions oriented along θ = 3π/4, which have in average smaller values of h but nearly similar values of normal velocity. The results of figure 13.e show a quite different picture. The density of stress due to the normal lubrication forces is actually slightly larger on [0, π/2] and it matches the shape of the normal velocity closely (with the difference that the sign is always positive due to the branch l in the diadic product F L n ⊗ l). This feature may be explained by the strong correlations between h and v n .

The contribution of the shear lubrication forces to the bulk stress is dominated by the interactions near θ = π/2, consistently with the evolution of average shear velocity (the fact that the contribution vanishes when θ = 0 is an effect of the product F L s ⊗ l with l nearly horizontal). Like v s , the density of stress reaches a pick slightly before θ = π/2.

Considering the cumulated contributions σ L = σ L n + σ L s , the non-symmetry of the density of stress becomes even more visible. The larger contribution is due to interactions in the range [π/4, π/2]. This is also the range in which there are nearly no solid contacts between the particles. Thus the non-symmetry may be due to a complex interplay between the solid contacts and the viscous interactions.

The density of contact stress (figure 13.f) is strongly anisotropic, as expected from the lack of contacts on [0, π/2]. When I v = 0 some contacts in [0, π/2] contribute negatively, but the contribution is small. For I v > 0 this contribution becomes negligible. On overall, the contacts contribute to the bulk stress mainly via repulsive forces in the direction near θ = 3π/4. The contribution of contacts to the bulk stress increases with I v . This feature was already observed in figure 10 on the macroscopic variables.

Conclusion

We described a complete modeling framework for the numerical simulation of dense suspensions. The model includes solid contacts between particles using a soft contact approach, short range hydrodynamic interactions defined by frame-invariant expressions of forces and torques in the lubrication approximation, and the poromechanical coupling solved using the DEM-PFV technique.

Numerical experiments of simple shear at imposed confining normal stress have been conducted in an attempt to reproduce recent rheometer experiments on beads. The simulations are in excellent agreement with the empirical data in terms of bulk shear stress and solid fraction at the steady state, for a range of dimensionless shear rate corresponding to 0 ≤ I v < 0.45.

The poromechanical coupling was shown to have a significant effect in the transient regime when the deformation starts. However, no significant effects of this coupling have been exhibited at steady state.

The results obtained by neglecting some of the lubrication terms leads to the following conclusions. First, the normal lubrication term has the larger contribution to the bulk stress. However, considering this term alone leads to underestimate the shear stress, and with such simplification the model is unable to reflect the change of solid fraction for increasing I v . Combining both normal and shear lubrication terms gives much better results in terms of stress and solid fraction. Further sophistication of the model by including the terms associated with rolling and twisting gives only marginal improvements.

The analysis of the various contributions to the bulk stress: contact forces, hydrodynamic forces and fluid pressure, has lead to the following conclusions. The contribution of contacts to the bulk shear stress in the permanent regime increases with I v . This result may be seen as counter-intuitive. First, as higher I v leads to lower solid fraction, one would expect contacts contribution to decrease progressively with I v and ultimately vanish. Second, the assumption that the lubrication effects would prevent solid contacts in suspensions has been the cornerstone of many theoretical and numerical models in the past. Instead, the simulations reported in the present paper suggest that both the contact stress and the lubrication stress increase monotonically in the range of I v investigated.

As roughness of particules is decreased, the contribution of lubrication forces is increased to some extent. However, it does not lead to a decrease of contact forces, which remain nearly unchanged. Again, it is against the idea that lubrication is preventing solid contacts. The reasoning leading to this idea may be simply flawn due to conceptual mistakes. Namely, lubrication forces are often perceived as repulsive force wereas they are neutral overall, inhibiting the opening and the closure of contacts almost equally as revealed by microstructural variables.

The distribution of micro-structural variables revealed a complex interplay between the contact fabric and the hydrodynamic interactions. The anisotropy of contact orientation appears to be more pronounced in suspensions as compared to dry granular materials, due to the effect of the hydrodynamic interactions. This can explain at least partly why the contribution of contacts to the bulk shear stress is increasing.

This contact stress may reach a maximum and ultimately vanish for larger I v . When does it happen remains an open question. Capturing this transition in numerical simulations is a great challenge for numerical models since the lubrication terms do not reflect appropriately all hydrodynamic interactions in more dilute regimes.

Acknowledgments

This work is supported by the PhD grant awarded to D. Marzougui by the University of Grenoble-Alpes.

Disclosures

This work has been supported by the PhD grant awarded to D. Marzougui by the University of Grenoble-Alpes. The content of the manuscript has not been published in previous publications. The study participants are the first author and co-authors and they consented to publish on June 12th, 2014. No institutional review board was required for this publication. The work as a whole has been approved by the IMEP-2 comitee of doctoral studies at Univ Grenoble Alpes on july 3rd 2011.

Fig. 1

 1 Fig. 1 Regular triangulation in 2D.

Fig. 4

 4 Fig.4Visco-elastic scheme of the interaction between two elastic-like particles.

Fig. 5

 5 Fig.5Evolution of the contact force and of the normalized lubrication force in the normal direction as a function of the gap between two particles. ε defines the roughess of the particles surfaces.

  setup The suspension is represented by a bi-periodic packing made of N = 1000 frictional spheres of average radius a = 0.025 ± 0.01 m. The physical properties are (unless stated otherwise for sensitivity analysis) roughness ε = 0.035 a, density ρ = 2500 kg/m 3 , normal contact stiffness k n /a = 5 × 10 5 Pa, shear stiffness k s = k n /2, and contact friction angle ϕ = 30

Fig. 6

 6 Fig. 6 Simulation cell.

4. 1 Fig. 7 (

 17 Figure 8 includes the evolution of different terms of equation (17) (component xy) as the suspension is sheared at I v = 0.21. The inertial stress σ Ixy (not represented here) is negligible compared with the total stress (σ I xy < 2.5% T x ), which indicates that the suspension is dominated by contacts and viscous interactions in this case (further discussion in the next paragraph). Second, the contact stress contributes to approximately half of the total stress (σ C xy ≈ 50% T x ) whereas the other half is due to the normal and shear lubrication forces (σ LN xy ≈ 30% T x and σ LS xy ≈ 20% T x ). The different contributions

Fig. 8

 8 Fig. 8 Decomposition of the total shear stress in contact stress, normal lubrication stress and shear lubrication stress.

Fig. 9

 9 Fig. 9 Normalized shear stress and solid fraction at steady state versus I v . In each series the change of I v is obtained by changing a different parameter: confining pressure ( ), viscosity (⋆), or shear rate ( ).

  on the basis of experiments. The values reported in this figure have been obtained at steady state by increasing the fluid viscosity while keeping the shear rate γ and the normal stress T y constant. Typical results for T x vs. γ are shown in the inset of the figure for different values of η. Series of simulations have been carried out including different combinations of the hydrodynamic effects. Starting with the simplest case where only solid contacts and normal lubrication forces are present, the other hydrodynamic terms are added one by one: shear lubrication force, rolling torque, twist torque and drag forces due to the poromechanical coupling.

Fig. 10

 10 Fig. 10 Stress ratio µ (top) and the solid fraction φ (bottom) at steady state versus I v . Each symbol represents a different combination of lubrication terms. The solid line is the phenomenological law of Boyer et al. [5]. Inset: the total shear stress for different values of fluid viscosity.

Fig. 12

 12 Fig. 12 Stress ratio µ and solid fraction φ at steady state for different values of the roughness parameter ǫ.

  only the values for M in the (Oxy) plane are plotted. The results are given for three different values of I v : I v = 0 (i.e. a dry medium -green line in the figure), I v = 0.025 (red line) and I v = 0.2 (blue line). All distributions are π-periodic, the comments hereafter refer to the interval [0, π].

Fig. 13

 13 Fig.[START_REF] Ding | Extension of the latticeboltzmann method for direct simulation of suspended particles near contact[END_REF] Distributions of normalized quantities in the (x,y) plane [top/left] for the dry case, I v = 0.025 and I v = 0.2: PDF of contact orientation, normal velocity, shear velocity, lubrication stress which is the sum of the normal and shear lubrication stress, and contact stress.