# Scaling behavior of immersed granular flows

Lhassan Amarsid, Jean-Yves Delenne, Patrick Mutabaruka, Yann Monerie,
Frédéric Perales, Farhang Radjai

► To cite this version:

Lhassan Amarsid, Jean-Yves Delenne, Patrick Mutabaruka, Yann Monerie, Frédéric Perales, et al. Scaling behavior of immersed granular flows. International Conference on Micromechanics of Granular Media (Powders & Grains), Jul 2017, Montpellier, France. 10.1051/epjconf/201714009044. hal-01594573

HAL Id: hal-01594573

https://hal.science/hal-01594573

Submitted on 26 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Scaling behavior of immersed granular flows

L. Amarsid<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,*, J.-Y. Delenne<fnref n="4" />, P. Mutabaruka<fnref n="5" />, Y. Monerie<fnref n="2" />,<fnref n="3" />, F. Perales<fnref n="1" />,<fnref n="3" />, and F. Radja<fnref n="2" />,<fnref n="5" />,<fnref n="3" />

**Abstract.** 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.

## 1 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–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 $\eta_f$, fluid density $\rho_f$, particle mean size $d$, particle density $\rho_s$, confining pressure (acting on the particle phase) $\sigma_s$ and shear rate $\dot{\gamma}$. In the absence of the fluid, it is well-known that the parameter space can be reduced to the inertial number $I = \dot{\gamma}d(\rho_s/\sigma_s)^{1/2}$ [10], and the rheology is described by effective friction coefficient $\mu$ and packing fraction $\Phi$ 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 = \eta_f \dot{\gamma}/\sigma_s$, and the rheology is described by $\mu$ and $\Phi$ 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 un-

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 = \rho_s/\rho_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.

## 2 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(\vec{r}, \vec{v}, t)$ of particle positions $\vec{r}$ and velocities $\vec{v}$. The spatio-temporal evolution of $f$ is governed by the Boltzmann equation:

$$ \left( \frac{\partial}{\partial t} + \vec{v} \cdot \frac{\partial}{\partial \vec{r}} + \frac{\vec{F}(\vec{r})}{m} \cdot \frac{\partial}{\partial \vec{v}} \right) f(\vec{r}, \vec{v}, t) = \Omega_{coll} \quad (1) $$

where $m$ is particle mass, $\vec{F}(\vec{r})$ represents external forces, and $\Omega_{coll}$ is the collision operator. The simplest collision

<footnote n="1">¹IRSN, PSN, CE Cadarache, BP3-13115 St Paul-Lez-Durance Cedex, France</footnote>

<footnote n="2">²LMGC, CNRS - University of Montpellier, 163 rue Auguste Broussonnet, 34090 Montpellier, France</footnote>

<footnote n="3">³Laboratoire MIST, IRSN-CNRS-University of Montpellier, France</footnote>

<footnote n="4">⁴IATE, UMR1208 INRA - CIRAD - University of Montpellier - SupAgro, 34060 Montpellier, France</footnote>

<footnote n="5">⁵⟨MSE⟩², UMI CNRS-MIT, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge 02139, USA</footnote>

<footnote>*e-mail: lhassan.amarsid@umontpellier.fr</footnote>
model is the BGK operator [13]:

$$ \Omega_{BGK} = -\frac{1}{\tau}(f - f^0) \qquad (2) $$

where $\tau$ 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–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 $\vec{v}_i$.

<figure><img src="image_2.png" /><figcaption>Figure 1. D2Q9 scheme.</figcaption></figure>

The discretized equations for different directions are solved in two steps:

$$ \begin{array}{ll} \text{Collision :} & f_i^{\text{out}}(\vec{r}, t) = f_i(\vec{r}, t) + \Omega_i \\ \text{Streaming :} & f_i(\vec{r} + \Delta t \vec{v}_i, t + \Delta t) = f_i^{\text{out}}(\vec{r}, t) \end{array} \qquad (3) $$

where $\Delta t$ is the time step, and the $f_i^{\text{out}}(\vec{r}, t)$ are the distribution functions after the collision step. The density $\rho(\vec{r}, t)$ and fluid velocity $\vec{u}(\vec{r}, t)$ are obtained as follows:

$$ \begin{aligned} \rho(\vec{r}, t) &= \sum_i f_i(\vec{r}, t) \\ \rho(\vec{r}, t) \vec{u}(\vec{r}, t) &= \sum_i f_i(\vec{r}, t) \vec{v}_i \end{aligned} \qquad (4) $$

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

$$ \mathbf{m} = \mathbf{Mf} \qquad (5) $$

where $\mathbf{m} = (m_0, m_1, \dots, m_8)^T$ is the moment vector, $\mathbf{f} = (f_0, f_1, \dots, f_8)^T$ and M is given by:

$$ \mathbf{M} = \begin{pmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ -4 & -1 & 2 & -1 & 2 & -1 & 2 & -1 & 2 \\ 4 & -2 & 1 & -2 & 1 & -2 & 1 & -2 & 1 \\ 0 & 1 & 1 & 0 & -1 & -1 & -1 & 0 & 1 \\ 0 & -2 & 1 & 0 & -1 & 2 & -1 & 0 & 1 \\ 0 & 0 & 1 & 1 & 1 & 0 & -1 & -1 & -1 \\ 0 & 0 & 1 & -2 & 1 & 0 & -1 & 2 & -1 \\ 0 & 1 & 0 & -1 & 0 & 1 & 0 & -1 & 0 \\ 0 & 0 & 1 & 0 & -1 & 0 & 1 & 0 & -1 \end{pmatrix} $$

Hence, the collision step is applied in the moment space, each moment $m_i$ being relaxed to its equilibrium state $m_i^{\text{eq}}$

with a relaxation time $s_i$. The moments corresponding to the density $\rho(\vec{r}, t)$ and the flux $\vec{j}(\vec{r}, t) = \rho(\vec{r}, t)\vec{u}(\vec{r}, t)$ are conserved. The moment vector $\mathbf{m}^{\text{out}}$ resulting from collision can be written as

$$ \mathbf{m}^{\text{out}} = \mathbf{m} - \mathbf{S}(\mathbf{m} - \mathbf{m}^{\text{eq}}) \qquad (6) $$

where $\mathbf{S} = \text{diag}(0, s_1, s_2, 0, s_4, 0, s_6, s_7, s_8)$ is a diagonal 9×9 matrix.

All relaxation times are proportional to $\tau^{-1}$ [20]. The equilibrium moment vector $\mathbf{m}^{\text{eq}}$ is given by

$$ \mathbf{m}^{\text{eq}} = \begin{bmatrix} \rho \\ -2\rho + 3 \left( j_x^2 + j_y^2 \right) / \rho \\ \rho - 3 \left( j_x^2 + j_y^2 \right) / \rho \\ j_x \\ -j_x \\ j_y \\ -j_y \\ \left( j_x^2 - j_y^2 \right) / \rho \\ j_x j_y / \rho \end{bmatrix} \qquad (7) $$

The distribution functions $f_i^{\text{out}}(\vec{r}, t)$ after the collision step are $\mathbf{f}^{\text{out}} = \mathbf{M}^{-1}\mathbf{m}^{\text{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$ ($\vec{v}_i + \vec{v}_{iopp} = \vec{0}$):

$$ f_{iopp}^{\text{in}}(\vec{r}, t + \Delta t \vec{v}_i) = f_i^{\text{out}}(\vec{r}, t) \qquad (8) $$

The equations of motion of the particles are integrated by means of the MD method [21, 22]. The normal force $\vec{F}_N$ between a pair $(i, j)$ of touching particles is governed by a viscoelastic law:

$$ \vec{F}_N = \begin{cases} (-k_n \delta_{ij} - \gamma_n \delta_{ij}) \vec{n}_{ij} & \text{if } \delta_{ij} < 0 \\ \vec{0} & \text{otherwise} \end{cases} \qquad (9) $$

where $\delta_{ij}$ is the gap or the overlap, $\delta_{ij}$ its derivative with respect to time, $k_n$ is the stiffness and $\gamma_n$ is a viscous damping parameter that controls restitution coefficient.

The friction force $\vec{F}_T$ obeys the dry Coulomb friction law:

$$ \vec{F}_T = - \min \{ \gamma_t \|\vec{v}_t\| ; \mu_s \|\vec{F}_N\| \} \vec{t}_{ij} \qquad (10) $$

where $\vec{v}_t$ is the tangential velocity at contact, $\gamma_t$ is the tangential viscosity parameter, and $\mu_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].
<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

## 3 Simulated system

The simulated system is shown in Fig. 2. The fluid fills a rectangular domain which is periodic in the horizontal direction. The particles are disks of mean diameter $d = 2.5 \times 10^{-3}$ m and distributed in the range [$d_{min}, d_{max}$] with $d_{min} = 0.6d_{max}$. A confining pressure $\sigma_s$ is applied on the top wall in contact with the granular column. This wall is permeable to the fluid filling a larger domain of constant volume. A velocity gradient $\dot{\gamma}$ is applied to the fluid nodes at the domain boundary to generate a plane Couette flow. The particles are initially distributed randomly in the box and assembled by downward motion of the top mobile wall under the action of the normal stress $\sigma_s$.

All simulations were performed with 1253 particles and a well-resolved fluid with lattice step $< 0.03d$. The friction coefficient $\mu_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 $\dot{\gamma}$ is varied in the range $[0.28, 5.6] s^{-1}$, confining stress $\sigma_p$ in the range $[20, 120]$ Pa, the relative density $r = \rho_s/\rho_f$ in the range $[0.5, 3]$ and the fluid viscosity $\eta_f$ is varied from that of water $\eta_w$ to 2500$\eta_w$. More than 70 simulations with a total CPU time of about $10^5$ hours were carried out.

## 4 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 $\sigma_s$ 2) the viscous stress $\sigma_v \sim \eta_f \dot{\gamma}$ and 3) the inertial stress $\sigma_i \sim \rho_s (d\dot{\gamma})^2$. The corresponding characteristic times are $t_s = d(\rho_s/\sigma_s)^{1/2}$, $t_v = d(\rho_s/\eta_f \dot{\gamma})^{1/2}$ and $t_i = \dot{\gamma}^{-1}$, respectively. We also have the Stokes time $t_{St} = \eta_f/\sigma_s = t_i(t_s/t_v)^2$. From these three time scales two independent dimensionless numbers can be built: $I = t_s/t_i = (\sigma_i/\sigma_s)^{1/2} = \dot{\gamma} d(\rho_s/\sigma_s)^{1/2}$, the inertial number, and $J = t_s/t_v = (\sigma_v/\sigma_s)^{1/2} = (\eta_f \dot{\gamma}/\sigma_s)^{1/2}$, which is the square root of $I_v$. The ratio $St = (t_v/t_i)^2 = \sigma_i/\sigma_v = I^2/J^2 = \rho_s d^2 \dot{\gamma}/\eta_f$ is the Stokes number.

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. Effective friction coefficient $\mu$ (a) and packing fraction $\Phi$ (b) as a function of the modified inertial number $I_m$ for all data points obtained with different values of system parameters. The solid lines are the fitting forms (12) and (13).</figcaption></figure>

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 $(\sigma_i + \alpha\sigma_v)/\sigma_s = I^2 + \alpha_v J^2 = I^2(1 + \alpha_v/St)$, where $\alpha_v$ is a constant to be determined. Taking the square root of this ratio, we get a modified inertial number

$$I_m = I \left(1 + \frac{\alpha}{St}\right)^{1/2} = \left\{ \frac{\rho_s d^2}{\sigma_s} \dot{\gamma}^2 + \alpha \frac{\eta_f}{\sigma_s} \dot{\gamma} \right\}^{1/2}, \quad (11)$$

In this way, the parameter space can be reduced only if the packing $\Phi$ and effective friction coefficient $\mu$ are unique functions of $I_m$ for a constant value of $\alpha$.

Figure 3 shows $\mu$ and $\Phi$ as a function of $I_m$ for $\alpha = 2.0$. We see that, all data points remarkably collapse on well-defined curves both excellently fit by the following functional forms:

$$\mu(I_m) = \mu_c + \frac{\delta\mu}{1 + b/I_m}, \quad (12)$$

$$\Phi(I_m) = \frac{\Phi_c}{1 + aI_m} \quad (13)$$
with $\mu_c = 0.280 \pm 0.002$ the quasi-static effective friction coefficient, $\Phi_c = 0.8123 \pm 0.0003$ the steady-state packing fraction, $b = 0.246 \pm 0.008$, $\delta\mu = 0.783 \pm 0.014$ and $a = 0.750 \pm 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 $St \approx \alpha \approx 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 definition (11) of the modified inertial number, when the shear rate $\dot{\gamma}$ tends to zero, the leading term is the viscous stress (the term in $\dot{\gamma}$) 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 \approx \sqrt{\alpha I_v}$. But the confining stress $\sigma_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.

## 5 Conclusion

Our extensive simulations using a sub-particle well-resolved 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 $\eta_f$ is replaced by $\eta_f + \rho_s d^2 \dot{\gamma} / \alpha$. 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 $\mu_c$, $\Phi_c$, $b$, $\delta\mu$ and $a$. This assumption needs to be checked by further simulations.

## References

[1] R.M. Iverson, Reviews of Geophysics **35**, 245 (1997)

[2] F. Legros, Engineering Geology **63**, 301 (2002)

[3] K. Hewitt, American Scientist **98**, 410 (2010)

[4] K.D. Nguyen, S. Guillou, J. Chauchat, N. Barbry, Advances in Water Resources **32**, 1187 (2009)

[5] C. Lareo, P. Fryer, M. Barigou, Food and Bioproducts Processing **75**, 73 (1997)

[6] M.J. Rhodes, *Introduction to particle technology* (John Wiley & Sons, Chichester, 1998)

[7] J.J. Stickel, R.L. Powell, Annual Review of Fluid Mechanics **37**, 129 (2005)

[8] F. Boyer, E. Guazzelli, O. Pouliquen, Phys. Rev. Lett. **107**, 188301 (2011)

[9] F. Blanc, F. Peters, E. Lemaire, Journal of Rheology **55**, 835 (2011)

[10] GDR-Midi, The European Physical Journal E **14**, 341 (2004)

[11] M. Trulsson, B. Andreotti, P. Claudin, Phys. Rev. Lett. **109**, 118305 (2012)

[12] S. Courrech du Pont, P. Gondret, B. Perrin, M. Rabaud, Physical Review Letters **90**, 044301 (2003)

[13] P.L. Bhatnagar, E.P. Gross, M. Krook, Phys. Rev. **94**, 511 (1954)

[14] C. Bardos, F. Golse, D. Levermore, Journal of Statistical Physics **63**, 323 (1991)

[15] G.R. McNamara, G. Zanetti, Phys. Rev. Lett. **61**, 2332 (1988)

[16] F.J. Higuera, J. Jiménez, EPL (Europhysics Letters) **9**, 663 (1989)

[17] Y.H. Qian, D. D'Humières, P. Lallemand, EPL (Europhysics Letters) **17**, 479 (1992)

[18] D. d'Humières, Progress in Aeronautics and Astronautics **159**, 450+ (1992)

[19] P. Lallemand, L.S. Luo, Phys. Rev. E **61**, 6546 (2000)

[20] A. Mussa, P. Asinari, L.S. Luo, J. Comput. Phys. **228**, 983 (2009)

[21] S. Luding, *Collisions and Contacts between two particles*, in *Physics of dry granular media - NATO ASI Series E350*, edited by H.J. Herrmann, J.P. Hovi, S. Luding (Kluwer Academic Publishers, Dordrecht, 1998), p. 285

[22] P.A. Cundall, O.D.L. Strack, Géotechnique **29**, 47 (1979)

[23] M. Bouzidi, M. Firdaouss, P. Lallemand, Physics of Fluids **13**, 3452 (2001)

[24] A.J.C. Ladd, Journal of Fluid Mechanics **271**, 285 (1994)