Additive rheology of complex granular flows

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.

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 σ p , or by collective particle motions during flow involving an inertial (or kinetic) pressure σ i $ ρ s hdi 2 _ γ 2 , where ρ s is the particle density, 〈d〉 is the mean particle diameter, and _ γ is the shear rate 27,37 . Hence, in the NPT statistical ensemble (with temperature T = 0 for a granular material), the apparent friction coefficient μ = σ t ∕σ p , where σ t is the shear stress, and the packing fraction Φ are expected to be uniquely dependent on the ratio I 2 ≡ σ i ∕σ 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 hdiðρ s =σ p Þ 1=2 under load vs. shear time _ γ À1 ), and it was found to unify experimental and numerical data in different flow geometries 26 .
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 σ c~fc ∕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 32,38,39 . As the stresses are additive, one may thus take the cohesive stress into account on the same footing as the confining pressure by replacing σ p by a linear combination σ n = σ p + ασ c , and therefore replacing I 2 by I 2 c σ i =σ n ¼ I 2 =ð1 þ αξÞ, where ξ = σ c ∕σ p is the cohesion index, and α is a material-dependent parameter 32 . In the same way, in submerged granular flows, where the viscous drag force σ v is present instead of cohesive stress, σ i may be replaced by a linear combination σ i + βσ v , leading to a visco-inertial number I 2 v ðσ i þ βσ v Þ=σ p ¼ I 2 ð1 þ β=StÞ, where St ≡ σ i ∕σ v is the Stokes number. In dense non-Newtonian suspensions, I v is found to be the control parameter for both μ and Φ 31,35 .
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 η, surface tension γ s , confining stress σ p and shear rate _ γ. 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.  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 σ p and moving at constant horizontal velocity v, and with periodic boundary conditions along lateral directions, The particle colors are proportional to their diameters. 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. 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 μ and Φ 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 The values of μ and Φ are normalized by μ c and Φ c , respectively, which design their quasi-static values (I m → 0) and vary linearly with ξ, as shown in the two insets to Fig. 3. The parameter values α ≃ 0.062 and β ≃ 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 α and β 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 σ p , inertial stress σ i , viscous stress σ v , and capillary stress σ c . The key variable for inertial flows is the shear rate. We thus distinguish the characteristic stresses that depend on the shear rate, i.e., σ i and σ v , from those that are independent of the shear rate, i.e., σ p and σ c . By virtue of stress additivity, the total shear-dependent stress is a linear combination σ i + βσ v of the former, and the total shearindependent stress is a linear combination σ p + ασ c of the latter. Hence, the flow variables (apparent friction coefficient and packing fraction) are expected to depend on the ratio (σ i + βσ v ) ∕ (σ p + ασ 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 These primary dimensionless parameters can be evaluated from the system control parameters. The order of magnitude of the viscous stress σ v is conveniently evaluated by replacing the average relative velocity ϑ $ _ γhdi induced by shearing in Eq. (15) and considering the dissipated power per unit volume f v δ rupt ∕ 〈d〉 3 , where δ rupt is the debonding distance, yielding σ v $ η_ γ. The capillary stress is of the order of the capillary force at contact (δ = 0) in Eq. (11) divided by the cross section 〈d〉 2 : σ c~γs ∕ d〉. Hence, the dimensionless parameters are given by Note that all other dimensionless variables can be expressed in terms of these three independent parameters. For example, the capillary number is given by Fitting functions for apparent friction and packing fraction. According to Eq. (1), I m → 0 only if I → 0. In this quasistatic limit, the flow variables may still depend on ξ, which is the only dimensionless variable in the absence of inertial and viscous stresses. Accordingly, the rheology is expected to be described by   6) and (7). The insets show the evolution of quasistatic values μ c and Φ c of the apparent friction coefficient and packing fraction, respectively, with the cohesion index ξ, and their linear fits (blue-solid lines).
the following general equations: These equations are based on functional distinction between the quasistatic limit (I m → 0) and shear-rate dependent behavior through I m . When I m → 0, we have μ → μ c and Φ → Φ c , and thus G μ → 1 and G Φ → 1. According to the simulation data displayed in Fig. 3, μ c and Φ c are linear functions of ξ: with a ≃ 0.095 and b ≃ 0.005. The limit values μ 0 ≃ 0.392 and Φ 0 ≃ 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 Φ 0 ≃ 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: with Δ μ ≃ 1.100, I μ ≃ 0.095 and I Φ ≃ 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 μ , I Φ and Δ μ , our simulation data are consistent with the experimental measurements of ref. 40 and ref. 28 in the dry case. The values of I μ and I Φ 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 α and β reflect the relative roles of viscous, inertial and cohesive forces in collective dissipation mechanisms whereas the values of I μ , I Φ , and Δ μ 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 ξ, G μ declines with increasing ξ (dynamic effect) whereas μ 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 ξ and I m . These roles are inversed for the packing fraction: Φ c declines whereas G Φ increases when ξ 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 η n and η t defined by σ n ¼ η n _ γ and σ t ¼ η t _ γ, respectively, where σ t = μσ n is the shear stress 29 . In this ensemble, the packing fraction Φ replaces pressure σ n as control parameter, and the rheology is characterized by the functions η n (Φ) and η t (Φ) 29 . This is the approach mostly used in experiments on suspensions. Although our simulations were carried out under NPT conditions, we may deduce η n (Φ) and η t (Φ) in the NVT ensemble from μ(I m ) and Φ(I m ). Since no external stress is imposed in NVT, the shear stress σ t is a dynamic variable that should scale with the internal shear-dependent stress σ i + βσ v . Moreover, the NPT and NVT points of view should be compared at the same normal stress state, i.e., σ n = σ p + ασ c . Hence, according to Eq. (1), σ n ¼ c n ðσ i þ βσ v Þ c n σ n I 2 m ¼ η n _ γ, implying c n ¼ 1=I 2 m ¼ η n =ðβη þ ρ s hdi 2 _ γÞ, and σ t = c t (σ i + βσ v ) with σ t ¼ μ=I 2 m ¼ η t =ðβη þ ρ s hdi 2 _ γÞ. In this way, in a volumecontrolled flow, c n and c t represent dimensionless viscosities with βη þ ρ s hdi 2 _ γ as reference viscosity; see Supplementary Note 3 for more detail. Figure 4 displays the effective dimensionless viscosities as a function of Φ. We see that all the data points collapse on a master curve when Φ is normalized by the critical packing fraction Φ c . Both viscosities diverge as Φ → Φ c and they are nicely fitted by the analytic expressions ; ð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 2 m =μ 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 respect to microstructural variables. Let us consider the shear plane defined by the directions of shearing _ γ and confining stress σ p (the flow being translationally invariant in the lateral direction). In this plane, the number density E(Θ) of bonds (both contacts and non-contact capillary bridges) per particle along a direction Θ can be approximated by 15,21,43 EðΘÞ which represents a truncated Fourier transform of E(Θ) with a period π 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 Θ b is the privileged bond orientation, and its value is close to π∕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 μ(I m ) and Φ(I m ), the data points collapse on a master curve within our statistical precision as a function of I m with the same values α ≃ 0.062 and β ≃ 0.075 as in μ(I m ) and Φ(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 ξ (as Φ c ) and A c is an increasing linear function of ξ (as μ c ); see Supplementary  Fig. 3. Moreover, the functional forms that fit A∕A c and Z∕Z c are the same as G μ and G Φ 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 (ξ 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][46][47][48][49][50][51][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 ρ s , each characteristic stress Σ corresponds also to a characteristic time T ¼ dðρ s =ΣÞ 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 μ 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 μ 0 (1 + aξ) 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 (σ p → 0), the modified inertial number is reduced to I m = {Ca(β + St)∕α} 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 γ s and particle-liquid-gas contact angle θ. We used the following expression 61 : for δ n <0; Àκ R e Àδ n =λ ; for 0 ≤ δ n ≤ d rupt ; 0; for δ n >d rupt ; is the geometrical mean radius and the pre-factor κ is The debonding distance d rupt is given by 53 The characteristic length λ in Eq. (11) is given by where R 0 ¼ 2R i R j =ðR i þ R j Þ and r ¼ maxfR i =R j ; R j =R i g are the harmonic mean radius and the size ratio between two particles, h(r) = r −1∕2 , and c ≃ 0.9. In all simulations, we set θ = 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][63][64] where η is the liquid viscosity and v n is the relative normal velocity, assumed to be positive when the gap δ n is decreasing. This force diverges when the gap δ 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 δ n0 corresponding to the size of asperities so that the lubrication force for δ n > 0 is given by For δ n < 0 (a contact between two particles), we assume that the lubrication force remains equal to its largest value: In our simulations, we set δ n0 = 5 × 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: 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 . ω i is the rotation vector of particle i, and m i , I i , and r i are its mass, inertia matrix, and position, respectively. n ij denotes the unit vector perpendicular to the contact plane between the particles i and j and points from j to i. 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.
The normal contact force f n is the sum of four contributions: where f e n is the elastic repulsion force, and f d n is the normal damping force. The elastic force f e n ¼ Àk n δ n is a linear function of the normal elastic deflection δ n , where k n is the normal stiffness, and the damping force f d n ¼ γ n _ δ n is proportional to the relative normal velocity _ δ n , where γ n is the normal viscous damping parameter. These elastic and damping forces occur only when two particles are in contact (δ n < 0). The tangential force f t is composed of an elastic force f e t ¼ Àk t δ t and a damping force where k t is the tangential stiffness, γ t is the tangential damping parameter, and δ t and _ δ t are the tangential displacement and velocity, respectively. According to the Coulomb friction law, the tangential force is below μf n , where μ is the friction coefficient 65-67 : 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 aglorithm 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 * = σ 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.