# Additive rheology of complex granular flows

Thanh Trung Vo, Saeid Nezamabadi, Patrick Mutabaruka, Jean-Yves

Delenne, Farhang Radjai

► To cite this version:

Thanh Trung Vo, Saeid Nezamabadi, Patrick Mutabaruka, Jean-Yves Delenne, Farhang Radjai. Additive rheology of complex granular flows. Nature Communications, 2020, 11 (1), 10.1038/s41467-020-15263-3. hal-02540059

HAL Id: hal-02540059

https://hal.inrae.fr/hal-02540059

Submitted on 10 Apr 2020

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
ARTICLE

https://doi.org/10.1038/s41467-020-15263-3

OPEN

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

Check for updates

# Additive rheology of complex granular flows

Thanh Trung Vo<fnref n="1" />,², Saeid Nezamabadi ID ²,<fnref n="3" />, Patrick Mutabaruka ID ², Jean-Yves Delenne<fnref n="3" /> & Farhang Radjai²,⁴✉

Granular flows are omnipresent in nature and industrial processes, but their rheological properties such as apparent friction and packing fraction are still elusive when inertial, cohesive and viscous interactions occur between particles in addition to frictional and elastic forces. Here we report on extensive particle dynamics simulations of such complex flows for a model granular system composed of perfectly rigid particles. We show that, when the apparent friction and packing fraction are normalized by their cohesion-dependent quasistatic values, they are governed by a single dimensionless number that, by virtue of stress additivity, accounts for all interactions. We also find that this dimensionless parameter, as a generalized inertial number, describes the texture variables such as the bond network connectivity and anisotropy. Encompassing various stress sources, this unified framework considerably simplifies and extends the modeling scope for granular dynamics, with potential applications to powder technology and natural flows.

<footnote n="1">¹Bridge and Road Department, Danang Architecture University, Da Nang, Vietnam. ²LMGC, University of Montpellier, CNRS, 34060 Montpellier, France.</footnote>

<footnote n="3">³IATE, UMR1208 INRAE-CIRAD-University of Montpellier-SupAgro, 34060 Montpellier, France. ⁴(MSE)², UMI MIT-CNRS, 77 Massachusetts Avenue, Cambridge, MA 02139, USA. ✉email: franck.radjai@umontpellier.fr</footnote>
The crucial role of granular flows in nature (landslides, debris avalanches, slope failure)¹⁻⁸ and industrial processes (handling powders and granulates, compaction, additive manufacturing)⁹⁻¹³ has been at the focus cross-disciplinary research for more than thirty years¹⁴,¹⁵. Recent progress in theoretical understanding of granular flows has been mainly inspired by collective effects such as force chains and jamming¹⁶⁻²⁵, and by searching for relevant dimensionless control parameters²⁶⁻³⁶.

Basic model granular media are composed of rigid particles with frictional contact interactions and, in contrast to interacting particles at the atomic scale or in colloids, rigid frictional particles are devoid of an intrinsic stress or time scale. The relevant scales are therefore set either externally, such as those arising from a confining pressure $\sigma_p$, or by collective particle motions during flow involving an inertial (or kinetic) pressure $\sigma_i \sim \rho_s \langle d \rangle^2 \dot{y}^2$, where $\rho_s$ is the particle density, $\langle d \rangle$ is the mean particle diameter, and $\dot{y}$ is the shear rate²⁷,³⁷. Hence, in the NPT statistical ensemble (with temperature $T=0$ for a granular material), the apparent friction coefficient $\mu = \sigma_f/\sigma_p$, where $\sigma_f$ is the shear stress, and the packing fraction $\Phi$ are expected to be uniquely dependent on the ratio $I^2 = \sigma_f/\sigma_p$, which accounts for the competing effects of particle inertia and confinement. The dimensionless number $I$ is the inertial number, defined as the ratio of two time scales (relaxation time $\langle d \rangle (\rho_s/\sigma_p)^{1/2}$ under load vs. shear time $\dot{y}^{-1}$), and it was found to unify experimental and numerical data in different flow geometries²⁶.

Most of time, however, the particle interactions are not purely frictional and involve characteristic forces that induce additional internal stresses. A well-known example is the cohesive contact force $f_c$ in fine powders. When a powder flows, the average action of the resulting cohesive stress $\sigma_c \sim f_c/d^2$ is similar to a confining stress, tending to prevent from dilation during flow, to enhance the contact forces and to reduce the relaxation time under load³²,³⁸,³⁹. As the stresses are additive, one may thus take the cohesive stress into account on the same footing as the confining pressure by replacing $\sigma_p$ by a linear combination $\sigma_n = \sigma_p + \alpha\sigma_o$ and therefore replacing $I^2$ by $I_c^2 = \sigma_i/\sigma_n = I^2/(1+\alpha\xi)$, where $\xi = \sigma_f/\sigma_p$ is the cohesion index, and $\alpha$ is a material-dependent parameter³². In the same way, in submerged granular flows, where the viscous drag force $\sigma_v$ is present instead of cohesive stress, $\sigma_i$ may be replaced by a linear combination $\sigma_i + \beta\sigma_v$, leading to a visco-inertial number $I_v^2 = (\sigma_i + \beta\sigma_v)/\sigma_p = I^2(1+\beta/St)$, where $St \equiv \sigma_f/\sigma_v$ is the Stokes number. In dense non-Newtonian suspensions, $I_v$ is found to be the control parameter for both $\mu$ and $\Phi$³¹,³⁵.

The above examples lead to the conjecture that granular flows are fundamentally governed by a single dimensionless parameter combining arbitrary particle interactions by virtue of stress additivity and the role of each interaction with respect to the shear rate and confining pressure. In this paper, we address this interesting issue by simulating wet granular flows such as unsaturated soils and powders at high relative humidity. The liquid bridges between particles induce both capillary and viscous (lubrication) forces whose effects on the flow behavior, together with the confining and inertial stresses, will be quantified for a broad range of parameter values, which are generally difficult to access by means of experiments.

Our results convincingly demonstrate the above conjecture not only for the apparent friction and packing fraction but also for the microstructural variables as a function of a generalized inertial number accounting for the confining, inertial, cohesive and viscous stresses. The results provide also direct evidence for the role of cohesive interactions in dense suspensions when properly interpreted in terms of effective viscosities. This work sets therefore the foundation for a unified description of complex granular flows both encompassing and extending previous work.

## Results

**Particle dynamics simulations and data collapse.** We performed extensive 3D long-shear simulations of granular samples composed of nearly 20,000 spherical particles by means of a particle dynamics method and with a broad range of the values of liquid viscosity $\eta$, surface tension $\gamma_s$, confining stress $\sigma_p$ and shear rate $\dot{y}$. Our results are based on the average values of the stress tensor, velocity fields, packing fraction and granular texture in steady flow of the particles in the simulation cell between the top and bottom walls with periodic boundary conditions in the other directions. Besides repulsive elastic force and friction force, the particle interactions include the approximate analytical expressions of the capillary force $f_c$ and viscous force $f_v$ acting between neighboring particles (see Methods section). The snapshots of Fig. 1 show the boundary conditions, compressive and tensile force chains, contact and non-contact forces (due to capillary bridges) and particle velocities in steady flow.

Figure 2 displays the apparent friction coefficient $\mu$ and packing fraction $\Phi$ as a function of the inertial number $I$ for all our 281 simulations. The confining pressure $\sigma_p$ was varied in the range [15, 1000] Pa, the cohesion index $\xi$ in the range [0, 3.0] (by varying $\gamma_s$ or $\sigma_p$), the liquid viscosity $\eta$ in the range [$\eta_w$, 800$\eta_w$], where $\eta_w$ is the water viscosity, and the shear rate in the wide range [0.31, 10.6]s⁻¹. As in dry flows²⁷, $\mu$ increases and $\Phi$ declines with increasing $I$ but with different values and at

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Simulated system of wet spherical particles. a Particles in the simulation cell composed of a rough immobile bottom wall, a rough top wall subjected to a constant confining stress $\sigma_p$ and moving at constant horizontal velocity $v$, and with periodic boundary conditions along lateral directions. The particle colors are proportional to their diameters.</figcaption></figure>

b Snapshot of compressive (gray) and tensile (blue) force chains in a thin layer parallel to the flow plane. Line thickness is proportional to normal force. c A zoomed-in view of compressive and tensile force chains.

d Snapshot of contact forces (violet) and non-contact capillary forces (green) in a thin window at the center of the flow plane. e A zoomed-in view of the particle velocity field.
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Simulation results for friction and packing fraction. Apparent friction coefficient $\mu$ (a) and packing fraction $\Phi$ (b) as a function of the inertial number $I$ for different sets of system parameters. The data points are average values over the steady state, and the error bars represent their standard deviation in each simulation during steady-state flow. 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.</figcaption></figure>

different rates depending on the viscosity and cohesion index. These differences are observed at both low values (quasi-static flow) and high values (inertial flow) of $I$, and the variability of $\mu$ and $\Phi$ with the variation of $I$ is of the same order of magnitude as with the variation of viscous and cohesive parameters.

The issue is whether all these different values of apparent friction and packing fraction can be expressed as a collapsed function of a single dimensionless number combining surface tension, liquid viscosity, confining pressure and shear rate. In other words, can $I$ be replaced by a more general inertial number $I_m$ that simultaneously accounts for the capillary, viscous and inertial forces? This is indeed what we observe in Fig. 3, displaying all the data points of Fig. 2 as a function of a modified inertial number defined by

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

The values of $\mu$ and $\Phi$ are normalized by $\mu_c$ and $\Phi_c$, respectively, which design their quasi-static values ($I_m \to 0$) and vary linearly with $\xi$, as shown in the two insets to Fig. 3. The parameter values $\alpha \approx 0.062$ and $\beta \approx 0.075$ were determined from two series of simulations but we see that they lead to data collapse for all other simulations, including those of dry cohesionless flows. This means that $\alpha$ and $\beta$ depend only on the material parameters (particle shape and size distribution, friction coefficient between particles) and not on the cohesive and viscous interactions.

**Visco-cohesive inertial number.** The physical argument behind the definition of $I_m$ is the following. There are four characteristic stresses of different origins governing the flow: confining stress $\sigma_p$, inertial stress $\sigma_i$, viscous stress $\sigma_v$, and capillary stress $\sigma_c$. The key

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Rescaled simulation data. Normalized apparent friction coefficient $\mu/\mu_c$ (a) and normalized packing fraction $\phi/\phi_c$ (b) as a function of the generalized inertial number $I_m$ defined by Eq. (1) with $\alpha \approx 0.062$ and $\beta \approx 0.075$. The error bars represent the standard deviation of the data over the steady-state flow. The symbols are the same as those in Fig. 2. The black-solid lines are the analytical expressions of Eqs. (6) and (7). The insets show the evolution of quasistatic values $\mu_c$ and $\Phi_c$ of the apparent friction coefficient and packing fraction, respectively, with the cohesion index $\xi$, and their linear fits (blue-solid lines).</figcaption></figure>

variable for inertial flows is the shear rate. We thus distinguish the characteristic stresses that depend on the shear rate, i.e., $\sigma_i$ and $\sigma_v$, from those that are independent of the shear rate, i.e., $\sigma_p$ and $\sigma_c$. By virtue of stress additivity, the total shear-dependent stress is a linear combination $\sigma_i + \beta\sigma_v$ of the former, and the total shear-independent stress is a linear combination $\sigma_p + \alpha\sigma_c$ of the latter. Hence, the flow variables (apparent friction coefficient and packing fraction) are expected to depend on the ratio $(\sigma_i + \beta\sigma_v)/(\sigma_p + \alpha\sigma_c)$, which simply represents the relative magnitude of the two groups of stresses. We define the generalized inertial number $I_m$ as the square root of this ratio, leading to the expression (1) by setting $I = (\sigma_i/\sigma_p)^{1/2}$, $\xi = \sigma_c/\sigma_p$ and $St = \sigma_i/\sigma_v$.

These primary dimensionless parameters can be evaluated from the system control parameters. The order of magnitude of the viscous stress $\sigma_v$ is conveniently evaluated by replacing the average relative velocity $\theta \sim \dot{y}\langle d \rangle$ induced by shearing in Eq. (15) and considering the dissipated power per unit volume $f_s\delta_{rupt}/\langle d \rangle^3$, where $\delta_{rupt}$ is the debonding distance, yielding $\sigma_v \sim \eta\dot{y}$. The capillary stress is of the order of the capillary force at contact ($\delta=0$) in Eq. (11) divided by the cross section $\langle d \rangle^2$: $\sigma_c \sim \gamma_s/d$. Hence, the dimensionless parameters are given by $I = \dot{y}\langle d \rangle (\rho_s/\sigma_p)^{1/2}$, $\xi = \gamma_s/(\sigma_p\langle d \rangle)$, and $St = \rho_s\langle d \rangle^2\dot{y}/\eta$. Note that all other dimensionless variables can be expressed in terms of these three independent parameters. For example, the capillary number is given by $Ca = \sigma_v/\sigma_c = I^2/(\xi St)$.

**Fitting functions for apparent friction and packing fraction.** According to Eq. (1), $I_m \to 0$ only if $I \to 0$. In this quasistatic limit, the flow variables may still depend on $\xi$, which is the only dimensionless variable in the absence of inertial and viscous stresses. Accordingly, the rheology is expected to be described by
the following general equations:

$$ \mu = \mu_c(\xi) G_{\mu}(I_m), \quad (2) $$

$$ \Phi = \Phi_c(\xi) G_{\Phi}(I_m), \quad (3) $$

These equations are based on functional distinction between the quasistatic limit ($I_m \to 0$) and shear-rate dependent behavior through $I_m$. When $I_m \to 0$, we have $\mu \to \mu_c$ and $\Phi \to \Phi_c$, and thus $G_{\mu} \to 1$ and $G_{\Phi} \to 1$. According to the simulation data displayed in Fig. 3, $\mu_c$ and $\Phi_c$ are linear functions of $\xi$:

$$ \mu_c \approx \mu_0(1 + a\xi) \quad (4) $$

$$ \Phi_c \approx \Phi_0(1 - b\xi) \quad (5) $$

with $a \approx 0.095$ and $b \approx 0.005$. The limit values $\mu_0 \approx 0.392$ and $\Phi_0 \approx 0.594$ are the values of the apparent friction coefficient and packing fraction in the absence of cohesive and viscous forces (dry limit), respectively. Remarkably, the limit value $\Phi_0 \approx 0.594$ obtained here by simulations is equal to the measured value of packing fraction in glass-bead flows$^{29}$.

The data points are nicely fitted by the following functional forms:

$$ \frac{\mu}{\mu_c} = G_{\mu} = 1 + \frac{\Delta_{\mu}}{1 + I_{\mu}/I_m} \quad (6) $$

$$ \frac{\Phi}{\Phi_c} = G_{\Phi} = \frac{1}{1 + I_m / I_{\Phi}} \quad (7) $$

with $\Delta_{\mu} \approx 1.100$, $I_{\mu} \approx 0.095$ and $I_{\Phi} \approx 2.010$. Interestingly, these functions are the same as those previously used for dry granular flows, with the new visco-cohesive inertial parameter $I_m$ replacing $I^{28,40}$. This shows that, up to the values of $I_{\mu}$, $I_{\Phi}$ and $\Delta_{\mu}$, our simulation data are consistent with the experimental measurements of ref. 40 and ref. 28 in the dry case. The values of $I_{\mu}$ and $I_{\Phi}$ are also very close to those obtained in the simulations of Roy et al.34 in a ring shear cell once re-expressed in terms of our definitions. While the functional forms are general$^{28,31,32,34,35}$, the fitting parameters depend on the space dimension and material properties of the granular materials such as particle size distributions, particle shape and friction coefficient between particles. The values of $\alpha$ and $\beta$ reflect the relative roles of viscous, inertial and cohesive forces in collective dissipation mechanisms whereas the values of $I_{\mu}$, $I_{\Phi}$, and $\Delta_{\mu}$ account for the effects of material parameters. Note also that, given the investigated range of values of $I_m$, Eq. (7) can be linearized with an error ~$10^{-3}$, in agreement with refs. 32,34,41.

The fitting forms reveal the double role played by cohesion. Since $I_m$ is a decreasing function of $\xi$, $G_{\mu}$ declines with increasing $\xi$ (dynamic effect) whereas $\mu_c$ increases (quasistatic effect). We can easily check from the parameter values that the quasistatic effect prevails although the dynamic effect becomes important at large values of $\xi$ and $I_m$. These roles are inversed for the packing fraction: $\Phi_c$ declines whereas $G_{\Phi}$ increases when $\xi$ is increased. In other words, the cohesive interactions lead to lower packing fraction (enhanced dilatancy due to cohesive forces) but the inertial effects tend to increase it.

**Transition to the NVT ensemble and effective viscosities.** The flow behavior can alternatively be described in the NVT ensemble (constant-volume shearing) in terms of effective normal and shear viscosities $\eta_n$ and $\eta_t$ defined by $\sigma_n = \eta_n \dot{\gamma}$ and $\sigma_t = \eta_t \dot{\gamma}$, respectively, where $\sigma_t = \mu\sigma_n$ is the shear stress$^{29}$. In this ensemble, the packing fraction $\Phi$ replaces pressure $\sigma_n$ as control parameter, and the rheology is characterized by the functions $\eta_n(\Phi)$ and $\eta_t(\Phi)^{29}$. This is the approach mostly used in experiments on

<figure><img src="image_6.png" /><figcaption>**Fig. 4 Effective viscosities in the NVT ensemble.** Normalized bulk viscosity $c_n$ (a) and shear viscosity $c_t$ (b) as a function of the normalized packing fraction $\Phi/\Phi_c$, where $\Phi_c$ is the quasistatic packing fraction depending linearly on the cohesion index; see Eq. (5). The error bars represent the standard variation over steady-state flow, and the symbols and their colors are the same as those in Fig. 2. The solid lines are the functional forms of Eq. (8) and Eq. (9), respectively.</figcaption></figure>

suspensions. Although our simulations were carried out under NPT conditions, we may deduce $\eta_n(\Phi)$ and $\eta_t(\Phi)$ in the NVT ensemble from $\mu(I_m)$ and $\Phi(I_m)$. Since no external stress is imposed in NVT, the shear stress $\sigma_t$ is a dynamic variable that should scale with the internal shear-dependent stress $\sigma_i + \beta\sigma_v$. Moreover, the NPT and NVT points of view should be compared at the same normal stress state, i.e., $\sigma_n = \sigma_p + \alpha\sigma_c$. Hence, according to Eq. (1), $\sigma_n = c_n(\sigma_i + \beta\sigma_v) \equiv c_n\sigma_n I_m^2 = \eta_n \dot{\gamma}$, implying $c_n = 1/I_m^2 = \eta_n/(\beta\eta + \rho_s\langle d \rangle^2\dot{\gamma})$, and $\sigma_t = c_t(\sigma_i + \beta\sigma_v)$ with $\sigma_t = \mu/I_m^2 = \eta_t/(\beta\eta + \rho_s\langle d \rangle^2\dot{\gamma})$. In this way, in a volume-controlled flow, $c_n$ and $c_t$ represent dimensionless viscosities with $\beta\eta + \rho_s\langle d \rangle^2\dot{\gamma}$ as reference viscosity; see Supplementary Note 3 for more detail.

Figure 4 displays the effective dimensionless viscosities as a function of $\Phi$. We see that all the data points collapse on a master curve when $\Phi$ is normalized by the critical packing fraction $\Phi_c$. Both viscosities diverge as $\Phi \to \Phi_c$ and they are nicely fitted by the analytic expressions

$$ c_n = \frac{1}{I_{\Phi}^2} \left( \frac{\Phi}{\Phi_c - \Phi} \right)^2, \quad (8) $$

$$ c_t = \mu c_n = \frac{1}{I_{\Phi}^2} \left( \frac{\Phi}{\Phi_c - \Phi} \right)^2 \left\{ 1 + \frac{\Delta_{\mu}}{1 + \frac{I_{\mu}}{I_{\Phi} \Phi_c - \Phi}} \right\}, \quad (9) $$

readily deduced from the expressions of $c_n$ and $c_t$ as a function of $I_m$ together with Eqs. (6) and (7). As in suspensions, $1/c_t = I_m^2/\mu$ represents a generalized fluidity parameter of granular flows$^{42}$.

**Scaling of coordination number and bond anisotropy.** Although $I_m$ provides a unified description of the rheology by capturing the effects of particle interactions on the apparent friction coefficient and packing fraction and, alternatively, the effective viscosities, it is essential to check its robustness with
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Scaling of microstructural variables. Number density Z of capillary bonds (a) and bond orientation anisotropy A (b) as a function of the generalized inertial number $I_m$ defined by Eq. (1) with $\alpha \approx 0.062$ and $\beta \approx 0.075$, for sets of parameters. The error bars represent the standard variation over steady-state flow, and the symbols and their colors are the same as those in Fig. 2. The solid lines are fitting forms similar to those used for $\mu/\mu_c$ and $\Phi/\Phi_c$ with different parameter values; see Supplementary Note 2.</figcaption></figure>

respect to microstructural variables. Let us consider the shear plane defined by the directions of shearing $\dot{y}$ and confining stress $\sigma_p$ (the flow being translationally invariant in the lateral direction). In this plane, the number density $E(\Theta)$ of bonds (both contacts and non-contact capillary bridges) per particle along a direction $\Theta$ can be approximated by$^{15,21,43}$

$$E(\Theta) = \frac{Z}{2\pi}\{1 + A \cos 2(\Theta - \Theta_b)\}, \quad (10)$$

which represents a truncated Fourier transform of $E(\Theta)$ with a period $\pi$ since no intrinsic polarity can be attributed to contact orientations. Higher-order terms are generally negligibly small in granular flows. The coordination number Z and bond orientation anisotropy A are the lowest-order descriptors of granular texture. The angle $\Theta_b$ is the privileged bond orientation, and its value is close to $\pi/4$ with respect to the flow direction.

Figure 5 shows Z and A normalized by their quasistatic values $Z_c$ and $A_c$, respectively, as a function of $I_m$. We see that, despite their larger variability as compared to $\mu(I_m)$ and $\Phi(I_m)$, the data points collapse on a master curve within our statistical precision as a function of $I_m$ with the same values $\alpha \approx 0.062$ and $\beta \approx 0.075$ as in $\mu(I_m)$ and $\Phi(I_m)$. The insets show the extent to which the same data points are dispersed when plotted as a function of I. We find that $Z_c$ is a decreasing linear function of $\xi$ (as $\Phi_c$) and $A_c$ is an increasing linear function of $\xi$ (as $\mu_c$); see Supplementary Fig. 3. Moreover, the functional forms that fit $A/A_c$ and $ZZ_c$ are the same as $G_\mu$ and $G_\Phi$ as a function of $I_m$ with different values of their free fitting parameters; see Supplementary Eqs. (3) and (4). This scaling of microstructural variables with $I_m$ through a functional dependence similar to flow variables clearly indicates that the shear strength is mainly due to the increasing aptitude of the particles to self-organize in an anisotropic network$^{41,44}$.

## Discussion

The additive rheology of granular materials, as demonstrated in this paper, is by no means self-evident. The expectation that a system involving several dimensionless parameters can ultimately be described by a single parameter combining those parameters is unusual. The deep reason behind such a behavior is the very nature of granular materials in which the particle interactions are concentrated at the contact points and the local dynamics is controlled by the shear rate. Hence, by careful distinction of stresses depending on the shear rate (first group) from those that are independent (second group), a single parameter can be defined by means of the stress additivity property. In this respect, the modified inertial number $I_m$ is a conceptual extension of the inertial number to arbitrary interactions between particles. Such a scaling works, however, by distinguishing the quasistatic limit and normalizing the packing fraction, apparent friction coefficient, coordination number and bond orientation anisotropy by their quasistatic limit values that depend on the dimensionless numbers of the second group ($\xi$ in our case). Encompassing dry and cohesive-viscous flows, quasistatic and dynamic states, flow variables and microstructural parameters, this scaling provides a general framework for complex granular flows.

This framework can be applied to quantify the effects of friction coefficient between particles and particle shape and size distributions as material parameters that can influence the relative roles of internal stresses in the collective flow behavior, and thus the fitting parameters. Previous results on the rheology of granular materials suggest that the functional forms fitting the master curves should not depend on the particle shapes, size distributions and friction coefficient, although their free parameters will certainly do$^{45-52}$. This is, however, a crucial step forward for application to different types of granular materials and for comparison with experiments.

It is also useful to point out that, for given particle size d and density $\rho_s$, each characteristic stress $\Sigma$ corresponds also to a characteristic time $T = d(\rho_s/\Sigma)^{1/2}$. This relation implies that, by virtue of stress additivity, the inverse quadratic times $1/T^2$ are additive. Hence, using the characteristic times and this rule, we arrive at the same expression of $I_m$ as in Eq. (1). Another interesting issue concerns an alternative rheology based on a ‘multiplicative’ expression of the flow variables, rather than additive combination of control parameters, as in ref. 34 applied to a cohesive-frictional granular material. In this approach, the apparent friction coefficient $\mu$ is expressed as a product of distinct functions of the primary dimensionless control parameters. This multiplicative partition works quite well for the cohesion index, and a prefactor similar to $\mu_0(1 + a\xi)$ in Eq. (4) is obtained. For the other functions, it is worth considering in more detail how they relate to the framework developed in this paper.

Let us finally recall that the inertial number was initially introduced in the context of cohesionless granular materials where the gravity or applied confining stress prevail. However, the cohesive stress may largely exceed the confining stress in fine powders, and therefore the inertial and viscous stresses must be advantageously compared to the cohesive stress rather than the confining stress. It is easy to see that, in this limit ($\sigma_p \to 0$), the modified inertial number is reduced to $I_m = \{Ca(\beta + St)/\alpha\}^{1/2}$. This is a simple expression that is expected to scale cohesive processes such as wet granulation and impact dynamics of cohesive aggregates. We thus propose to use impact experiments as a convenient means to investigate this scaling.

## Methods

**The model of capillary bridge.** We assume that the liquid inside the agglomerate is in the ‘pendular’ state with a uniform distribution of capillary bridges between particles$^{38,53-59}$. This distribution may be a consequence of mixing the liquid with
the particles, drainage of a saturated packing, or capillary condensation from a vapor. For a separation distance above a debonding distance $d_{rupt}$, the bridge breaks and its liquid is shared between the two particles proportionally to their sizes$^{53,60}$.

The capillary force $f_c$ between two particles depends on the liquid volume $V_b$ of the bond, liquid-vapor surface tension $\gamma_s$ and particle-liquid-gas contact angle $\theta$. We used the following expression$^{61}$:

$$f_c = \begin{cases} -\kappa R, & \text{for } \delta_n < 0, \\ -\kappa R e^{-\delta_n/\lambda}, & \text{for } 0 \le \delta_n \le d_{rupt}, \\ 0, & \text{for } \delta_n > d_{rupt}, \end{cases} \quad (11)$$

where $R = \sqrt{R_i R_j}$ is the geometrical mean radius and the pre-factor $\kappa$ is

$$\kappa = 2\pi\gamma_s \cos\theta. \quad (12)$$

The debonding distance $d_{rupt}$ is given by$^{53}$

$$d_{rupt} = \left(1 + \frac{\theta}{2}\right) V_b^{1/3}. \quad (13)$$

The characteristic length $\lambda$ in Eq. (11) is given by

$$\lambda = c h(r) \left(\frac{V_b}{R'}\right)^{1/2}, \quad (14)$$

where $R' = 2R_i R_j / (R_i + R_j)$ and $r = \max\{R_i/R_j; R_i/R_j\}$ are the harmonic mean radius and the size ratio between two particles, $h(r) = r^{-1/2}$, and $c \approx 0.9$. In all simulations, we set $\theta = 0$ with the assumption that the particles are covered by a layer of the wetting liquid.

Note that the simulated system is an idealized model of wet granular materials in the pendular state. Nevertheless, we believe that our results can be extended to higher amounts of liquid since the liquid can easily flow in an unsaturated material to wet larger particle areas with a lower Laplace pressure. This leads to a nearly constant cohesive stress as far as the material is not fully saturated$^{54,55}$. Hence, the leading effect of increased liquid volume is simply the increase of debonding distance when Eq. (11) is used.

**The normal viscous force.** The normal lubrication force $f_{vis}$ due to the effect of liquid bridges between two smooth spherical particles is given by$^{62–64}$

$$f_{vis} = \frac{3}{2} \pi R^2 \eta \frac{v_n}{\delta_n}, \quad (15)$$

where $\eta$ is the liquid viscosity and $v_n$ is the relative normal velocity, assumed to be positive when the gap $\delta_n$ is decreasing. This force diverges when the gap $\delta_n$ tends to zero. But for slightly rough particles, the characteristic size of the asperities allows for collision in finite time. Hence, we introduce a characteristic length $\delta_{n0}$ corresponding to the size of asperities so that the lubrication force for $\delta_n > 0$ is given by

$$f_{vis} = \frac{3}{2} \pi R^2 \eta \frac{v_n}{\delta_n + \delta_{n0}} \quad \text{for } \delta_n > 0 \quad (16)$$

For $\delta_n < 0$ (a contact between two particles), we assume that the lubrication force remains equal to its largest value:

$$f_{vis} = \frac{3}{2} \pi R^2 \eta \frac{v_n}{\delta_{n0}} \quad \text{for } \delta_n \le 0. \quad (17)$$

In our simulations, we set $\delta_{n0} = 5 \times 10^{-4} d_{min}$, where $d_{min}$ is the smallest particle diameter. This value is sufficiently small to allow the lubrication force to be effective without leading to its divergence at contact.

**Simulation method.** For the simulations, we used the molecular dynamics (MD) method with frictional contact interactions modeled by linear elastic repulsion along the normal direction and linear spring with a Coulomb threshold along the tangential direction, together with the previous approximate expressions of capillary force and viscous force acting between neighboring particles.

The particle displacements are calculated by step-wise resolution of Newton's second law:

$$m_i \frac{d^2 \mathbf{r}_i}{dt^2} = \sum_j \left[ (f_n^{ij} + f_c^{ij} + f_{vis}^{ij}) \mathbf{n}^{ij} + f_i^{ij} \mathbf{t}^{ij} \right], \quad (18) \\ \mathbf{I}_i \frac{d\omega_i}{dt} = \sum_j f_i^{ij} c^{ij} \times \mathbf{t}^{ij},$$

where particle $i$ is assumed to interact with its neighbors $j$ via normal contact forces $f_n$, tangential contact forces $f_t$, capillary forces $f_c$, and viscous forces $f_{vis}$. $\omega_i$ is the rotation vector of particle $i$, and $m_i$, $\mathbf{I}_i$, and $\mathbf{r}_i$ are its mass, inertia matrix, and position, respectively. $\mathbf{n}^{ij}$ denotes the unit vector perpendicular to the contact plane between the particles $i$ and $j$ and points from $j$ to $i$. $\mathbf{t}^{ij}$ is the unit vector in the contact plane pointing in the direction opposite to the relative tangential displacement of the two particles. $c^{ij}$ is the vector joining the center of particle $i$ to the contact point with particle $j$.

<table><caption>**Table 1 Constant simulation parameters.**</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Number of particles</td><td>$N_p$</td><td>19,628</td><td>μm</td></tr><tr><td>Smallest particle diameter</td><td>$d_{min}$</td><td>800</td><td>kg.m<sup>-3</sup></td></tr><tr><td>Particle density</td><td>$\rho_s$</td><td>2600</td><td>kg.m<sup>-3</sup></td></tr><tr><td>Friction coefficient</td><td>$\mu_s$</td><td>0.4</td><td></td></tr><tr><td>Normal stiffness</td><td>$k_n$</td><td>10<sup>6</sup></td><td>N/m</td></tr><tr><td>Tangential stiffness</td><td>$k_t$</td><td>8.10<sup>5</sup></td><td>N/m</td></tr><tr><td>Normal damping</td><td>$\gamma_n$</td><td>0.5</td><td>Ns/m</td></tr><tr><td>Tangential damping</td><td>$\gamma_t$</td><td>0.5</td><td>Ns/m</td></tr><tr><td>Contact angle</td><td>$\theta$</td><td>0</td><td>degree</td></tr><tr><td>Time step</td><td>$\delta t$</td><td>3.10<sup>-7</sup></td><td>s</td></tr></tbody></table>

The normal contact force $f_n$ is the sum of four contributions:

$$f_n = f_n^e + f_n^d + f_c + f_{vis}, \quad (19)$$

where $f_n^e$ is the elastic repulsion force, and $f_n^d$ is the normal damping force. The elastic force $f_n^e = -k_n \delta_n$ is a linear function of the normal elastic deflection $\delta_n$, where $k_n$ is the normal stiffness, and the damping force $f_n^d = \gamma_n \delta_n$ is proportional to the relative normal velocity $\delta_n$, where $\gamma_n$ is the normal viscous damping parameter. These elastic and damping forces occur only when two particles are in contact ($\delta_n < 0$). The tangential force $f_t$ is composed of an elastic force $f_t^e = -k_t \delta_t$ and a damping force $f_t^d = \gamma_t \delta_t$, where $k_t$ is the tangential stiffness, $\gamma_t$ is the tangential damping parameter, and $\delta_t$ and $\gamma_t$ are the tangential displacement and velocity, respectively. According to the Coulomb friction law, the tangential force is below $\mu_{f,n0}$, where $\mu$ is the friction coefficient$^{65–67}$:

$$f_t = - \min \left\{ -k_t \delta_t + \gamma_t \dot{\delta}_t, |\mu f_n| \right\} \operatorname{sgn}(\dot{\delta}_t). \quad (20)$$

The tangential lubrication force was neglected as it is one order of magnitude below the normal lubrication force.

The equations of motion were integrated by a step-wise velocity-Verlet aglrophim$^{68}$. The constant physical parameters were set to typical values of fine granular materials composed of hard particles. We used a weak size polydispersity with a uniform distribution of particle volumes and a ratio 2 between the largest and smallest particle diameters. All the constant physical and numerical parameter values are given in Table 1. Note that the relative pressure $p^* = \sigma_p / ((d)k_n)$, representing the ratio of contact deflection to particle diameter, is of the order of $10^{-5}$. Our simulations correspond therefore with a high precision to the ideal limit of perfectly rigid particles.

## Data availability

All relevant data are available upon request from the authors.

## Code availability

The simulation code is available at sourcesup.renater.fr/www/cfgd3d and upon reasonable request from the authors.

Received: 6 September 2019; Accepted: 27 February 2020; Published online: 19 March 2020

## References

1. Savage, S. B. & Hutter, K. The dynamics of avalanches of granular materials from initiation to runout. Part I: analysis. *Acta Mechanica* **86**, 201–223 (1991).

2. Campbell, C. S., Cleary, P. W. & Hopkins, M. Large-scale landslide simulations: global deformation, velocities and basal friction. *J. Geophys. Res. Solid Earth* **100**, 8267–8283 (1995).

3. Cruden, D. M. & Varnes, J. D. Chapter: Landslide types and processes. In *Landslides: Investigation and Mitigation* Vol. 247, 36–75 (National Research Council, Washington, D.C., 1996).

4. Tegzes, P., Vicsek, T. & Schiffer, P. Avalanche dynamics in wet granular materials. *Phys. Rev. Lett.* **89**, 094301 (2002).

5. Rematre A., Malet, J. P., Maquaire, O., Ancey, C. & Locat, J. Flow behaviour and runout modelling of a complex debris flow in a clay-shale basin. *Earth Surf. Proc. Landforms* **30**, 479–488 (2005).

6. Staron, L. & Hinch, E. The spreading of a granular mass: role of grain properties and initial conditions. *Granul. Matter* **9**, 205–217 (2007).
7. Mutabaruka, P., Delenne, J.-Y., Soga, K. & Radjai, F. Initiation of immersed granular avalanches. *Phys. Rev. E* **89**, 052203 (2014).

8. Delannay, R., Valance, A., Mangeney, A., Roche, O. & Richard, P. Granular and particle-laden flows: from laboratory experiments to field observations. *J. Phys. D* **50**, 053001 (2017).

9. Mikami, T., Kamiya, H. & Horio, M. Numerical simulation of cohesive powder behavior in a fluidized bed. *Chem. Eng. Sci.* **53**, 1927-1940 (1998).

10. Iveson, S., Beathe, J. & Page, N. The dynamic strength of partially saturated powder compacts: the effect of liquid properties. *Powder Technol.* **127**, 149-161 (2002).

11. Herminghaus, S. Dynamics of wet granular matter. *Adv. Phys.* **54**, 221-261 (2005).

12. Rognon, P., Roux, J.-N., Naaim, M. & Chevoir, F. Dense flows of cohesive granular materials. *J. Fluid Mech.* **596**, 21-47 (2008).

13. Rondet, E., Delalonde, M., Ruiz, T. & Desfoursb, J. P. Fractal formation description of agglomeration in low shear mixer. *Chem. Eng. J.* **164**, 376-382 (2010).

14. Jaeger, H. & Nagel, S. Granular solids, liquids and gases. *Rev. Mod. Phys.* **68**, 1259-1273 (1996).

15. Radjai, F., Roux, J.-N. & Daouadji, A. Modeling granular materials: century-long research across scales. *J. Eng. Mech.* **143**, 04017002 (2017).

16. Radjai, F., Wolf, D. E., Jean, M. & Moreau, J. Bimodal character of stress transmission in granular packings. *Phys. Rev. Lett.* **80**, 61-64 (1998).

17. Mueth, D. M. et al. Signatures of granular microstructure in dense shear flows. *Nature* **406**, 385-389 (2000).

18. Radjai, F. & Roux, S. Turbulentlike fluctuations in quasistatic flow of granular media. *Phys. Rev. Lett.* **89**, 064302 (2002).

19. Corwin, E. I., Jaeger, H. M. & Nagel, S. R. Structural signature of jamming in granular media. *Nature* **435**, 1075-1078 (2005).

20. Richard, P., Nicodemi, M., Delannay, R., Ribière, P. & Bideau, D. Slow relaxation and compaction of granular systems. *Nat. Mat.* **4**, 121-128 (2005).

21. Radjai, F. & Richefeu, V. Bond anisotropy and cohesion of wet granular materials. *Philos. T. R. Soc. A* **367**, 5123-5138 (2009).

22. Staron, L., Lagrée, P.-Y., Josserand, C. & Lhuillier, D. Flow and jamming of a two-dimensional granular bed: Toward a nonlocal rheology? *Phys. Fluids* **22**, 113303 (2010).

23. Miller, T., Rognon, P., Metzger, B. & Einav, I. Eddy viscosity in dense granular flows. *Phys. Rev. Lett.* **111**, 058002 (2013).

24. Luding, S. So much for the jamming point. *Nat. Phys.* **12**, 531-532 (2016).

25. Kou, B. et al. Granular materials flow like complex fluids. *Nature* **551**, 360-363 (2017).

26. GDR-MiDi. On dense granular flows. *Eur. Phys. J. E* **14**, 341-365 (2004).

27. da Cruz, F., Emam, S., Prochnow, M., Roux, J.-N. & Chevoir, F. Rheophysics of dense granular materials: discrete simulation of plane shear flows. *Phys. Rev. E* **72**, 021309 (2005).

28. Jop, P., Forterre, Y. & Pouliquen, O. A constitutive law for dense granular flows. *Nature* **441**, 727-730 (2006).

29. Boyer, F., Guazzelli, E. & Pouliquen, O. Unifying suspension and granular rheology. *Phys. Rev. Lett.* **107**, 188301 (2011).

30. Kamrin, K. & Koval, G. Nonlocal constitutive relation for steady granular flow. *Phys. Rev. Lett.* **108**, 178301 (2012).

31. Trulsson, M., Andreotti, B. & Claudin, P. Transition from the viscous to inertial regime in dense suspensions. *Phys. Rev. Lett.* **109**, 118305 (2012).

32. Berger, N., Azéma, E., Douce, J.-F. & Radjai, F. Scaling behaviour of cohesive granular flows. *Eur. Phys. Lett.* **112**, 64004 (2015).

33. Khamseh, S., Roux, J.-N. & Chevoir, FmC Flow of wet granular materials: a numerical study. *Phys. Rev. E* **92**, 022201 (2015).

34. Roy, S., Luding, S. & Weinhart, T. A general(ized) local rheology for wet granular materials. *New J. Phys.* **19**, 043014 (2017).

35. Amarsid, L. et al. Viscoinertial regime of immersed granular flows. *Phys. Rev. E* **96**, 012901 (2017).

36. Badetti, M., Fall, A., Chevoir, F. & Roux, J.-N. Shear strength of wet granular materials: macroscopic cohesion and effective stress. *Eur. Phys. J. E* **41**, 68 (2018).

37. Bagnold, R. A. Experiments on a gravity-free dispersion of large solid spheres in a newtonian fluid under shear. *Proc. Royal Soc. London* **225**, 49-63 (1954).

38. Richefeu, V., El Youssoufi, M. & F. Radjai, F. Shear strength properties of wet granular materials. *Phys. Rev. E* **73**, 051304 (2006).

39. Rognon, P. G., Roux, J.-N., Wolf, D., Naaim, M. & Chevoir, F. Rheophysics of cohesive granular materials. *Eur. Phys. Lett.* **74**, 644-650 (2006).

40. Pouliquen, O., Cassar, C., Jop, P., Forterre, Y. & Nicolas, M. Flow of dense granular material: towards simple constitutive laws. *J. Stat. Mech.* **2006**, 07020-0702 (2006).

41. Azéma, E. & Radjai, F. Internal structure of inertial granular flows. *Phys. Rev. Lett.* **112**, 078001 (2014).

42. Bocquet, L., Colin, A. & Ajdari, A. Kinetic theory of plastic flow in soft glassy materials. *Phys. Rev. Lett.* **103**, 036001 (2009).

43. Rothenburg, L. & Bathurst, R. J. Analytical study of induced anisotropy in idealized granular materials. *Géotechnique* **39**, 601-614 (1989).

44. Radjai, F., Delenne, J.-Y., Azéma, E. & Roux, S. Fabric evolution and accessible geometrical states in granular materials. *Granul. Matter* **14**, 259-264 (2012).

45. Voivret, C., Radjai, F., Delenne, J.-Y. & Youssoufi, M. S. E. Multiscale force networks in highly polydisperse granular media. *Phys. Rev. Lett.* **102**, 178001 (2009).

46. Azéma, E., Radjai, F. & Saussine, G. Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles. *Mech. Mat.* **41**, 721-741 (2009).

47. Estrada, N., Azéma, E., Radjai, F. & Taboada, A. Identification of rolling resistance as a shape parameter in sheared granular media. *Phys. Rev. E* **84**, 011306 (2011).

48. Cegeo, B. et al. Particle shape dependence in 2D granular media. *Eur. Phys. Lett.* **98**, 44008 (2012).

49. Shaebani, M. R., Madadi, M., Luding, S. & Wolf, D. E. Influence of polydispersity on micromechanics of granular materials. *Phys. Rev. E* **85**, 011301 (2012).

50. Ogarko, V. & Luding, S. Prediction of polydisperse hard-sphere mixture behavior using tridisperse systems. *Soft Matter* **9**, 9530-9534 (2013).

51. Kumar, N., Imole, O. I., Magnanimo, V. & Luding, S. Effects of polydispersity on the micro-macro behavior of granular assemblies under different deformation paths. *Particuology* **12**, 64-79 (2014).

52. Nguyen, D.-H., Azéma, E., Radjai, F. & Sornay, P. Effect of size polydispersity versus particle shape in dense granular media. *Phys. Rev. E* **90**, 012202 (2014).

53. Lian, G., Thornton, C. & Adams, M. A theoretical study of the liquid bridge forces between two rigid spherical bodies. *J. Colloid Interf. Sci.* **161**, 138-147 (1993).

54. Scheel, M. et al. Liquid distribution and cohesion in wet granular assemblies beyond the capillary bridge regime. *J. Phys. Cond. Matter* **20**, 494236 (2008).

55. Delenne, J.-Y., Richefeu, V. & Radjai, F. Liquid clustering and capillary pressure in granular media. *J. Fluid Mech.* **762**, R5 (2012).

56. Vo, T.-T., Mutabaruka, P., Delenne, J.-Y., Nezamabadi, S. & Radjai, F. Strength of wet agglomerates of spherical particles: effects of friction and size distribution. *EPJ Web Conf.* **140**, 08021 (2017).

57. Than, V. D. et al. Basic mechanical properties of wet granular materials: a DEM study. *J. Eng. Mech.* **143**, C4016001 (2017).

58. Vo, T.-T. et al. Mechanical strength of wet particle agglomerates. *Mech. Res. Com.* **92**, 1-7 (2018).

59. Vo, T.-T. Agglomeration of wet particles in dense granular flows. *Eur. Phys. J. E* **42**, 127 (2019).

60. Willett, C., Adans, M., Johnson, S. & Seville, J. Capillary bridges between two spherical bodies. *Langmuir* **16**, 9396-9405 (2000).

61. Richefeu, V., Radjai, F. & Youssoufi, M. S. E. Stress transmission in wet granular materials. *Eur. Phys. J. E* **21**, 359-369 (2007).

62. Brenner, H. The slow motion of a sphere through a viscous fluid towards a plane surface. *Chem. Eng. Sci.* **16**, 242-251 (1961).

63. Happel, J. & Brenner, H. *Low Reynolds Number Hydrodynamics* (Martinus Nijhoff Publishers, Dordrecht, 1983).

64. Lefebvre, G. & Jop, P. Erosion dynamics of a wet granular medium. *Phys. Rev. E* **8**, 032205 (2013).

65. Schäfer, J., Dippel, S. & Wolf, D. E. Force schemes in simulations of granular materials. *J. Phys. I France* **6**, 5-20 (1996).

66. Dippel, S., Batrouni, G. G. & Wolf, D. E. How transversal fluctuations affect the friction of a particle on a rough incline. *Phys. Rev. E* **56**, 3645-3656 (1997).

67. Luding, S. Collisions and contacts between two particles. (eds. Herrmann, H. J., Hovi, J.-P. & Luding, S.) In *Physics of Dry Granular Media-NATO ASI Series E350*, 285 (Kluwer Academic Publishers, Dordrecht, 1998).

68. Radjai, F. & Dubois, F. (eds.) *Discrete Numerical Modeling of Granular Materials* (Wiley-ISTE, New-York, 2011).

## Acknowledgements

The authors acknowledge financial support by the Ministry of Education and Training in Vietnam, Campus France, and Labex NUMEV (University of Montpellier). This work was realized with the support of MESO@LR-Platform at the University of Montpellier.

## Author contributions

F.R., S.N. and J.Y.D. conceived the research idea. P.M. developed the particle dynamics program CFGD++3D used in this work and provided support for the simulations. T.T.V. carried out all the simulations, analyzed the data, and prepared all the figures. All authors discussed the results and contributed to the writing of the paper.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41467-020-15263-3.

**Correspondence** and requests for materials should be addressed to F.R.

**Peer review information** *Nature Communications* thanks Stefan Luding, M. Saadatfar and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020