Power-law scaling in granular rheology across flow geometries
PPower-law scaling in granular rheology across flow geometries
Seongmin Kim and Ken Kamrin Harvard John A. Paulson School of Engineering and Applied Sciences,Harvard University, Cambridge, Massachusetts 02138, USA Department of Mechanical Engineering, MIT, Cambridge, Massachusetts 02139, USA (Dated: May 25, 2020)Based on discrete element method simulations, we propose a new form of the constitution equationfor granular flows independent of packing fraction. Rescaling the stress ratio µ by a power ofdimensionless temperature Θ makes the data from a wide set of flow geometries collapse to a mastercurve depending only on the inertial number I . The basic power-law structure appears robust tovarying particle properties (e.g. surface friction) in both 2D and 3D systems. We show how thisrheology fits and extends frameworks such as kinetic theory and the Nonlocal Granular Fluiditymodel. Granular materials exhibit complex mechanical behav-iors: depending on the situation, they can either sustainloads like solids or flow like fluids. Diverse attempts havebeen made to build a continuum model for granular flows.The µ ( I ) rheology [1–3], a phenomenological model, sug-gests a one-to-one relation between two local dimension-less variables, the shear-to-normal stress ratio µ = τ /P and the inertial number I ≡ ˙ γ/ (cid:112) P/ρ s d for 3D spheresand I ≡ ˙ γ/ (cid:112) P/m for 2D disks where τ is the shearstress, P is the pressure, ˙ γ is the shear rate, ρ s is theparticle density, and d and m are respectively the meanparticle diameter and mass. In this model, the shear ratevanishes if µ is smaller than a bulk friction coefficient µ s ,and I monotonically increases as µ increases for µ > µ s .However, this one-to-one relation between µ and I loses accuracy outside of homogeneous shear flows. Ingeneral, nonlocal phenomena deviate flow from the µ ( I )rheology [1, 4–12]. To reconcile this deviation, nonlocalmodels such as the Nonlocal Granular Fluidity (NGF)model [13–16] have been proposed. Inspired by a non-local model for emulsion flows [17, 18], the NGF modelassumes that a scalar “fluidity” field g enters the flow rulethrough ˙ γ = gµ and follows a phenomenological reaction-diffusion differential equation where the fluidity is gener-ated by shearing and diffuses in space. Recently, through3D discrete element method (DEM) simulations, Zhangand Kamrin [19] have found that the fluidity field canbe represented kinematically by the velocity fluctuations δv and the packing fraction φ : g = ˙ γ/µ = F ( φ ) δv/d .Since g , the single evolving state field of the NGF model,appears to arise from two kinematically observable statefields, we are motivated to seek further possible reduc-tions.Interestingly, kinetic theory, which mathematically de-rives the constitutive equations using the Chapman-Enskog method, predicts a similar relation between µ ,˙ γ , δv , and φ . Introducing a granular temperature T ≡ δv /D where D is the spatial dimensions, kinetic theorypredicts the pressure as P = ρ s F ( φ ) T and the shearstress as τ = ρ s F ( φ ) √ T ˙ γ d where F ( φ ) and F ( φ ) de-pend on the radial distribution function [20–23]. Thus,kinetic theory asserts µ = ( F ( φ ) /F ( φ )) ˙ γd/ √ T which becomes identical to Zhang’s relation if F ( φ ) /F ( φ ) = √ F ( φ ). According to kinetic theory, since φ can be sub-stituted by a function of dimensionless granular tempera-ture Θ ≡ ρ s T /P , µ/I should be expressible as a functionof Θ. Although the assumptions of standard kinetic the-ory become less accurate near the jammed state, we areintrigued to consider whether some generic µ ( I, Θ) rela-tion continues to exists into the dense regime, effectivelyremoving rheological dependence on φ . The notion ofexpanding the µ ( I ) model by dimensionless temperaturehas also been considered in [24], which we shall discusslater.To explore a potential µ ( I, Θ) relation, we take ahint from the power-law dependencies of thermodynamicquantities in many complex systems which exhibit con-tinuous phase transitions. Near the critical tempera-ture T c where the microscopic entities are highly corre-lated, the macroscopic fields follow scaling forms charac-terized by a power function of the reduced temperature( T − T c ) /T c [25]. Although granular systems are ather-mal, the velocity fluctuations created by shearing may actlike the temperature. Moreover, previous studies haveobserved more correlated motion of grains as a granularmaterial approaches the jammed state [26–29]. It is thusnatural to suspect power-law scaling in a µ ( I, Θ) relationas a possible unifying principle in granular rheology.Inspired by critical scaling, in this Letter we show thatrescaling µ by a simple power of Θ collapses data frommany DEM strongly onto a master curve that dependsonly on I . In doing so, we identify and validate a generalrelation of the form µ ( I, Θ) that holds across geometriesand flow regimes.We use LAMMPS to simulate granular flows of 3Dspheres and 2D disks. The average diameter and the den-sity of particles are denoted as d and ρ s which gives thecharacteristic mass m = ρ s πd / m = ρ s πd / . d to1 . d . For the contact forces, we use the standard spring-dashpot model with the Coulomb friction as in previousstudies [2, 7, 13, 14, 19, 30, 31]. In order to simulatehard particles, we choose the normal elastic constant highenough to keep the average overlapping distance smaller a r X i v : . [ c ond - m a t . s o f t ] M a y FIG. 1: Planar shear geometries tested: (a) simpleshear, (b) shear with gravity, (c) chute flows ( θ = 60 ◦ and 90 ◦ ), and (d) concave flows. The dashed lines areschematic velocity profiles.than 10 − d . The tangential elastic constant is set to be2 / . µ ( I ) rheologywhile shear flows with gravity (Fig. 1b), flows in a verti-cal chute (Fig. 1c with θ = 90 ◦ ), flows in a tilted chute(Fig. 1c with θ = 60 ◦ ), and “concave” flows (Fig. 1d)exhibit nonlocality. Concave flows are so-named afterthe shape of the shearing profile, which arises from anoutward external force (cid:126)F z ∝ ( m/d )( z − z )ˆ z for z themidpoint of the system. The gravity G is constant foreach case. The simulated domain is cuboid ( L x = 20 d and L y = 16 d ; L α is the system length in the α -direction)for 3D systems and rectangular ( L x = 160 d ) for 2D. Thehorizontal boundaries are periodic. We employ a widelyused feedback scheme to assert top-wall pressure P wall [2, 7, 14, 19, 31]. The horizontal wall velocity V wall isconstant. We use different P wall and V wall combinationsto generate varied flow profiles. In total, we ran 105different simulations, spanning two surface friction co-efficients ( µ p = 0 . .
1) and two grain shapes (3Dspheres and 2D disks). The total number of particles ineach simulation varies from around 6 . × to 2 . × .See Supplemental Material [32] for more details.When steady state is reached, the averaged continuumfields are calculated by coarse-graining. Following pre-vious studies [7, 14, 31], we calculate the instantaneousvelocity field by (cid:126)v ( z k , t ) = (cid:80) i A ki (cid:126)v i ( t ) / (cid:80) i A ki where v i is the velocity of the i th particle and A ki is the cross-sectional area (length in 2D) between the i th particleand the plane of z = z k . The interval of z k is kept lessthan 0 . d . We define the instantaneous granular temper-ature tensor as T ( t ) = (cid:80) i A ki ( δ(cid:126)v i ( t ) ⊗ δ(cid:126)v i ( t )) / (cid:80) i A ki where δ (cid:126)v i ( t ) = (cid:126)v i ( t ) − (cid:126)v ( z k , t ). When we calculate thevelocity fluctuations, we use the instantaneous velocity field as in [19]. The instantaneous stress is given by σ ( r k , t ) = σ K ( r k , t ) + (cid:80) i A ki σ i ( t ) /A where σ i is theparticle-wise stress from contacts, A is the area of thehorizontal plane ( L x in 2D), and σ K = − ρ s φ T is thekinetic stress [33]. The granular temperature is cho-sen as T xx = δv x because the diagonal components areslightly different each other possibly due to rigid-wall ef-fects. Similarly, we choose P as − σ zz , τ as | σ xz | , and ˙ γ as | ∂ z v x | . All the fields are then averaged over time. Forwell-averaged steady flow data within a limited numberof snapshots excluding wall effects, we cut off the datawhere total local shear is less than 1, φ < .
4, or thedistance from the walls is less than 3 d .The relations between the coarse-grained fields areshown in Fig. 2. As many previous studies have ob-served, µ and I are not one-to-one in inhomogeneousflows (Fig. 2a and Fig. 2e). Also, Θ is not determinedonly by I (Fig. 2b, Fig. 2f). However, there is a certaintrend. For a given I , smaller µ corresponds to larger Θ asif heating softens the material. In the spirit of the power-law scaling in continuous phase transitions, we have triedmultiplying either µ or I by a power of Θ, which are thesimplest cases, changing the exponent p to achieve thebest data collapse. Surprisingly, all the 3D sphere datawith µ p = 0 . µ is multiplied by Θ / : µ Θ / = f ( I ) (Fig. 2c). Rescal-ing I does not give a better data collapse than rescaling µ . The same exponent p = 1 / µ p = 0 . µ with a power of Θ also producesa well-collapsed master curve for disks, but the best ex-ponent p is about 1 / µ p = 0 . . µ Θ p = f ( I ) (1)where p depends on the spatial dimensions and f ( I ) de-pends as well on particle information. See the Supple-mental Material for the fitting functions in Fig. 2.Next, we run simulations on inclined chute flows wherethe velocity depends on two spatial coordinates, y and z to check the predictive values of our µ ( I, Θ) rheology ina complex geometry. For easy calculations, we imposethe no-slip boundary condition by setting two identicalgranular systems flowing in opposite directions periodi-cally neighboring each other as in [34] (see Supplemen-tal Video). We perform DEM simulations in a cuboiddomain( L x = 120 d and L y = 40 d ) using the same ma-terial used for the planar shear flows of 3D spheres with µ p = 0 . . × particles aresimulated in total. The continuum fields are averagedalong 300 lines (50 y coordinates × z coordinates)parallel to the x axis. The overlap lengths between thelines and the particles are used for the weighting in thecoarse-graining. We use a basis aligned with the localshearing planes, per [35], so that µ , I , and Θ are definedthe same way as before. The same cut off standards areused. Figure 3c shows that µ Θ / = f ( I ) still holds in the (a) (b) (c) (d)(e) (f) (g) (h) FIG. 2: Relations between µ , I , and Θ in various planar shear configurations. Non-collapse of µ vs I in 3D (a) and2D (e), and similar non-collapse of Θ vs I (b, f). The dashed lines in (a) and (b) are µ loc ( I ) and Θ loc ( I ) in Eq. (3)respectively. Multiplying µ by Θ p makes the scattered points collapse into a master curve f ( I ) (dashed trend linesfor µ p = 0 . µ p = 0 . p = 1 / p = 1 / θ = tan − (0 .
5) andtan − (0 . µ ( I, Θ) model and the µ ( I ) model using the steady-stateCauchy momentum equation ∂ j σ ij + ρ s φG i = 0. Forthe weight density term we fix φ = 0 .
60, inferred fromthe mean height and DEM floor pressure. We assumethe stress deviator and the strain-rate tensor are co-directional. The boundary conditions are traction-free onthe free surface and (cid:126)v = 0 on the other three boundaries.Rather than assume a fluctuation energy balance relationto model the temperature field, we use Θ( y, z ) extractedfrom the DEM data (see Fig. 3b). See Supplemental Ma-terial for simulation details. The steady-state velocityprofile predicted by the µ ( I, Θ) relation is almost identi-cal to the DEM data in Fig. 3d and Fig. 3f. However, the µ ( I ) rheology, which assumes vanishing shear rate where µ < µ s , disagrees with the DEM data as shown in Fig. 3dand Fig. 3e.The connection between our rheology and the well-known µ ( I ) rheology becomes clearer when Eq. (1) is rewritten as µ ( I, Θ) = (cid:18) Θ loc ( I )Θ (cid:19) p µ loc ( I ) (2)where µ loc ( I ) and Θ loc ( I ) are µ and Θ, respectively, lo-cally determined by I in simple shear flows. The µ ( I )rheology is retrieved when Θ = Θ loc ( I ). Equation (2)indicates the model can be calibrated entirely from sim-ple shear tests, if p is indeed universal and known fora family of materials. Additionally, Eq. (2) reflects thekey physical idea that Θ produces fluidization; higher Θscales down the flow strength at fixed I . The Θ field pro-duces fluidization while presumably spreading diffusivelydue to an underlying fluctuation energy balance law gov-erning the temperature [20, 21, 36]; this bears a strongsimilarity with the dynamics/role of the NGF fluidityfield, furthering the possibility of a connection betweenNGF’s fluidity diffusion equation and fluctuation energybalance [37].Another consequential relation identified in our DEMsimulations is a one-to-one relation between φ and µ (Fig. 4a) at steady state. Contrary to the standard ki-netic theory where φ is determined by Θ, it is not Θ but µ that collapses our φ data the best. In 3D developed flow, (a) (b)(c) (d) . . . . . . . . (e) . . . . . . . . (f) FIG. 3: (a) Inclined chute geometry with no-slip sidesand a rough floor; DEM velocity pictured. Red particlesfixed. (b) The distribution of log Θ. (c) Comparing the µ ( I, Θ) trend line previously obtained in planar sheartests to this geometry. (d) Comparing velocity fromcontinuum model solutions to DEM data, viewed downthe chute’s center-plane (fixed y ). (e) Comparing DEMvelocity contours to that obtained from solving the µ ( I )relation, and (f) the µ ( I, Θ) relation. The unit ofvelocity is ( | (cid:126)G | H ) / where H = 20 d . All figures are fortan θ = 0 . θ = 0 . φ ( µ ) ≈ φ − aµ for φ (cid:38) . φ = 0 .
69 and a = 0 .
27. The same formula appliesin 2D with φ = 0 .
87 and a = 0 .
22 for φ (cid:38) .
78. Theeffect of particle surface friction on the φ ( µ ) relation isnot large, confirming previous observations [2].This φ ( µ ) relation explains how our µ ( I, Θ) relationis connected to Zhang’s fluidity expression ˜ g ≡ gd/δv = F ( φ ), which has been observed to hold in previous studies[19, 38, 39]. First, we divide the range of I into three (a) (b) FIG. 4: DEM data in various planar shear flows of 3Dspheres with µ p = 0 .
4. (a) φ decreases linearly with µ for φ (cid:38) . φ ≈ . − . µ (dashed line). (b) Ourfindings predict ˜ g = gd/δv ∼ (0 . − φ ) (dashedparabola) for 10 − . (cid:46) I (cid:46) − and ˜ g ∼ constant for10 − (cid:46) I (cid:46) µ Θ / = f ( I ) ∼ I / for 10 − (cid:46) I (cid:46) − . I / for 10 − . (cid:46) I (cid:46) − I / for 10 − (cid:46) I (cid:46) . (3)In the 1 / I fol-lowing Θ ∼ I / (Fig. 2b). Combining this fact withEq. (3) and the fact that ˜ g can be rewritten as I/µ √ g ∼ constant. This plateau regime is in linewith kinetic theory where ˜ g = F ( φ ) / √ F ( φ ) becomesalmost constant for φ (cid:38) .
49 [20–23, 38, 40]. In the1 / g = I/µ √ g ∼ µ , which canbe further re-expressed under the linear φ ( µ ) collapse as˜ g ∼ ( φ − φ ) . Therefore, in the 1/3 regime, ˜ g decreasesquadratically in φ . Merging this regime’s behavior withthe plateau of the 1/2 regime, as shown in Fig. 4, deliv-ers the basic large-˜ g behavior of the ˜ g − φ relationshipapparent in our data and observed in [19]. However, inthe 1 / g (cid:46) µ ( I, Θ) relation, onthe other hand, remains well-collapsed and explains thespread in Zhang’s representation as due to ˜ g gaining ad-ditional Θ dependence; in the 1/4 regime, Eq. (3) and˜ g = I/µ √
3Θ imply ˜ g ∼ ( φ − φ ) Θ / .Gaume and coworkers [24] have also treated µ , I , andΘ as independent variables to attempt a relation betweenthem. They have suggested Θ ∝ I h ( µ ) where h ( µ ) lin-early changes with µ . Although this formula approxi-mately fits their DEM data in annular shear flow, ourdata does not match this trend and it appears their for-mula cannot be carried accurately to large I ; µ is notdetermined at I = 1, and µ increases as Θ increases for I >
1. By comparison, advantages of our model include aform motivated by power-law scaling in phase transitions,covering up to higher I and producing a strong data col-lapse over a wide array of geometries. Our model alsoreveals a potentially universal scaling exponent p , which,once identified, allows model fitting solely from simpleshear data using Eq. (2). Additionally, our model offersa connection to and expansion from existing approaches,namely kinetic theory and the NGF model, while clearlyencapsulating, through Eq. (2), the physical role of heat-softening.Using many DEM simulations, we have found a gen-eral constitutive equation for simple granular materials,which relates three dimensionless variables: µ , I , andΘ. The granular rheology can be expressed as a power-law scaling form µ Θ p = f ( I ) where the exponent p isabout 1 / / f ( I ) has certain general behaviors but details depend on the ma-terial properties. Our calibrated relation can be used togenerate the velocity field in inclined chutes where flowdepends on two spatial coordinates. We also observe aone-to-one relation between φ and µ , which allows usto reconcile our model with φ -dependent constitutive re-lations proposed by both the empirical and theoreticalapproaches. Kinetic theory, NGF modeling, and our cur-rent work all point strongly to the idea that the diffusingfield responsible for granular nonlocality is directly re-lated to the temperature. A clear next step is to explorethe inclusion of a fluctuational energy balance law ac-curate into the dense regime; this would provide Θ andcomplete the rheological model. [1] G. D. R. MiDi, Eur Phys J E Soft Matter , 341 (2004).[2] F. da Cruz, S. Emam, M. Prochnow, J. N. Roux, andF. Chevoir, Phys Rev E Stat Nonlin Soft Matter Phys , 021309 (2005).[3] P. Jop, Y. Forterre, and O. Pouliquen, Nature , 727(2006).[4] O. Pouliquen, Physics of Fluids , 542 (1999).[5] T. S. Komatsu, S. Inagaki, N. Nakagawa, and S. Nasuno,Phys Rev Lett , 1757 (2001).[6] P. Jop, Y. Forterre, and O. Pouliquen, Physics of Fluids , 088102 (2007).[7] G. Koval, J. N. Roux, A. Corfdir, and F. Chevoir, PhysRev E Stat Nonlin Soft Matter Phys , 021306 (2009).[8] K. Nichol, A. Zanin, R. Bastien, E. Wandersman, andM. van Hecke, Phys Rev Lett , 078302 (2010).[9] K. A. Reddy, Y. Forterre, and O. Pouliquen, Phys RevLett , 108301 (2011).[10] E. Wandersman and M. van Hecke, EPL (EurophysicsLetters) , 24002 (2014).[11] E. Martinez, A. Gonzalez-Lezcano, A. J. Batista-Leyva,and E. Altshuler, Phys Rev E , 062906 (2016).[12] Z. Tang, T. A. Brzinski, M. Shearer, and K. E. Daniels,Soft Matter , 3040 (2018).[13] K. Kamrin and G. Koval, Phys Rev Lett , 178301(2012).[14] K. Kamrin and G. Koval, Computational Particle Me-chanics , 169 (2014).[15] K. Kamrin and D. L. Henann, Soft Matter , 179 (2015).[16] D. L. Henann and K. Kamrin, Proc Natl Acad Sci U SA , 6730 (2013).[17] J. Goyon, A. Colin, G. Ovarlez, A. Ajdari, and L. Boc-quet, Nature , 84 (2008).[18] L. Bocquet, A. Colin, and A. Ajdari, Phys Rev Lett ,036001 (2009).[19] Q. Zhang and K. Kamrin, Phys Rev Lett , 058001(2017).[20] J. T. Jenkins and S. B. Savage, Journal of Fluid Mechan-ics , 187 (1983).[21] C. K. K. Lun, S. B. Savage, D. J. Jeffrey, and N. Chep-urniy, Journal of Fluid Mechanics , 223 (1984).[22] V. Garz and J. W. Dufty, Physical Review E , 5895(1999).[23] J. T. Jenkins and D. Berzi, Granular Matter , 151(2010). [24] J. Gaume, G. Chambon, and M. Naaim, Physical ReviewE , 051304 (2011).[25] M. Kardar, Statistical Physics of Fields (Cambridge Uni-versity Press, Cambridge, 2007).[26] F. Radjai and S. Roux, Physical Review Letters ,064302 (2002).[27] L. Staron, J.-P. Vilotte, and F. Radjai, Physical ReviewLetters , 204302 (2002).[28] O. Pouliquen, Physical Review Letters , 248001 (2004).[29] L. E. Silbert, A. J. Liu, and S. R. Nagel, Physical ReviewLetters , 098301 (2005).[30] P. A. Cundall and O. D. L. Strack, Gotechnique , 47(1979).[31] D. Liu and D. L. Henann, Soft Matter , 5294 (2018).[32] See Supplemental Material for more discussion on DEMand continuum simulation methods, additional DEMdata and continuum solutions, and f ( I ) fit functions.[33] T. Weinhart, R. Hartkamp, A. R. Thornton, and S. Lud-ing, Physics of Fluids , 070605 (2013).[34] P. Chaudhuri, V. Mansard, A. Colin, and L. Bocquet,Physical Review Letters , 036001 (2012).[35] M. Depken, W. van Saarloos, and M. van Hecke, PhysicalReview E , 031302 (2006).[36] L. Bocquet, W. Losert, D. Schalk, T. C. Lubensky, andJ. P. Gollub, Physical Review E , 011307 (2001).[37] K. Kamrin, Frontiers in Physics (2019),10.3389/fphy.2019.00116.[38] D. Berzi and J. T. Jenkins, Physical Review Fluids ,094303 (2018).[39] F. Qi, S. K. de Richter, M. Jenny, and B. Peters, PowderTechnology , 722 (2020).[40] D. Berzi, J. T. Jenkins, and P. Richard, Journal of FluidMechanics , A27 (2020). Supplemental Material for“Power-law scaling in granular rheology across flow geometries”
S1. SIMULATION CONDITIONS
We use LAMMPS, which implements the discrete element method (DEM), to simulate granular flows of 3D spheresand 2D disks. For the contact forces, we use the standard spring-dashpot model where the normal force is F n = k n δ n − γ n v n and the tangential force is F t = k t δ t where δ n and δ t are the normal and tangential components of thecontact displacement respectively and v n is the normal component of the relative velocity. The tangential elasticconstant k t is set to be 2 / k n . The restitution coefficient (cid:15) is chosen to be0 .
24. The damping coefficient is then given by γ n = (cid:112) mk n / (1 + ( π/ ln (cid:15) ) ) [2, 31]. The simulation time step is setto be 6% of the binary collision time τ c = π (cid:113) m k n (1 + (ln (cid:15)/π ) ). The external body force in the concave flows is (cid:126)F z = ( mG/d )( z − z )ˆ z where G is a constant and z is the midpoint of the system.Table S1 to S4 summarize the simulation conditions. N is the total number of particles except wall particles. Theunit of pressure P is 3 . × − k n /d in 3D and 3 . × − k n in 2D. The unit of acceleration G is (1 / P /ρ s in 3Dand (1 / P /ρ s d in 2D. The unit of velocity V is 8 . (cid:112) P /ρ s in 3D and 4 . (cid:112) P /ρ s in 2D. We output data every∆ n steps to obtain total N out snapshots.TABLE S1: Simulation conditions for the planar shear tests (3D spheres with µ p = 0 . Geometry
N P wall /P G/G V wall /V ∆ n N out Simple shear 18327 4 0.003125 160000 1800Simple shear 18327 4 0.0125 80000 1800Simple shear 18327 4 0.05 40000 1800Simple shear 18327 1 0.1 40000 1800Simple shear 18327 1 0.2 40000 1800Simple shear 18327 1 0.4 40000 1800Simple shear 18327 1 0.8 40000 1800Simple shear 18327 1 1.6 40000 1800Simple shear 18327 1 3.2 20000 3600Simple shear 18327 1 6.4 20000 3600Shear with gravity 18327 1 16 3.2 20000 7200Shear with gravity 18327 1 2 12.8 20000 7200Shear with gravity 18327 1 32 1.6 40000 3600Shear with gravity 18327 1 4 6.4 40000 7200Shear with gravity 18327 1 8 1.6 40000 3600Shear with gravity 18327 4 8 0.1 40000 3600Chute flows ( θ = 90 ◦ ) 18327 8 12 40000 3600Chute flows ( θ = 90 ◦ ) 18327 8 16 40000 3600Chute flows ( θ = 90 ◦ ) 18327 8 20 40000 3600Chute flows ( θ = 60 ◦ ) 18327 8 16 40000 3600Chute flows ( θ = 60 ◦ ) 18327 8 20 40000 3600Chute flows ( θ = 60 ◦ ) 18327 8 24 40000 3600Concave flows 18327 16 3.5 0.00625 80000 3600Concave flows 18327 16 3 0.025 40000 3600Concave flows 18327 16 3 0.4 40000 3600Concave flows 18327 16 3 1.6 20000 7200 µ p = 0 . Geometry
N P wall /P G/G V wall /V ∆ n N out Simple shear 18327 4 0.0015625 160000 1800Simple shear 6923 4 0.0015625 160000 3600Simple shear 6923 4 0.003125 160000 1800Simple shear 6923 4 0.0125 80000 1800Simple shear 6923 4 0.05 40000 1800Simple shear 6923 1 0.1 40000 1800Simple shear 6923 1 0.2 40000 1800Simple shear 6923 1 0.4 40000 1800Simple shear 6923 1 0.8 40000 1800Simple shear 6923 1 1.6 40000 1800Simple shear 6923 1 3.2 40000 3600Simple shear 6923 1 6.4 40000 3600Shear with gravity 6923 1 1 6.4 20000 14400Shear with gravity 18327 1 16 0.1 40000 3600Shear with gravity 18327 1 16 1.6 40000 3600Shear with gravity 18327 1 4 1.6 40000 3600Chute flows ( θ = 90 ◦ ) 18327 8 10 40000 3600Chute flows ( θ = 90 ◦ ) 18327 8 12 40000 3600Chute flows ( θ = 90 ◦ ) 18327 8 14 40000 3600Chute flows ( θ = 60 ◦ ) 18327 8 10 40000 3600Chute flows ( θ = 60 ◦ ) 18327 8 12 40000 3600Chute flows ( θ = 60 ◦ ) 18327 8 14 40000 3600Concave flows 18327 16 3.5 0.025 40000 3600Concave flows 18327 16 3 0.4 40000 3600Concave flows 18327 16 3 1.6 20000 7200 TABLE S3: Simulation conditions for the planar shear tests (2D disks with µ p = 0 . Geometry
N P wall /P G/G V wall /V ∆ n N out Simple shear 6739 1 0.0015625 160000 3600Simple shear 6739 1 0.00625 80000 3600Simple shear 6739 1 0.025 20000 3600Simple shear 6739 1 0.1 20000 3600Simple shear 6739 1 0.2 20000 3600Simple shear 6739 1 0.4 20000 3600Simple shear 6739 1 0.8 20000 3600Simple shear 6739 1 1.6 20000 3600Simple shear 6739 1 3.2 10000 7200Simple shear 6739 1 6.4 10000 7200Simple shear 6739 1 12.8 40000 7200Simple shear 6739 1 25.6 40000 7200Shear with gravity 6739 1 1 12.8 20000 14400Shear with gravity 6739 1 1 6.4 20000 14400Shear with gravity 19810 4 16 4.8 40000 3600Shear with gravity 19810 4 2 0.075 40000 3600Shear with gravity 19810 4 2 4.8 40000 3600Shear with gravity 19810 4 4 0.3 40000 3600Chute flows ( θ = 90 ◦ ) 19810 4 2 40000 3600Chute flows ( θ = 90 ◦ ) 19810 4 3 40000 3600Chute flows ( θ = 90 ◦ ) 19810 4 4 40000 3600Chute flows ( θ = 60 ◦ ) 19810 4 3 40000 3600Chute flows ( θ = 60 ◦ ) 19810 4 4 40000 3600Chute flows ( θ = 60 ◦ ) 19810 4 6 40000 3600Concave flows 19810 16 2/3 0.15 40000 3600Concave flows 19810 16 2/3 1.2 40000 3600Concave flows 19810 16 2/3 4.8 40000 3600Concave flows 19810 4 1/6 2.4 40000 3600 µ p = 0 . Geometry
N P wall /P G/G V wall /V ∆ n N out Simple shear 6739 1 0.0015625 160000 3600Simple shear 6739 1 0.00625 80000 3600Simple shear 6739 1 0.025 40000 3600Simple shear 6739 1 0.1 20000 3600Simple shear 6739 1 0.2 20000 3600Simple shear 6739 1 0.4 20000 3600Simple shear 6739 1 0.8 20000 3600Simple shear 6739 1 1.6 20000 3600Simple shear 6739 1 3.2 20000 3600Simple shear 6739 1 6.4 20000 3600Simple shear 6739 1 12.8 20000 3600Shear with gravity 6739 1 16 0.1 40000 3600Shear with gravity 6739 1 1 12.8 20000 14400Shear with gravity 6739 1 1 6.4 20000 7200Shear with gravity 6739 1 8 0.4 40000 3600Shear with gravity 6739 4 4 0.025 40000 3600Shear with gravity 6739 4 4 0.1 40000 3600Chute flows ( θ = 90 ◦ ) 19810 4 2 40000 3600Chute flows ( θ = 90 ◦ ) 19810 4 3 40000 3600Chute flows ( θ = 90 ◦ ) 19810 4 4 40000 3600Chute flows ( θ = 60 ◦ ) 19810 4 2 40000 3600Chute flows ( θ = 60 ◦ ) 19810 4 3 40000 3600Chute flows ( θ = 60 ◦ ) 19810 4 4 40000 3600Concave flows 19810 16 2/3 0.075 40000 3600Concave flows 19810 16 2/3 0.3 40000 3600Concave flows 19810 16 2/3 4.8 40000 3600 TABLE S5: Simulation conditions for the inclined chute flows (3D spheres with µ p = 0 . Geometry
N G/G tan θ ∆ n N out Inclined chute flows 115619 64 0.5 40000 3600Inclined chute flows 115619 64 0.6 40000 3600
S2. FITTING FUNCTIONS IN FIG. 2
We approximate the master curves as f ( I ) ≈ . I . + 0 . I . for 3D spheres with µ p = 0 . . I . + 0 . I . for 3D spheres with µ p = 0 . . I . + 0 . I . for 2D disks with µ p = 0 . . I . + 0 . I . for 2D disks with µ p = 0 . µ loc ( I ) ≈ .
35 + 0 . I . − . I . , and the one in Fig. 2b isdetermined by Θ loc ( I ) = ( f ( I ) /µ loc ( I )) . S3. SUPPLEMENTAL FIGURES
We provide additional figures from the DEM simulations. Fig. S1 to S4 show the relations between µ , I , Θ, and φ obtained from the planar shear flows. Fig. S5 and Fig. S6 show supplemental DEM data and solutions to the Cauchymomentum equation in the inclined chute geometry.4 (a) (b) (c) (d) FIG. S1: DEM data from planar shear tests. Non-collapse of µ vs I (a) and Θ vs I (b) for 3D spheres with µ p = 0 . µ vs I (c) and Θ vs I (d) for 2D disks with µ p = 0 . (a) (b) (c) FIG. S2: DEM data from planar shear tests. One-to-one relationship between φ and µ for 3D spheres with µ p = 0 . µ p = 0 . µ p = 0 . φ = 0 . − . µ . Dashed lines in (b) and (c) are the same: φ = 0 . − . µ . (a) (b) (c) (d) FIG. S3: DEM data from planar shear tests. Non-collapse of φ vs Θ for 3D spheres with µ p = 0 . µ p = 0 . µ p = 0 . µ p = 0 . (a) (b) (c) (d) FIG. S4: DEM data from planar shear tests. Non-collapse of φ vs I for 3D spheres with µ p = 0 . µ p = 0 . µ p = 0 . µ p = 0 . (a) (b) (c) FIG. S5: DEM data from inclined chute flows with no-slip sides. (a) The distribution of log Θ for tan θ = 0 .
6. (b)Non-collapse of µ vs I (b) and Θ vs I (c). Trend lines are from the planar shear tests. (a) . . . . . . (b) . . . . . . (c) FIG. S6: Simultaneous solutions of the Cauchy momentum equation with the µ ( I, Θ) relation in inclined chute flowswith no-slip sides for tan θ = 0 . θ = 0 . y ). (b) Comparing DEM velocity field to that obtained fromsolving the µ ( I ) relation and (c) the µ ( I, Θ) relation.6
S4. CONTINUUM SIMULATION METHOD
We use the finite difference method to solve the Cauchy momentum equation in the inclined chute flows. Thevelocity field is calculated on a 41 ×
50 grid representing the yz plane. The stress field is staggered, located on cellcenters (a 40 ×
49 grid of locations). Θ is interpolated to the grid of stress. For regularization, µ loc ( I ) is modifiedto gradually vanish from I = 10 − to I = 10 − which prevents numerical errors by giving a finite but insignificantshear rate for µ < µ s (Fig. S7). We assume that the stress deviator σ (cid:48) and the strain-rate tensor D are co-directional: D / | D | = σ (cid:48) / | σ (cid:48) | . We apply σ xz , σ yz = 0 at the surface by assuming the surface is flat and imposing an imaginarystress of opposite sign mirrored across the surface. We know analytically that with co-directional flow rules, steadyflows always develop lithostatic pressure, which we exploit by pre-setting σ zz = ρ s φG cos θ ( H − z ). We update thevelocity field putting either µ = µ loc ( I ) or µ = Θ − / f ( I ) in the momentum equation until the Frobenius norm of thevelocity change becomes small enough. We have checked that the final results are independent of the initial velocity. -8 -6 -4 -2 FIG. S7: Modified µ loc ( I ) for regularization. S5. VIDEO
Download “Inclined Chute Flows tan05.avi” to watch the motion of particles in the inclined chute flow with tan θ =0 . (cid:126)G = G sin θ ˆ x − G cos θ ˆ z , while the other half receives (cid:126)G = − G sin θ ˆ x − G cos θ ˆ zz