# Time scales and rheology of visco-cohesive granular flows

Farhang Radjai

► To cite this version:

Farhang Radjai. Time scales and rheology of visco-cohesive granular flows. EPJ Web of Conferences, 2021, 249, pp.03044. 10.1051/epjconf/202124903044. hal-03453949

HAL Id: hal-03453949

https://hal.science/hal-03453949

Submitted on 29 Nov 2021

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>

Distributed under a Creative Commons Attribution 4.0 International License
# Time scales and rheology of visco-cohesive granular flows

Farhang Radjai¹,*,

¹CNRS, LMGC, University of Montpellier, 163 rue Auguste Broussonnet, F-34095 Montpellier, France

**Abstract.** In the presence of viscous and cohesive interactions between particles, a granular flow is governed by several characteristic time and stress scales that determine its rheological properties (shear stress, packing fraction, effective viscosities). In this paper, we revisit and extend the scaling arguments used previously for dry cohesionless granular flows and suspensions. We show that the rheology can be in principle described by a single dimensionless control parameter that includes all characteristic times. We also briefly present simulation results for 2D sheared suspensions and 3D wet granular flows where the effective friction coefficient and packing fraction are consistently described as functions of this unique control parameter.

## 1 Introduction

Scaling arguments have been successfully used for a long time in the field of granular materials. For example, as perfect rigid particles and friction law involve no intrinsic stress scale, the shear stresses in quasi-static flows are contained in a Coulomb cone, which bears on the description of the shear strength in terms of stress ratios [1]. Hence, the ratio of shear stress to the normal stress is an effective friction coefficient $\mu$ whose evolution with shear strain describes the quasi-static rheology together with that of dilatancy, which due to rate independence in the quasi-static limit, is also a ratio of volumetric strain to shear strain. We refer to such states as quasi-static flow states (QSFS). The scaling argument was also used to rationalize the behavior of inertial flows by considering that a dry cohesionless flow driven by a shear rate $\dot{\gamma}$ involves two time scales, and only their ratio is relevant for the evolution of $\mu$ and $\Phi$ [2]: 1) shear time $t_i = \dot{\gamma}^{-1}$ and 2) relaxation time $t_p = d(\rho_s/\sigma_p)^{1/2}$ of a particle of mass density $\rho_s$ and diameter $d$ under the action of a confining stress $\sigma_p$. The ratio is the inertial number $I = t_p/t_i = \dot{\gamma}d(\rho_s/\sigma_p)^{1/2}$, and the rheology is described by the functions $\mu(I)$ and $\Phi(I)$. In a similar vein, in a suspension of particles immersed in a viscous fluid, if the viscous drag stress $\eta_f\dot{\gamma}$, where $\eta_f$ is the fluid dynamic viscosity, is considered together with a confining stress $\sigma_p$ exerted on the granular phase, the rheology is described by the functions $\mu(J)$ and $\Phi(J)$ of the dimensionless number $J = \eta_f\dot{\gamma}/\sigma_p$ [3].

In the above examples the use of dimensionless numbers stems from a general dimensional analysis (unit independence) or objectivity principle (frame independence requiring $\mu$ to be defined from stress invariants) without referring to the underlying physical mechanisms, which are complex and often characterized by broad statistical distributions of the particle momenta, contact forces, collision times...). However, when more characteristic times

or stresses come into play, a closer examination of such mechanisms becomes necessary. For example, for dense suspensions of particles in a liquid, both particle inertia and liquid viscosity need to be considered in addition to the confining stress. Numerical simulations suggest that the rheology can still be described by a single dimensionless control parameter involving the ratio of a linear combination of viscous and kinetic stresses to the confining stress [4, 5]. On the other hand, in the case of visco-cohesive granular materials, the cohesive stress is added to the confining stress and a single control parameter is defined from four different stress scales [6]. Besides understandable arguments such as delineating the shear-rate dependent internal stresses from shear-rate independent stresses, the rationale behind such an “additive rheology” needs to be clarified from a particle-scale standpoint [7].

In the following, I first discuss a physical picture of steady granular flows as a *continuous perturbation/relaxation process* with a reference state (or ground state) defined by the confining stress. In this framework, I show why mean contact forces can be used to express the particle relaxation times. Then, I introduce the concept of cohesive stress arising from active adhesion forces at the contact points between particles, and argue that the cohesive stress in dense granular flows plays the same role as the confining stress, defining together a time scale that characterizes the reference flow state. Next, I discuss particle transition times for kinetic and viscous stresses and the meaning of a combined visco-inertial transition time. Finally, I show how a single dimensionless number $I_m$ is naturally defined as a number characterizing the ratio of the transition rate to the relaxation rate, and briefly present simulation results in which $\mu$ and $\Phi$ collapse on a master curve when plotted as a function of $I_m$.

<footnote>*e-mail: franck.radjai@umontpellier.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. (a) Simple shear geometry; (b) Local particle environment.</figcaption></figure>

## 2 Perturbation/relaxation process

Let us focus on the simple shear geometry of Fig. 1(a). The granular flow is subjected to periodic boundary conditions along the flow (x direction) and driven by a shear strain rate $\dot{\gamma}$. The vertical stress $\sigma_p$ is the confining pressure, the only stress applied externally on the particles (not on the liquid) and perpendicular to the flow (y direction). The memory of the initial packing state is lost upon shearing, and the steady flow state implies the balance of external and internal stresses. Hence, $\sigma_{yy} = \sigma_p$, and $\sigma_{yx} = \sigma_{xy} = \tau$ is the shear stress. The ratio $\mu = \tau/\sigma_p$ characterizes the steady stress state. Both $\mu$ and $\Phi$ depend on the material properties (friction coefficient between particles, particle shapes and size distributions) and the control parameters $\sigma_p$ and $\dot{\gamma}$. The particles are assumed to be perfectly rigid so that the only material parameters at the particle scale are the mean particle diameter $d$ and the particle density $\rho_s$ or mean particle mass $\rho_s d^3$. Here, we are only interested in the role of control parameters in combination with fixed material parameters.

A strict QSFS means that the particles follow exactly the displacements imposed by the uniform shear rate (affine displacement). In other words, for a particle $i$ at a position $\vec{r}(i)$, the velocity components are $v_x(i) = \dot{\gamma}r_y(i)$ and $v_y(i) = 0$. Such a motion is obviously forbidden due to steric exclusions between particles. Hence, nonaffine displacements $\vec{s}(i)$ are developed to accommodate particle motions [8]. The motion can be described as the sum of affine and nonaffine components: $v_x(i) = \dot{\gamma}r_y(i) + \dot{s}_x(i)$ and $v_y(i) = \dot{s}_y(i)$, with $\langle \dot{s}_x(i) \rangle = \langle \dot{s}_y(i) \rangle = 0$. The nonaffine displacement field is observed to be highly correlated in space, but disappears after a short time before a new pattern appears. This process clearly shows that the overall shear deformation is composed of periods of the buildup of correlated motions and short periods of relaxation. At the scale of a single particle, the nonaffine displacements induced by the motions of all particles implies the variation of the contact forces acting on the particles by its neighbors. During its motion, the particle probes the force network, which is highly inhomogeneous. The force acting on the particle are balanced during the buildup of collective motions, but becomes unstable whenever the changes in nonaffine components of the displacements are no more

accommodated. The transition between two stable configurations is governed by the shear rate whereas relaxation to a new configuration is controlled by the confining stress.

For the estimation of the relaxation time, we need an evaluation of the force fluctuations. The total force $\vec{F}(i)$ on particle $i$ changes as a result of the evolution of contact forces $\vec{f}_{ij}$ exerted by contact neighbors $j$. As we seek the connection with the imposed stress $\sigma_p$, let us consider the $yy$-component of the tensor moment $M_{yy}(i) = \sum_j f_{y,ij} r_{y,ij}$ of particle $i$, where $r_{y,ij}$ is the $y$-component of the position vector of the contact $(ij)$. The tensor moment inside a volume $V$ is simply the sum of the tensor moments of the particles belonging to $V$ and it can be shown that its density tends to the Cauchy tensor as the number of particles increases $1/V \sum_i M_{yy}(i) \rightarrow \sigma_{yy}$ [9]. Near equilibrium, where the forces are nearly balanced, the tensor moment $M(i)$ of each particle is independent of the choice of the origin for the contact position vectors $\vec{r}_{ij}$. If the origins of the contact vectors for each particle are taken at their centers, the particle tensor moments are $M_{yy}(i) = Z(i)\langle f_y c_y \rangle_i$ where $\vec{c}_{ij} = \vec{r}_i - \vec{r}_{ij}$ is the contact vector (drawn from particle center to the contact point), the average $\langle \cdots \rangle_i$ denotes summation over all contacts of particle $i$, and $Z(i)$ is the number of its contacts. The Cauchy stress tensor (here $yy$-element) can therefore be expressed as

$$ \sigma_{yy} = n_p \langle M_{yy}(i) \rangle = n_p \langle Z(i) \langle f_y c_y \rangle_i \rangle = 2n_c \langle f_y c_y \rangle \quad (1) $$

where $n_p = N_p/V$ is the number density of particles and $n_c = N_c/V = Zn_p/2$ is the number density of contacts. This expression means that the stress tensor is an ensemble average (NPT ensemble) over single particle tensor moments weighted by their numbers of contacts $Z(i)$ (local coordination numbers). It can be described statistically by considering one particle and the probability distributions of forces, coordination numbers and contact directions.

Let us now consider the fluctuations $\delta f_y$ of the vertical force component. To simplify the derivation, we assume spherical particles of average radius $R$ although it can easily be extended to more general particle shapes. We thus have $c_y = R \sin \theta$. Note that $f_y = (f_n \cos \theta + f_i \sin \theta)$ is a signed variable and $\sum_j f_{y,ij} \approx 0$ due to equilibrium. Hence $\langle f_y \rangle = 0$ and $\delta f_y = \langle f_y^2 \rangle^{1/2}$. The tensor moment can be evaluated from the components under the general assumption that the forces and contact orientations are not correlated with the contact directions $\theta$, implying an isotropic contact network. It can be shown that the contact network anisotropy is of second order. Therefore $\langle f_y c_y \rangle = R \langle f_n \sin^2 \theta + f_i \sin \theta \cos \theta \rangle = \frac{R}{2} \langle f_n \rangle$. We may alternatively evaluate it from the Pearson correlation $K_{fc} = \langle f_y c_y \rangle / (\delta f_y \delta c_y)$. We have $\delta c_y = \langle R^2 \cos^2 \theta \rangle^{1/2} = R/\sqrt{2}$. Hence,

$$ \delta f_y = \frac{1}{\sqrt{2}K_{fc}} \langle f_n \rangle = \frac{1}{\sqrt{2}K_{fc}} \frac{\sigma_p}{n_c R} = kd^2 \sigma_p \quad (2) $$

with

$$ k = \frac{\sqrt{2\pi}}{3Z\Phi K_{fc}}. \quad (3) $$

This prefactor is of the order 0.4 and its variations are small. Equation (2) shows that the force fluctuations are
proportional to the average normal force and confining stress. This is consistent with the nearly exponential distribution of force components, which implies that standard deviation is equal to the mean. Since $f_y$ is signed, its standard deviation is two times $\langle|f_y|\rangle$. These two viewpoints correspond respectively either to set of particles inside a granular packing or to a moving particle experiencing local force variations.

For the evaluation of the relaxation time one can thus use $d^2\sigma_p$ instead of $\delta f_y$. For an acceleration equal to $a = \delta f_y/m$, the particle moves a distance $at^2/2$ during a time $t$. A coarse-grained relaxation time can be defined assuming fluctuations of the order of $\delta f_y$ and distances of the order of particle diameter $d$. The relaxation time is the time for a particle to move a distance equal to its own diameter [2]:

$$t_x = \left( \frac{md}{\delta f_y} \right)^{1/2} \qquad (4)$$

Hence, for an external confining pressure $\sigma_p$, the force fluctuation is $d^2\sigma_p$ and relaxation time is given by $t_p = \left(\frac{m}{\sigma_p d}\right)^{1/2} = d\left(\frac{\rho_s}{\sigma_p}\right)^{1/2}$, where we set $m \approx \rho_s d^3$.

The particles are driven by the shear rate $\dot{\gamma}$, and the transition time is the time required for a particle to move a distance equal to its own diameter under the action of the kinetic or inertial stress $\sigma_i = \rho_s(d\dot{\gamma})^2$. Replacing $\delta f_y$ by $d^2\sigma_i$ in equation (4), we get the kinetic time $t_i$

$$t_x = \left( \frac{md}{d^2\rho_s(d\dot{\gamma})^2} \right)^{1/2} = \dot{\gamma}^{-1} \qquad (5)$$

which coincides with shear time. The rheology is therefore characterized by the ratio of the transition time to the relaxation time:

$$I = \frac{t_p}{t_i} = \dot{\gamma} \left( \frac{\rho_s}{\sigma_p} \right)^{1/2} \qquad (6)$$

In dense flows $t_p$ is always small compared to $t_i$. The inertial effects begin to affect $\mu$ and $\Phi$ for values of $I$ as low as $10^{-3}$. This means that the stress fluctuations due to structural relaxation under the confining stress are much higher than the kinetic stresses induced by shearing. For this reason, the shear process may be viewed as a continuous perturbation/relaxation process with the unperturbed static packing as the reference state. This process is analogous to stick-slip motion, but with the major difference that the particles never stick. In view of the development of nonaffine displacements, we may also describe the motion as a prediction/correction process controlled by shear rate and confining pressure, respectively.

### 3 Effect of adhesion and viscous forces

Let us add an active normal adhesion force $f_c$ to all contacts in the packing. This is equivalent to applying normal forces $f_c \vec{n}_{ij}$ and $f_c n_{ji} = -f_c \vec{n}_{ij}$ on the particles $i$ and $j$ at their contact point $(ij)$. Replacing $f_y$ by $f_c$ in equation (1), the adhesion forces induce an internal cohesive stress

$$\sigma_c = 2n_c f_c \langle n_y c_y \rangle = n_c d f_c \qquad (7)$$

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Effective friction coefficient in simulated wet materials in quasi-static flow as a function of cohesion number.</figcaption></figure>

Replacing $\delta f_y$ by $\sigma_c d^2$ in equation (4), we get the relaxation time

$$t_c = d \left( \frac{\rho_s}{\sigma_c} \right)^{1/2} \qquad (8)$$

However, the cohesive relaxation time does not make sense alone in the presence of a confining stress $\sigma_p$. Both $\sigma_p$ and $\sigma_c$ tend to counterbalance the effect of shear-induced perturbations. By stress additivity, we therefore consider $\sigma_n = \sigma_p + \sigma_c$ as the reference stress, and the normal relaxation time is given by

$$t_n = d \left( \frac{\rho_s}{\sigma_c + \sigma_p} \right)^{1/2} = \frac{t_p}{(1 + \xi)^{1/2}} \qquad (9)$$

The relative cohesion is defined by the cohesion index $\xi = \sigma_c/\sigma_p$. For the same reason, the effective friction coefficient should now be defined with respect to this reference normal stress:

$$\mu = \frac{\tau}{\sigma_p(1 + \xi)} \qquad (10)$$

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Effective friction $\mu$ and packing fraction $\Phi$ as a function of the control parameter $I_m$ (see text) in a 2D sheared dense suspension simulated by coupled DEM-LBM for a broad range of parameters. In each series of simulations the values of all parameters are kept constant except those with the color and symbol indicated [4].</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Effective friction $\mu$ and packing fraction $\Phi$ as a function of $I_m$ in a sheared packing of wet particles simulated for a broad range of parameters. For each set of simulations, the symbols and their colors correspond to the parameters that are varied with their ranges, all other parameters being kept constant [7].</figcaption></figure>

It should be noted that the cohesive stress is often defined by its order of magnitude. For example, for capillary forces, we may set $\sigma'_c = \frac{\gamma_s}{d} = f_c/(\pi d^2)$, where $\gamma_s$ is the surface tension of the liquid. This differs from the exact expression of $\sigma_c$ given by equation (7). For this reason, in equation (10) $\xi$ should be replaced by $\alpha\xi$. Comparing the two expressions, we get $\alpha = 1/(3Z\Phi)$. In QSFS, we get $\alpha \approx 0.09$ by setting $Z \approx 6$ and $\Phi \approx 0.6$. Hence,

$$ \mu_c = \frac{\tau}{\sigma_p} = \mu(1 + \alpha\xi) \quad (11) $$

This is what we observed from 3D simulations of packings of spherical particles with capillary forces, as shown in Fig. 2 with $\alpha \approx 0.09$ [7]. The viscous forces may be either due to the action of a saturating liquid or lubrication forces acting at the gaps between particles in wet granular materials. The viscous stress is $\sigma_v = \eta_f \dot{\gamma}$. Using equation (4) and replacing $\delta f_y$ by $\eta_f \dot{\gamma}$, we get the corresponding perturbation time

$$ t_v = d \left( \frac{\rho_s}{\eta_f \dot{\gamma}} \right)^{1/2} . \quad (12) $$

However, the perturbation time can not be arbitrarily separated between the kinetic and viscous effects when both effects are simultaneously present. We use again stress additivity to sum the two stresses. As $\eta_f \dot{\gamma}$ and $\rho_s d^2 \dot{\gamma}^2$ represent only the orders of magnitude of the viscous and kinetic stresses, we use a linear weighted combination $\sigma_s = \sigma_i + \beta\sigma_v = \sigma_i(1 + \beta/St)$, where $St = \sigma_i/\sigma_v$ is the Stokes number. The perturbation time is given by

$$ t_s = d \left( \frac{\rho_s}{\sigma_i(1 + \beta/St)} \right)^{1/2} = \frac{t_i}{(1 + \beta/St)^{1/2}}. \quad (13) $$

It can be shown that the value of $\beta$ is 2 in the presence of a saturating liquid and $8/(27Z) \approx 0.08$ for wet materials.

## 4 Additive rheology

We can define a single dimensionless number $I_m$ that characterizes the perturbation/relaxation process in the presence of the four stresses:

$$ I_m = \frac{t_n}{t_s} = I \left( \frac{1 + \beta/St}{1 + \alpha\xi} \right)^{1/2} \quad (14) $$

This control parameter should be considered as a generalization of the inertial number to cohesive and viscous flows. It contains the inertial, viscous and cohesion numbers. Hence, the rheology is expected to be fully described by unique functions $\mu(I_m)$ and $\Phi(I_m)$.

We performed three series of simulations for 1) 2D cohesive-inertial granular flows [6], 2) 2D sheared suspensions [4], and 3) 3D sheared viscous and cohesive flows [7]. In all cases, for broad range of values of liquid viscosity, cohesion number and inertial number all the friction and pacing fraction data points collapse on a master curve when plotted as a function of $I_m$. Figures 4 and 3 show the collapsed data for 3D visco-cohesive flows and 2D suspensions. Interestingly, the data are nicely fit by the same functional forms as in dry granular flows, thus unifying cohesive and cohesionless granular flows, on the one hand, and viscous sheared suspensions and granular flows, on the other hand.

I warmly thank my dear colleagues J.-Y. Delenne and S. Nezamabadi, as well as my ex-students P. Mutabaruka, L. Amarsid, T.T. Vo and N. Berger for insightful discussions on the topics presented in this paper. The simulation data were adapted from cited papers.

## References

[1] F. Radjai, H. Troadec, S. Roux, Key features of granular plasticity, in *Granular Materials: Fundamentals and Applications*, edited by S. Antony, W. Hoyle, Y. Ding (RS.C, Cambridge, 2004), pp. 157–184

[2] GDR-MiDi, Eur. Phys. J. E **14**, 341 (2004)

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

[4] L. Amarsid, J.Y. Delenne, P. Mutabaruka, Y. Monerie, F. Perales, F. Radjai, Phys. Rev. E **96**, 012901 (2017)

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

[6] N. Berger, E. Azéma, J.F. Douce, F. Radjai, EPL-Europhysics Letters **112** (2015)

[7] T. Vo, S. Nezamabadi, P. Mutabaruka, J.Y. Delenne, F. Radjai, Nature Communications **11**, 1 (2020)

[8] F. Radjai, S. Roux, Phys Rev Lett **89**, 064302 (2002)

[9] L. Staron, F. Radjai, J. Vilotte, Eur. Phys. J. E **18**, 311 (2005)