Rheology of inelastic hard spheres at finite density and shear rate

Considering a granular fluid of inelastic smooth hard spheres we discuss the conditions delineating the rheological regimes comprising Newtonian, Bagnoldian, shear thinning, and shear thickening behavior. Developing a kinetic theory, valid at finite shear rates and densities around the glass transition density, we predict the viscosity and Bagnold coefficient at practically relevant values of the control parameters. The determination of full flow curves relating the shear stress $\sigma$ to the shear rate $\dot\gamma$, and predictions of the yield stress complete our discussion of granular rheology derived from first principles.

Predicting and understanding granular flow is desirable for safety and efficiency. Many geophysical flows from avalanches to landslides involve macroscopic particles and potentially threaten lives all around the planet [1,2]. A large fraction of the raw materials handled in industry comes in granular form [3]. With the advent of 3d printing technologies, this fraction will further increase [4,5]. Hence, efficient handling of granular flows promises considerable advantages like energy savings [6]. The crucial property of granular particles-namely that they are of macroscopic size-makes a theoretical description challenging [7][8][9]. Firstly, dissipative collisions break time-reversal symmetry and place granular flows firmly in the realm of far-from-equilibrium physics [10,11]. Secondly, the macroscopic mass of an individual granular particle makes its thermal excitations negligible [7] and necessitates a driving force that is continuously acting on the granular assembly to keep it flowing [12][13][14].
At small volume fractions ϕ 1 and infinitesimal shear ratesγ → 0 standard procedures starting from the Boltzmann or Enskog equation predict a Newtonian rheology for (smooth) granular particles. Various approximate expressions exist for the viscosity as a function of inelasticity (often quantified by a coefficient of restitution ε) and volume fraction [15][16][17][18][19][20]. However, most natural and industrial granular flows occur at considerable volume fractions all the way up to close packing densities. Relevant shear ratesγ are also often substantial compared to microscopic time scales.
Experimental as well as numerical studies have found a wealth of phenomena at finite densities and shear rates [21][22][23][24][25][26]. One of the earliest results of granular physics by Bagnold [21]-now commonly referred to as Bagnold scaling-is the observation that granular fluids do not follow a Newtonian rheology, but that the shear stress shows a quadratic dependence on the shear rate, σ = Bγ 2 , instead. Phenomenologically, Bagnold scaling is a manifestation of shear thickening as the shear rate dependent viscosity η(γ) ≡ σ/γ ∼γ increases with shear rate. Theoretical predictions of the corresponding Bagnold coefficient B are rare [15,27,28] and for low density or low shear rate FIG. 1. (a) Theoretical dynamic state diagram of a granular fluid with coefficient of restitution ε = 0.5 depending on packing fraction ϕ (φ := ϕ + 5.5%) and shear rateγ. Rheological behavior (R) color coded as indicated. Critical Péclet numbers Pe * (dashed) where shear heating becomes important and Peα (dash-dotted), whereγ, matches the intrinsic relaxation rate, τ −1 α . (b) Schematic of the shear geometry with the shear profile overlayed (red). (c) Maximal Péclet number, Pe∞, (orange dashed, left axis), and Bagnold coefficient, B, (blue solid, right axis) as function of volume fraction ϕ comparing Bagnold's measurements [21] (symbols) and the kinetic theory presented here (lines).
only. Later studies also show the opposite behavior, i.e., shear thinning, in granular fluids. Taken together, these observations imply that the Newtonian rheology is valid only in a limited part of the parameter space and capturing the full rheology in a single theoretical framework remains a challenge.
In this letter we will show that, indeed, granular flows display all three regimes: Newtonian, shear thinning, and arXiv:1710.04452v2 [cond-mat.soft] 11 Jul 2018 shear thickening (Fig. 1a). We will derive the conditions to observe any of these behaviors. Based on this classification, we then calculate the relevant material properties, namely the Newtonian viscosity η, the Bagnold coefficient B, and more generally, flow curves, i.e., σ(γ), or, equivalently, η(γ).
Formally, the stationary shear stress σ at a constant finite shear rate can be expressed via a Generalized Green-Kubo Relation, as the time integral over the shear stress auto-correlation whose evolution depends on the flow [29]. Here σ xy is the Kirkwood shear stress defined in terms of the particle positions and momenta [30], V is the volume and T the temperature of a reference state. The average · ref is performed with respect to the unsheared reference system. Although the relation was originally derived for a reference system in thermal equilibrium [29] it ultimately relates stationary states including out-of-equilibrium reference systems [31]. Our choice allows to specify the perturbing stress in Eq. (1) microscopically. The auto-correlation function can be controlled by identifying a dominant decay channel. For thermal colloidal suspensions the slow structural relaxations close to the glass transition provides such a clearly defined, slow decay [32]. Approximating the stress auto-correlation function in Eq. (1) in terms of the density correlator Φ q (t) lies at the heart of the Integration Through Transients (ITT) formalism [29,33], Here q(−t) ≡ (1 + kt) · q denotes the advected wave vector and the velocity gradient tensor k αβ :=γδ αy δ βx prescribes simple shear (Fig. 1b). Note that the advected wave vector's dependence onγ effects a nonlinear stressshear rate relation in Eq. (2). The coupling constant V qq(−t) can be calculated explicitly [see Refs. [31,33] and Eq. (4) below]. The ITT framework has led to a wealth of qualitative and quantitative predictions regarding the rheology of thermalized colloidal suspensions [34][35][36][37] Power Balance-The temperature T of an overdamped colloidal suspension is assumed to be controlled by a heat bath. In particular, one assumes that the work performed on the suspension by the shear force does not increase the temperature but that instead the temperature can be chosen freely. Formulating ITT for underdamped, Newtonian dynamics, a thermostat has to be included explicitly with the final results depending on the precise choice of artificial thermostat [38,39]. In granular flows the balance between shear heating and dissipation occurs naturally and actually controls the qualitative behavior of the granular fluid. The system we have in mind in the following is a sheared fluidized bed [40]. To be specific, let's consider a fluid composed of monodisperse smooth hard spheres of diameter d, and mass m = 1 with a coefficient of normal restitution ε. We assume that, initially, the fluid is prepared at a given density n, or packing fraction ϕ = πnd 3 /6 and a random fluidization force is applied throughout the system with a characteristic power per particle P D to mimic fluidization [41] [42]. Then the initial granular temperature T 0 results from the power balance √ T , increases with density and temperature and χ(ϕ) denotes the pair correlation function at contact [20].
Once shear is applied with a prescribed shear rateγ, shear heating has to be included in the power balance of the steady state, resulting in a higher temperature T > T 0 (Protocol H in Ref. [31]). This clearly defines two regimes: A fluidization dominated regime, where σγ nP D and a shear dominated regime, where σγ nP D including purely shear driven systems (P D ≡ 0). Choosing the particle diameter, d, as our length scale and the inverse collision frequency in the stationary state, ω −1 c (T ), as our time scale the packing fraction, the coefficient of restitution, and the Péclet number, Pe =γ/ω c (T ), alone determine the system's state. The temperature only controls the overall timescale. The shear stress scales as σ = nTσ with a dimensionless functionσ =σ(Pe, ϕ, ) because nT provides the only scale for an energy density. The crossover between the regimes is expected where σγ ∼ nP D . This implicitly determines a critical Péclet number Pe * (ϕ, ε) (cf. Tab. I).
In the fluidization dominated regime, the temperature remains independent of the shear rate such that we can linearizeσ ∼γ to find Newtonian behavior. In the shear dominated regime the power balance reads nTσγ = nT Γω c [43], i.e., the collision frequency is proportional to the shear rate and the corresponding Péclet number, Pe = Pe ∞ (ϕ, ε) > 0 (cf. Tab. I), is independent of the shear rate. Fromγ ∝ ω c ∝ √ T we obtain T ∝γ 2 and with that, σ ∝ T ∝γ 2 : Bagnold rheology is observed in the shear dominated regime where the shear rate is the only time scale and the Péclet number, Pe ∞ , is fixed by the packing fraction and the coefficient of restitution. Upon increasing the shear rateγ towards the Bagnold regime, the temperature T (Pe) = T 0 · (1 − Pe / Pe ∞ ) −2/3 diverges. This implies that Pe ∞ is the maximal Péclet number which can not be exceeded in a granular fluid. Glassy dynamics and yield stress-Granular fluids have been found to undergo a glass transition [14] at a critical packing fraction ϕ c (ε) which increases with increasing dissipation [44]. Upon approaching the (granular) glass transition, ϕ ϕ c , the characteristic correlation time for density fluctuations, τ α , diverges. Forγ τ −1 α the rheology is Newtonian, as the granular fluid can respond immediately to the slow shear deformations. However, the viscosity diverges as ϕ ϕ c . For higher prescribed shear rates,γ τ −1 α [i.e., Pe Pe α (ϕ, ε), cf. Tab. I], the glass is forcibly molten.
To lowest order it can be assumed that Φ q (t → ∞) ∝ e −γt for Pe > Pe α . Consequently the Green-Kubo-Integral, Eq. (2), becomes independent of the shear rate, i.e, σ(γ) ≈ const. In terms of the viscosity η(γ) ∼ 1/γ we expect to observe shear thinning. For even higher shear rates, Pe > Pe * , shear heating will dominate and eventually bend the flow curve to the shear thickening Bagnold regime. For densities above the glass transition, ϕ > ϕ c , Pe α → 0 and the Newtonian regime vanishes altogether. Instead, a finite (dynamical) yield stress σ y := σ(γ → 0) > 0 emerges which has to be overcome to melt the amorphous granular glass.
Slow shear -For small shear rates Pe → 0 in the linear response regime, the glass transition divides the state space into two qualitatively different phases. For relatively low densities, ϕ < ϕ c , gITT provides corrections of order ϕ 2 to the low density, Enskog predictions [19] for the viscosity. For densities approaching the glass transition, the divergence of the relaxation time, τ α , results in a divergent viscosity, η(ϕ . The critical exponent γ(ε) ∼ 2.5-2.3 weakly decreases with increasing dissipation [46] and compares well with the experimental value γ ≈ 2.35 [48]. Experimental values for ϕ c (ε) [49][50][51] are compatible with our interpretation as a granular glass transition as well. At the glass transition, ϕ ≡ ϕ c (ε), we find a finite critical yield stress (cf. Fig. 3d), σ c y , in the range 6nT -9nT , compatible with experiments [50], and comparable to theoretical predictions (σ c y /nT ∼ 6 [52]) and measurements (σ c y /nT ∼ 10-15 [53]) for colloidal suspensions. Thereby placing granular fluids firmly in the realm of soft matter.  B (cf. Fig. 3c). The latter increases with density and elasticity (larger ε) as both trends makes the temperature more sensitive to changes in the shear rate. In particular, B diverges as ε → 1.
As a first quantitative application of gITT, let us compare our predictions to Bagnold's original data [21]. To this end, we extract the Bagnold coefficient, B, and the maximal Péclet number, Pe ∞ , from his measurements (Fig. 1c) [54]. Our kinetic theory proves to recover the qualitative trends of both B and Pe ∞ as a function of packing fraction with no adjustable parameters. Considering that Bagnold's measurements have been shown to leave room for improvement [55] the prediction of gITT compare favorably. Qualitatively the strongly agitated flow curves of Refs. [56,57] agree well with our discussion showing all three regimes. However, note that the packing fraction in the shear zone is unknown and uncontrolled in these experiments in contrast to what is assumed here. More tailored and careful measurements are needed to assess gITT's quantitative accuracy.
Depending on experimental conditions additional effects may become relevant that go beyond the model considered in this letter. On earth (but not, e.g., on the moon [58]) all granular flows are actually two-phase flows of granular particles together with an interstitial fluid (mostly air or water). At sufficiently low packing fraction, the effective viscosity of the molecular fluid becomes relevant [59]. This will introduce another Newtonian regime at small shear rates that crosses over to Bagnold rheology when the stress induced by the granular particles dominates over the viscous stress of the interstitial fluid. For high shear rates, the Bagnold regime will obviously not extend toγ → ∞. At some point typical inter-particle forces are so high that interactions can no longer be regarded as hard-core. A finite interaction time compared to the shear rate appears as a new time scale and destroys Bagnold scaling [60]. For very high densities, approaching random close packing, the rheology will be dominated by the imminent jamming transition which is not accounted for in the present model. We prescribe a linear shear profile and a homogeneous constant density and therefore necessarily obtain monotonous flow curves. Non-monotonic, unsta-ble flow curves and the associated discontinuous shear thickening are possible in inhomogeneous systems only [26,56,57,61].
Conclusion-To summarize, we have discussed that the rich rheology of a granular fluid is controlled by three critical Péclet numbers (Tab. I): (i) Newtonian rheology prevails for Pe < Pe α . (ii) For intermediate shear rates, Pe α < Pe < Pe * , shear thinning reflects that the shear rate is faster than the intrinsic relaxation rate of the fluid which eventually results in a finite dynamic yield stress above the glass transition density. The latter can be substantially lower than the jamming transition commonly located at random close packing ϕ rcp ≈ 0.64 [62]. For low densities Pe α ∼ 1 and at the same time Pe * 1 in the elastic limit ε → 1. Under these conditions, Pe α > Pe * and the thinning regime vanishes altogether. (iii) For Pe * < Pe < Pe ∞ strong shear heating leads to shear thickening behavior which ultimately entails Bagnold scaling as the Péclet number approaches its maximum, Pe → Pe ∞ . This constitutes yet another shear thickening mechanism different from other mechanisms discussed in the literature, namely clustering [63,64], dilation [65,66], friction [67,68], or steric effects [69]. The fact that Pe ∞ is finite implies that a kinetic theory to predict the Bagnold coefficient B must be applicable at finite shear rates.
To support our arguments and to make them quantitative, we presented a kinetic theory based on the ITT formalism that recovers the phenomenology, covering many orders of magnitude in shear rate and shear stress, or viscosity, respectively. Earlier attempts at formulating ITT for inelastic soft spheres [70] retain no dissipative effects on the same level of approximation. For the inelastic hard sphere fluid considered here, besides the implicit dependence of Φ k (t) on the coefficient of restitution ε [46], Eq. (4) also includes an explicit dependence on ε. Thereby we extended quantitative predictions for the transport coefficients of a sheared granular fluid beyond the low density and low shear rate regime amenable to standard kinetic theories.
The results presented here will be useful as constitutive equations for modeling and simulating large scale granular flows which demand a continuum description. In addition, the Bagnold coefficient is needed for a recent kinetic theory [71]. We also hope that the availability of a kinetic theory for granular shear flow in a range of practically relevant parameters will spur quantitative experiments.
We acknowledge illuminating discussions with Hisao Hayakawa, Koshiro Suzuki, Thomas Voigtmann, and Claus Heussinger. We thank Philip Born for critically reading the manuscript and the DFG for partial funding through FOR 1394 and KR 4867/2-1.
We read off shear stress and pressure data from Figs. 3 and 4 in Ref. [21]. Calculation of the apparent Bagnold coefficient is straight forward, B app (γ) = σ/γ 2 . Fig. 4a shows that B app (γ) converges exponentially to what we take as the proper Bagnold coefficient B ≡ B app (γ → ∞). The particle mass m ≈ 1.2 mg we obtain from the particle diameter d = 1.32 mm and the fact that they are density matched with water. Bagnold reports packing densities in terms of the linear density λ which can directly be converted into a packing fraction [55]. Bagnold did not measure the coefficient of restitution of his particles so we assume ε = 0 [76].
To estimate the Péclet number from the shear rate, we need the collision frequency or the granular temperature. Neither of which have been measured by Bagnold. We estimate the temperature from the pressure by inverting the Woodcock equation of state [74]. The result is shown in Fig. 4b. We observe that we reasonably recover Bagnold scaling T ∝γ 2 and that the highest packing fraction seems to behave differently from the rest [55]. Although the publication has a different focus, Bagnold coefficients may be extracted from Fig. 2a in Fall et al. [59]. From the given particle diameter d ≈ 40 µm and the polystyrene particle's mass density ρ s = 1.05 g · cm 3 , we can derive the appropriately non-dimensionalized apparent Bagnold coefficient B app . The coefficient of restitution ε is not reported. From the measurements of much larger polystyrene spheres at higher impact velocities [75], we infer that the smaller particles at typically slower impact speeds have a coefficient of restitution ε 0.95.
Applying the same exponential fit as for Bagnold's data, we find that only the lowest three densities allow to extract a reliable value for the asymptotic Bagnold coefficient B. Treating ε as a fit parameter, good agreement between theory and experiment for ε = 0.97 can be obtained (see Tab. II). Note that the nearly elastic particles of Fall et al. result in much higher Bagnold coefficients than the strongly inelastic particles employed by Bagnold as expected. Exp. gITT 0.568 16 · 10 3 25 · 10 3 0.576 25 · 10 3 32 · 10 3 0.581 42 · 10 3 36 · 10 3