Scaling behavior of immersed granular flows

The shear behavior of granular materials immersed in a viscous fluid depends on fluid properties (viscosity, density), particle properties (size, density) and boundary conditions (shear rate, confining pressure). Using computational fluid dynamics simulations coupled with molecular dynamics for granular flow, and exploring a broad range of the values of parameters, we show that the parameter space can be reduced to a single parameter that controls the packing fraction and effective friction coefficient. This control parameter is a modified inertial number that incorporates viscous effects.


Introduction
Immersed granular materials occur in many natural processes such as sediment transport, landslides and submarine avalanches, as well as in many applications involving particle-laden fluids in powder technology and food and pharmaceutical industries [1][2][3][4][5][6]. The central issue in the field of dense suspensions of non-Brownian particles is how the suspended particles affect the rheology [7] whereas in the granular community the query is how the inertial, viscous, capillary and lubrication effects in the presence of an interstitial fluid modify the shear behavior and packing fraction in granular flows [8,9].
Fluid-grain flows involve a large parameter space: fluid viscosity η f , fluid density ρ f , particle mean size d, particle density ρ s , confining pressure (acting on the particle phase) σ s and shear rateγ. In the absence of the fluid, it is well-known that the parameter space can be reduced to the inertial number I =γd(ρ s /σ s ) 1/2 [10], and the rheology is described by effective friction coefficient μ and packing fraction Φ as functions of I. On the other hand, it has been suggested that, when particle-inertial effects can be neglected, the parameter space is reduced to the viscous number I v = η fγ /σ s , and the rheology is described by μ and Φ as a function of I v [8].
Recently, it was suggested that a general control parameter combining I and I v may account for both inertial and viscous effects in a unified framework [11]. Trulsson et al. simulated sheared granular materials by means of Molecular Dynamics (MD) and the effect of fluid was accounted for through the application of a drag force on all particle centers. However, since the fluid is introduced through its drag force effect on particles, it is une-mail: lhassan.amarsid@umontpellier.fr clear whether the unified approach accounts for other major fluid parameters such as density and volume effects and lift forces induced by shear. In particular, the relative density r = ρ s /ρ f may also influence the flow regime [12]. Moreover, the fluid in the pore space carries dynamic pore pressures, which may come into play in a dense suspension. For all these reasons, the visco-inertial flow regime needs to be investigated by using computational fluid dynamics coupled with granular dynamics.
In this paper, we use MD simulations for particle dynamics coupled with the Lattice Boltzmann Method (LBM) for the dynamics of the fluid phase to investigate the rheology of dense granular flows immersed in a viscous fluid. The simulations cover a broad range of parameter values. As we shall see, the effective flow properties are well described by a single modified inertial number, in agreement with the approach suggested by Trulsson et al. [11]. We first briefly present the numerical approach and system parameters. Then, we describe our main results and conclude with the scopes of this work.

NUMERICAL METHOD
We employed the MD method interfaced with the Lattice Boltzmann Method (LBM) for the simulations. The fluid is modeled by a time-dependent distribution function f ( r, v, t) of particle positions r and velocities v. The spatio-temporal evolution of f is governed by the Boltzmann equation: where m is particle mass, F( r) represents external forces, and Ω coll is the collision operator. The simplest collision model is the BGK operator [13]: where τ is a relaxation time and f 0 is the Maxwell-Boltzmann distribution. The Boltzmann equation together with the BGK collision operator yields the Navier-Stokes equations [14].
In LBM [15][16][17] the velocity vector space is discretized into a finite number of directions. We used the D2Q9 lattice (two dimensions with nine velocity directions), as shown in Fig. 1. An independent distribution function f i corresponds to each velocity direction e i . The discretized equations for different directions are solved in two steps: where Δt is the time step, and the f out i ( r, t) are the distribution functions after the collision step. The density ρ( r, t) and fluid velocity u( r, t) are obtained as follows: The BGK operator is simple but leads to fluctuating velocity fields. In our simulations, we used a Multi-Relaxation Time (MRT) collision approach [18,19]. Nine moments are attributed to every fluid node, corresponding to the nine distribution functions, through a matrix M such that . . , f 8 ) T and M is given by: Hence, the collision step is applied in the moment space, each moment m i being relaxed to its equilibrium state m eq i with a relaxation time s i . The moments corresponding to the density ρ( r, t) and the flux j( r, t) = ρ( r, t) u( r, t) are conserved. The moment vector m out resulting from collision can be written as All relaxation times are proportional to τ −1 [20]. The equilibrium moment vector m eq is given by The distribution functions f out i ( r, t) after the collision step are f out = M −1 m out . Finally, the streaming step is applied in the velocity space.
The no-slip boundary conditions are implemented using the bounce-back rule, which consists in reflecting back the incoming distribution functions at a boundary node to their original fluid nodes in the opposite direction iopp ( e i + e iopp = 0): The equations of motion of the particles are integrated by means of the MD method [21,22]. The normal force F N between a pair (i, j) of touching particles is governed by a viscoelastic law: where δ i j is the gap or the overlap,δ i j its derivative with respect to time, k n is the stiffness and γ n is a viscous damping parameter that controls restitution coefficient. The friction force F T obeys the dry Coulomb friction law: where v t is the tangential velocity at contact, γ t is the tangential viscosity parameter, and μ s is the friction coefficient.
The particles are meshed on the lattice grid and represented by solid nodes. The interaction between particles and fluid occurs at their interface. The solid nodes are considered as moving boundaries over which the no-slip condition is imposed [23]. The hydrodynamic forces acting on particles are calculated by the momentum exchange method proposed in [24]. γ σ s Figure 2. A snapshot of the suspension and its boundary conditions (left); a snapshot of the force network and negative and positive fluid pressures in blue and red, respectively.

Simulated system
The simulated system is shown in Fig. 2 All simulations were performed with 1253 particles and a well-resolved fluid with lattice step < 0.03d. The friction coefficient μ s is set to 0.4 between the particles and to 0 with the top wall. The bottom wall is made rough by sticking a layer of particles to the bottom of the box. The shear rateγ is varied in the range [0.28, 5.6] s −1 , confining stress σ p in the range [20,120] Pa, the relative density r = ρ s /ρ f in the range [0.5, 3] and the fluid viscosity η f is varied from that of water η w to 2500η w . More than 70 simulations with a total CPU time of about 10 5 hours were carried out.

Results and discussion
In order to reduce the parameter space, we consider the characteristic stresses which, in contrast to characteristic times, have the advantage of being additive. They are of three different origins: 1) the static stress σ s 2) the viscous stress σ v ∼ η fγ and 3) the inertial stress σ i ∼ ρ s (dγ) 2 . The corresponding characteristic times are t s = d(ρ s /σ s ) 1/2 , t v = d(ρ s /η fγ ) 1/2 and t i =γ −1 , respectively. We also have the Stokes time t S t = η f /σ s = t i (t s /t v ) 2 . From these three time scales two independent dimensionless numbers can be build: I = t s /t i = (σ i /σ s ) 1/2 =γd(ρ s /σ s ) 1/2 , the inertial number, and J = t s /t v = (σ v /σ s ) 1/2 = (η fγ /σ s ) 1/2 , which is the square root of  The total shear stress in steady flow is the sum of static, inertial and viscous stresses, and we expect that the rheology is governed by the additive effects of viscous and inertial stresses on particles compared to the static stress. Hence, owing to this additivity, we can build a unique dimensionless number where α v is a constant to be determined. Taking the square root of this ratio, we get a modified inertial number In this way, the parameter space can be reduced only if the packing Φ and effective friction coefficient μ are unique functions of I m for a constant value of α. Figure 3 shows μ and Φ as a function of I m for α = 2.0. We see that, all data points remarkably collapse on well-defined curves both excellently fit by the following functional forms: with μ c = 0.280 ± 0.002 the quasi-static effective friction coefficient, Φ c = 0.8123 ± 0.0003 the steady-state packing fraction, b = 0.246 ± 0.008, δμ = 0.783 ± 0.014 and a = 0.750±0.003. Hence, all system parameters, including the relative fluid-particle density, affect the rheology only through I m . The transition from viscous to inertial regime is governed by the Stokes number. According to the definition of I m , this transition occurs for S t α 2. This result is in agreement and extends the scope of a unique framework introduced by Trulsson et al. [11] to a more general parameter space. Note that, according to the defintion (11) of the modified inertial number, when the shear rateγ tends to zero, the leading term is the viscous stress (the term inγ) irrespective of the values of other parameters. For this reason, at low shear rates, we generally expect a viscous behavior, where the control parameter is I m √ αI v . But the confining stress σ s plays the same role with respect to viscous and inertial forces and its variation does not lead to a transition between inertial and viscous regimes.

Conclusion
Our extensive simulations using a sub-particle wellresolved computation of a viscous fluid coupled to rigid particles suggest that the visco-inertial flow regime of immersed granular materials can be described by a modified inertial number combining the inertial number with the Stokes number. This behavior can also be expressed by using a viscous description of the flow. In this description, both shear viscosity and bulk viscosity can be described by unique function of the packing fraction provided the fluid viscosity η f is replaced by η f + ρ s d 2γ /α. This remarkable scaling reflects the joint effects of particle inertia and fluid viscous forces on the granular microstructure. A detailed analysis of local evolutions of pores pressures and contact networks will be presented elsewhere.
Finally, it is important to note that the unified picture discussed in this work is with regard to the functional dependence of flow properties on those parameters that control particle dynamics. It is thus implicitly assumed that parameters such as particle size distribution and friction coefficient between particles enter the rheology only through their effect on the model parameters μ c , Φ c , b, δμ and a. This assumption needs to be checked by further simulations.