Numerical Tests of Constitutive Laws for Dense Granular Flows

0.67 MB·English

by  Gregg Lois

Numerical Tests of Constitutive Laws for Dense Granular Flows
Gregg Lois(1), Ana¨el Lemaˆıtre(1,2), and Jean M. Carlson(1)
(1) Department of physics, University of California, Santa Barbara, California 93106, U.S.A. and
(2) L.M.D.H. - Universite Paris VI, UMR 7603,
4 place Jussieu - case 86, 75005 Paris - France
(Dated: February 2, 2008)
5 We numerically and theoretically study macroscopic properties of dense, sheared granular mate-
0 rials. In this process we first introduce an invariance in Newton’s equations, explain how it leads
0 to Bagnold’s scaling, and discuss how it relates to the dynamics of granular temperature. Next we
2 implementnumericalsimulationsofgranularmaterialsintwodifferentgeometries–simple shearand
flow down an incline–and show that measurements can be extrapolated from one geometry to the
n
other. Then we observe non-affine rearrangements of clusters of grains in response to shear strain
a
J andshowthatfundamentalobservations,whichservedasabasisfortheShearTransformationZone
(STZ)theoryofamorphoussolids[1,2],canbereproducedingranularmaterials. Finallywepresent
3
constitutiveequations for granular materials as proposed in [3], based on thedynamics of granular
2
temperature and STZ theory, and show that they match remarkably well with our numerical data

from both geometries.
]
t
f
o
s I. INTRODUCTION els. Some models posit that frustrated rotation plays
t. a preeminent role in jamming [54, 55]. Many mod-
a
els propose extensions of kinetic theory by using either
m Historically reservedto the engineering community [4,
granular temperature [56], introducing strongly density-
- 5, 6, 7], granular materials recently emerged as a new
d dependent viscosities [57, 58, 59], or incorporating a
field of study for physicists [8, 9, 10, 11, 12, 13, 14]: a
n quasi-static stress at the internal friction angle of the
state of matter that is not classifiable according to the
o material [60, 61, 62, 63, 64, 65]. Other models are based
c traditional ternary–solid,liquid or gas– and requires sci- onthe introductionofchain-forces[44, 45, 46], activated
[ entists to rethink the foundations of statistical physics
processes [67], “granular eddies” [47], coexisting liquid
and thermodynamics [15, 16]. The first theories of gran-
1 and solid microphases in a Ginzburg-Landau formula-
v ular materials were motivated primarily by the need
tion [68, 69, 70], or ’spots’ of free-volume associated to
5 to predict the creep motion of soils and their stability
cooperative diffusion [71, 72].
3 properties. Inspired by continuum theories of plastic-
5 ity, these approaches were restricted to quasi-static de- Several of the above-mentioned models conjecture
1
formation and incipient failure [4, 5, 6, 7, 17, 18, 19]. jamming mechanisms–frustrated rotations, chain forces,
0
5 In physics, initial progress came in understanding dilute granulareddies–whichare so specific to granularmateri-
0 granular materials with the development of kinetic the- als that they do not allow for connection with jamming
/ ory[20, 21,22,23,24, 25,26,27,28,29, 30,31,32]. But in other materials. Our approach rests on the viewpoint
t
a kinetic theory is based on the assumption that interac- that, since jamming is observedin manyamorphoussys-
m
tions between grains occur during instantaneous binary tems,itislikelythatithasacommonorigin. Itwasthus
- collisions,aconditionthatisrealizedbyveryfew“rapid” proposed[3]thatarescalingofthe dynamicsofperfectly
d
n flows outside the laboratory. hard granular materials would make it possible to map
o theirpropertiesontomoretypicalglassformers. Building
Traditionalapproachesthus leaveus with little under-
c on the analogy between granular materials and metallic
standingoftheelementarymechanismsofdeformationin
:
v glasses then provides a local mechanism of jamming on
granular materials. It is no wonder that the elaboration
i thebasisoftheShearTransformationZone(STZ)theory
X of physically inspired constitutive equations for granu-
of plasticity [1, 2].
r larmaterialsis currentlyfacing acollectionofcontrover-
a sialissueswhich challengeits mostfundamental aspects. The goal of the present paper is both to discuss some
Common points of contention include: (i) the relevant of the microscopic assumptions underlying the construc-
features of the grain-graininteraction–e.g. Hertzian ver- tion of constitutive equations for granular materials and
sus hard-sphere repulsion– [14, 33, 34]; (ii) the domain test the specific predictions of the STZ theory formula-
of applicability of kinetic theories, [29, 31, 35, 36, 37]; tion of such equations [3]. Our approachis based on the
(iii) the observation and micromechanical origin of Bag- followingexpectations concerningthe fourpoints ofcon-
nold scaling in dense flows [38, 39, 40, 41, 42, 43]; (iv) tention identified above: (i) major properties of dense
the plausible need to formulate “non-local” constitutive granular flows can be captured by perfectly hard grains;
equations [44, 45, 46, 47, 48] due to the existence of (ii) jamming involves variations of the frictional stress
chain forces [49, 50, 51, 52, 53]. Uncertainties about inregimeswherethecollisionalcontributions–thosepre-
these questions have fueled a wealth of theoretical mod- dicted by kinetic theory– are negligible; (iii) for flows of
2
perfectly hard grains,Bagnold’s scaling is relevant at all discuss why the results for perfectly hard grains are ap-
densitiesandarisesfromaninvarianceintheequationsof plicable to natural and experimental granular materials.
motion;and(iv)therheologyofdensegranularmaterials InSectionIIIweconstructournumericaltest,provideal-
is local in a sense to be defined further. gorithmic details on the ContactDynamics method, and
Indeed, (i) a recent work by Campbell has shown– compare the rheology of a dense granular materials in a
convincingly to us–that most natural and experimental periodic cell and down an incline plane. Finally, in Sec-
flowsoccurinaregimeofstrainrateswheregrainscanbe tion IV, we present Falk and Langer’s STZ theory, show
considered as hard-spheres [73]. Although this question how it adapts to granular materials, and conclude with
maystillbedebated,wewilltakethisforgrantedandre- fits of our data.
strictourstudytoamaterialcomposedofperfectlyhard
grains. (ii) The breakdown of kinetic theory in dense
flows has now been characterized in several numerical II. FUNDAMENTAL RESULTS FOR PERFECTLY
and experimental studies [74, 75, 76, 77, 78, 79] and in- HARD GRAINS
dicates that jamming is associated with changes in the
contributionoflong-lastingcontactsto the stresstensor, When grains are dry–so that no water bridges induce
not the collisional contribution. We will take this idea attraction–and of size larger than the micrometer scale–
for grantedin the presentwork,butwe do believe thata sothatnoelectrostaticinteractionintervenes–theirinter-
quantitativeanalysisofthisstatementisneededandwill action is purely repulsive. The interaction results from
devoteafuturepapertothisquestion. (iii)Wewillshow, the elastic deformation of grains at contact and the dis-
following [3] that in the case of perfectly hard grains, sipation of energy via friction and collisions. The com-
Bagnold’sscalingholds for anydensity. Assumption(iv) plexity of this interaction motivates our first question:
is strongly inspired by the works by Aranson, Tsimring whichpropertiesofthegrain-graininteractioncontribute
andcoworkers,whohaveshownthataGinzburg-Landau to any particular macroscopic observation? In some in-
formulationof a fluid-solidmixture [68, 69, 70] couldac- stances, details of the grain-grain interaction seem criti-
count for important properties of granular flows. This cal: for example, the Hertzian repulsion [33] is essential
suggests that a local formulation of granular rheology, tounderstandtheacousticpropertiesofgranularmateri-
coupled to hydrodynamic equations, can be sufficient to als [80, 81]. Numerical implementations of granular ma-
account for a large part of the phenomena observed in terials have thus relied on more or less elaborate models
dense granular flows. However, no numerical study has of the grain-graininteraction [82, 83].
been specifically devoted to this question. Here we are concerned with dense flows of granular
Therefore, we will first study whether dense granular materialsand,inparticular,flowsdowninclinesasfound
flows can be described as a local phenomenon, governed in the experiments by Pouliquen [84] and numerics by
by local hydrodynamic equations, or whether granular Silbert and coworkers [41, 42]. For these dense granu-
flows must be treated non-locally, relying on the emer- lar flows a recent study by Campbell helps us assess the
gence of long-range correlations through mechanisms importance of the elastic (soft) part of the repulsive po-
such as force chains. For this purpose we have imple- tential [73] versus the limit in which grains appear as
mented Contact Dynamics simulations of sheared gran- perfectly hard spheres. Campbell presented a detailed
ular materials in two different configurations: (i) simple analysis of the different flow regimes obtained in a three
shear in a periodic cell and (ii) thick flows of granular dimensional simple shear simulation of dense granular
material down an incline plane. We will show that the flows when varying the stiffness k of the repulsion, the
measurementsobtainedineitherconfigurationcanbeex- shear rate γ˙, and the mass density φ. He found that
trapolated to the other. the dimensionless parameter Υ k , where D is the
≡ φD3γ˙2
Next we follow Falk and Langer [1, 2] and test for the grain size, dictates the character of the flow. This quan-
consistency of the STZ picture of material deformation tity is directly related to a Mach number which involves
advocated by these authors. Our observation provides the ratio of the shear velocity Dγ˙ over the sound speed
further grounding for the analogy between granular ma- c : M = Dγ˙/c = 1/√Υ. The hard-sphere limit cor-
s s
terials and common glass-formers and direct support for responds to the situation where sound waves travel very
the relevance of STZ theory to granularmaterials. A re- fast compared to the motion induced by the shear flow.
viewoftheSTZconstitutiveequationsforgranularflows This is the limit of very small Mach number, or small
willfollowandfinallyafitofthe theorywithournumer- shear rates. Specifically, in the numerics by Campbell,
ical data. this limit is reached for Mach numbers below 10 2.
−
∼
The organization of our paper is as follows. In Sec- Since the sound speed in granular materials is of order
tion II we present equations of motion in the mathemat- 100m/s, if we assume a grain size of order 1mm, the
ical limit of perfectly hard grains. We discuss the in- Mach number is expressible as M = 10 5γ˙. Therefore,
−
variance properties of the equations of motion and how in orderto be inthe limit where grainsbehave as if they
thesepropertiesarerelatedtoBagnold’sscaling. Wealso are perfectly stiff, it suffices to restrain oneself to shear
3
ratesbelow1000s 1. Mostexperimentalandnaturalsit- collision are each proportional to the shear rate, there-
−
uations occur atshear rates far below this limiting value fore the stress is proportional to the square of the shear
andwecanconclude thatmostflowsofgranulatesarein rate.
the limit where the soft part of the repulsion is entirely Recent observations have complicated Bagnold’s orig-
maskedbythe stericexclusion. Tostudythese flowsitis inal assessment. On the one hand, the observations by
sufficient to consider properties of perfectly hard grains. Bagnoldhavebeencriticized: theymayhavearisenfrom
In the following sections we explore how the mathemati- a secondary instability of the granular flow in his shear
cal limit of perfectly hard spheres gives insight into fun- cell [40]. On the other hand, Bagnold’s scaling has been
damentalprocesseswhichrelatetomostallexperimental directly observed by measuring shear stress and strain
shear flows. rate profiles in numerical simulations of granular flows
downinclines [41, 42],andis foundto be consistentwith
experimentalobservationsofthe averageflow ratein the
Equations of motion and Hard Sphere conditions same geometry [84].
Although kinetic theory predicts an exponent which
The motion of N spherical grains in a d-dimensional agreeswith the powerlaw dependence of Bagnold’sscal-
granular material is determined by Newton’s equations ing, the derivation relies on the assumption that all in-
for the positions q , angular orientations θ , momenta p teractionsbetweengrainsoccurthroughbinarycollisions.
i i i
and angular velocities ω : This assumption is applicable to dilute flows but is not
i
upheld for dense flows, such as those originally studied
dq p dp
i = i , i = F +Fext, (1) byBagnold. Anewargumentmustbeformulatedforthe
ij
dt mi dt dense flow regime.
j
X
We contend that in both the dense and “rapid” flow
regimes, Bagnold’s scaling arises from a very fundamen-
dθi dωi 1 tal invariance of Newton’s equations in the hard sphere
=ω , = R nˆ F , (2)
i i ij ij
dt dt I × limit. Thisinvariancedoesnotrequireanyoftheassump-
i
j
X tions of kinetic theory to hold. Namely, for a granular
where Fext represents an external force such as gravity, materialfreefromexternalforces,thetimeevolutionwill
F represents a contact force on grain i by grain j, nˆ obey equations (1, 2) with Fext = 0. If we now rescale
ij ij
is the unit normal vector pointing from grain i to grain the contact forces by a scalar value Fij Fij/A and si-
→
j, R is the radius of grain i, and I is the moment of multaneously rescale the time t t√A, then Newton’s
i i →
inertia. equations are transformed to read:
These equations must be complemented with a pre-
dq pnew dpnew
scription for the contact forces. In the hard sphere limit i = i , i = F , (3)
ij
dt m dt
these contact forces are determined self-consistently by i
j
X
the conditions that (i) there is no penetration between
grain–a force is instantaneously created upon contact to
dθ dωnew 1
impede penetration and remains non-zero until the con- i =ωnew, i = R nˆ F , (4)
tact is broken–and (ii) by the friction law which couples dt i dt Ii i ij × ij
j
X
to rotational degrees of freedom.
Important properties of granular materials arise di- where pnew =p /√A and ωnew =ω /√A. This form for
i i i i
rectly from an invariance of the equations of mo- Newton’s equations is identical to (1) and (2) with new
tion (1, 2). We now spend some time studying these values for the momenta and angular velocities.
properties and assessing their consequences for macro- Undertherescalingofcontactforcesandtime,the po-
scopic observations, in particular Bagnold’s scaling. sitionsandangularorientationsremainunchanged,while
the velocities are changed in accordance with the time
rescaling. If we were to watch a movie of one granular
Bagnold’s scaling flow where the grains have initial velocities p , ω and
i i
watch another movie at half the speed where the ini-
Thesuccessofkinetictheorycameinalargepartfrom tial velocities are doubled p 2p , ω 2ω , the two
i i i i
→ →
itsabilitytoaccountforthescalingbetweenstressσ and movies would look exactly the same in the hard sphere
strainrateγ˙ (σ γ˙2)firstobservedbyBagnoldindense limit. The difference in the dynamics is that the contact
∼
granular materials [38]. In Bagnold’s experiment, the forces measured in the second movie would be 4 times
strainratewassetandthestressmeasured. Bagnoldjus- larger than those in the first.
tifiedthisbehaviorbyassumingthatthestressmeasured This invariance is a property of perfectly hard grains
in his experiments resulted solely from binary collisions: which must hold in the inertial regime. This includes
the frequency of collisions and momentum change per theregimeinfinitely closetojamming,wheremulti-body
4
interactions dominate collisional terms and the basic as- andtimescalesthentheinvarianceinNewton’sequations
sumptions of kinetic theory fail. In an experiment, this guarantees that the dynamics will remain unaltered, re-
scaling breaks down only when it is no longer appropri- sulting only in new velocities that are determined from
ate to model the experimental system by perfectly hard thetimerescaling. However,undertherescaling,quanti-
grains. Relying on the arguments of Campbell intro- tiessuchasthepressureandstrainratewilltakedifferent
ducedearlier[73],wecanconcludethatmostexperimen- values. The quantities of interest are those that are in-
talgranularflowsareintheregimewhereitisappropriate variant to the transformation. These quantities must be
to model the system by hard grains. formed of ratios of force scales, ratios of time scales, or
This invariance also holds for any value of the resti- ratios of force to time scales.
tution and friction coefficients. The separation of time Denote p as the pressure, σ the shear stress, γ˙ the
scales for dissipation and shear is not a necessary con- strainrate,φ the massdensity, andD the averagediam-
dition for the invariance to be upheld. The reason is eterofgrains. Someimportantinvariant(anddimension-
that dissipation occurs at contacts between grains: dis- less)quantitiesare: p , σ , pDd, σDd, T ,
φD2γ˙2 φD2γ˙2 T T φDd+2γ˙2
sipative processes are correlatedto the motion of grains, and σ. These must be single valued functions of den-
p
hence the rate of dissipative processes scales with γ˙.
sity only (independent of the strain rate). This observa-
tion,establishedbytheinvarianceinNewton’sequations
whichholdsforgranularflowsinthedenseandcollisional
Granular Temperature
regimes, automatically predicts Bagnold’s scaling: since
σ isafunctionofdensityonly,itfollowsthatσ γ˙2.
Because the positions and angular orientations of φD3γ˙2 ∝
This leads to a simple yet remarkable property: at a
grainsareinvarianttoasimultaneouschangeintimeand
given density, changing the strainrate does not lead the
force scales, the path that a granular material takes in
systemtowardanydifferentregimeor“phase”. The ma-
configurationspaceisalsoinvariant: onlythespeedalong
terial does not go closer to jamming because the strain
this path is altered. Therefore, the path that a granular
rate is scaled down. The system is always exploring the
material takes in configuration space can be separated
same trajectories in phase space, but at a slower speed.
from the rate at which events occur along that path.
In short, there is no quasi-static limit for perfectly hard
Thisobservationleadstoanaturaldefinitionofgranu-
grains. This is apparent in the jamming phase diagram
lar temperature T: the rate at which microscopic events
of granularmaterials [85] where changing the density al-
occur on the configuration space path is defined as
lows the granular material to jam, while changing the
√T/ R (the averagegrainsize R is insertedso thatT
h i h i strain rate does not. This property is also directly re-
has units of velocity squared). These microscopic events
lated to the observation by Campbell that “there is no
arenotnecessarilycollisions,butcanalsobe understood
path between inertial (rapid) flow and quasi-static flow
asforcefluctuations orminute displacementofgrainsal-
by varying the shear rate at fixed concentration” [73].
lowing for propagationof “hard-sphere” noise.
This is no surprise once we understand the invariance of
The invariance is directly related to the structure of
Newton’s equations for perfectly hard grains: only two
Newton’s equations and the existence of inertial terms–
regimesareaccessible,jammedorinertial,anddensityis
or accordingly,to the fact that kinetic energy,no matter
the only parameterwhichcontrolsthe statisticalproper-
howsmall,iswelldefined. Thegranulartemperaturede-
ties of a flow of perfectly hard grains.
finedaboveasa“velocity”alongphasespacetrajectories
is thus naturally proportional to kinetic energy:
III. CONSTRUCTION OF A NUMERICAL TEST
m I
T = h i( v2 v 2)+ h i( ω2 ω 2), (5)
k
2 h i−h i 2 h i−h i
In constructing a numerical test our goals are to mea-
where m is the mass, v the velocity, I the moment of sure stress-strain relations when granular temperature
inertia, ω the angular velocity, and brackets denote an and stresses are appropriately scaled, show that they
average over grains. Since the granular and kinetic tem- compare well with the standard response of yield stress
peratures are defined up to a constant factor, we will liquids, and show that the local rheology measured in
equate them. This means that at any time, the kinetic simple shear flow matches the bulk rheology of a granu-
energy provides an estimate of the frequency of elemen- lar flow down an incline.
tary events in a multicontact system. In orderto address these issues, we implement numer-
icalsimulationsofgranularmaterialsintwo differentge-
ometries:
What quasi-static limit?
WeimplementsimpleshearflowinacellwithLees-
•
Consider a uniform granular material in simple shear Edwards(LE)boundaryconditions. Inthis config-
flow. If we imagine simultaneously changing the force uration,thedensityandtheshearrateisprescribed
5
FIG. 1: Snapshot of a granular material simulation in the FIG. 2: Snapshot of a granular material simulation in the
simpleshearconfiguration. Eachgrainhasanaveragevelocity incline flow configuration. Fixed grains (indicated by filled
inthex-directiongivenbyγ˙y,whereγ˙ isthestrainrate. The circles) create a stationary incline at angle θ on which the
centerof thecell is defined as x=y=0. flowing grains are accumulated and allowed to flow. Gravity
drivesthe motion and is directed vertically downward.
andthe simulationcellis,byconstruction,transla-
tionally invariant. This grants direct access to av- axis(d)inthisfigurecorrespondstothedistancebetween
eraged quantities of the granular temperature and a pair of grains, divided by the sum of their radii. This
stresstensor. Usingthisconfigurationwecanchar- normalizes the figure so that d = 1 corresponds to con-
acterize the steady state relation between stresses, tacting grains. Other than this peak, the function has
granular temperature and strain-rate and extract some small variation (which implies a correlation) be-
numerically the parameters of a constitutive law tween d=1 and d=3. However, there is no correlation
for granular materials. A screenshot of this shear- beyondd=3. Becausethereisnolargescalecorrelation,
ing geometry is shown in Figure 1. this implies that the granular material is amorphous.
Inthispaperwepresentdatafornormalandtangential
Weimplementgranularflowdownaninclinedplane coefficients of restitution given by e = e = 0. It was
• n t
made of stationary grains. The simulation cell is shown by Chevoir et al [86] that the regimes reached by
periodic in the direction (x) parallel to the plane granularmaterialsarealmostindependentofthe restitu-
andtheflowisinhomogeneousintheperpendicular tion coefficients below some threshold around e = 0.7.
n
(y)direction. Inthisconfiguration,the stressesare Different friction coefficient have been used: frictionless
prescribed by the angle of the incline. We perform grains and µ=0.4.
x-averaged, y-dependent measurements of granu-
OursimulationsofgranularmaterialsrelyontheCon-
lar temperature, velocity profiles and strain-rate.
tact Dynamics algorithm [87, 88, 89, 90, 91]. We refer
Large heights of the granular layer grant access to
totheliteraturefortechnicaldetailsandpresentbelowa
the bulk rheology of the flow. This permits us to
brief overview of the method, supplemented by a single
checktheexistenceofawell-definebulkrheologyin
addition we made in order to construct a Lees-Edwards
the large height limit, and to compare it with the
simulation cell for frictional granular materials.
measurements in simple shear. A picture of this
shearing geometry is shown in Figure 2.
Contact Dynamics
In orderto makea quantitative comparisonofthe two
simulations,weusethesamematerial: atwo-dimensional
polydisperse mixture of constant density grains with the A Contact Dynamics algorithm was constructed to
radii drawn from a flat distribution with average radius carry out numerical simulation of spheres interacting
R andwidth σ. For allofthe simulations inthis paper through the enforcement of hard sphere conditions [87].
h i
we set σ/ R = 0.5, using R = 0.7. This distribu- Whenacontactoccursthereisannon-continuousforce
h i h i
tionpreventscrystallizationandproducesanamorphous created that prevents the contacting grains from pene-
granular material, as can be seen from measurement of trating. The magnitude of this force is chosen to ensure
the pair correlationfunction in Figure 3. The horizontal that the final relative velocity u of the grains is related
6
It was recognized in early implementations of LE
boundary conditions that when deformation is applied
6
through the image cells, the information needs time to
propagatefromthecellboundariestoitscenter. Inorder
4
to ensure rapid propagationof this information and pre-
2 vent the boundaries between cells from making unphys-
on ical contributions to the motion, it is necessary to mod-
elati 0 1 1.5 2 2.5 3 3.5 4 4.5 5 ify Newton’s equations by introducing so-called SLLOD
orr
c terms. Thesetermscanbeunderstoodasasortof“shear
air
p 102 bath”,withallparticlesinthecellbeingdirectlycoupled
to the overall deformation [93]. In practice, the SLLOD
terms introduce a mechanical perturbation to the equa-
100
tions of motion that gives eachgrain an averagevelocity
consistentwithsimpleshearflow. Ifweseparatethe mo-
10−2
mentum p of each grain i into the average part m γ˙y
1 1.5 2 2.5 3 i i i
d and fluctuating part p˜i, so that pi =mγ˙yi+p˜i, then the
SLLOD equations read:
FIG. 3: Thepair correlation function for a frictionless gran-
ularmaterialathighdensityinlinear(top)andlog(bottom) dqi p˜i dp˜i
= +xˆγ˙(q yˆ), = F xˆγ˙(p˜ yˆ). (7)
i ij i
scales. This function is representative of the pair correlation dt m · dt − ·
i
j
functions for other densities and for granular materials with X
friction between grains. There is a large peak correspond-
The equation for the position q is simply the result of
ing to contacting grains (d=1) and some variation between i
writing the momentum in terms of an average and fluc-
d = 1 and d = 3. For d > 3 there is no correlation between
grains. This implies that thematerial is amorphous. tuating part. The equation for p˜i contains a new term
xˆγ˙(p˜ yˆ)whichforcestheshearflow. Sinceeverygrainin
i
·
the primitive cell is actedupon by this mechanicalforce,
to the initial relative velocity u0 via the equations the constant strain rate is imposed on all of the grains
simultaneously at the beginning of the simulation. Fur-
u = e u0, u =e u0, (6) thermore it can be proven that, in the LE geometry, the
n − n n t t t
SLLOD equations give an exact representationof simple
whereen,et areconstantrestitutioncoefficientsthatwill shear flow arbitrarily far from equilibrium [93, 94].
depend on the shape and consistency of the grains, and For a granular material with non-zero friction coeffi-
the n, t subscripts represent the normal and tangential cient µ, the equations of motions should incorporate ro-
parts of the relative velocitywith respectto the contact. tations of the grains (for µ = 0 the tangential contact
At each time step, the algorithm computes the contact force is always zero and there is no rotation). It is ex-
forces by first ensuring that these relations hold at each pected that a SLLOD termshould arisein the equations
contact. Toincludefriction,the ContactDynamicsalgo- of motion for the angular velocity since, in the linear ve-
rithmensuresthattheresultingtangentialforceFt isless locity profile indicative of simple shear flow, the top and
than or equal to µFn where µ is the friction coefficient bottomofeverygrainshouldbemovingwithslightlydif-
between grains and Fn is the normal force. If this con- ferent velocities. This will give eachgrainanaveragero-
straint does not hold, then the algorithm sets Ft = µFn tationofγ˙/2whichmustbeincorporatedinequation(2)
in order to comply with Coulomb friction. justastheaveragevelocityxˆγ˙(q yˆ)wasincorporatedin
i
·
In this way the Contact Dynamics algorithm calcu- equation (7). This leads to the following equations:
lates, at each time step, contact forces that are consis-
tent with Newton’s equations and the hard sphere con- dθi γ˙ dω˜i 2
=ω˜ + , = R nˆ F , (8)
tact law. dt i 2 dt m r2 i ij × ij
i i j
X
where ω˜ denotes the fluctuating part of the angular ve-
i
SLLOD equations for simple shear flow locityandwehaveinsertedthemomentofinertiaofcon-
stant density disks in two dimensions.
Lees-Edwards (LE) boundary conditions permit us to Equations(7)and(8)nowgiveanexactrepresentation
prescribethedeformationofamaterialbycontrollingthe of simple shear flow for a frictional granularmaterialar-
positions of the image cells [92]. In allofthe simulations bitrarily far from equilibrium.
presented here, we impose a constant strain rate γ˙ so The primary interest of this procedure is that it per-
that a grain at position y has an average velocity of γ˙y mits us to simulate a sheared granular material with a
in the x-direction (see Figure 1). homogeneous shear rate. Experimental procedures, e.g.
7
in a Couette cell, do not guarantee that the strain rate 300
is homogeneous: the existence of walls induce a non- γ˙ =10−2
2
uniformity of the flow and possibly localization of the −200
˙γ
deformation. Our protocol grants direct access to the p
rheology of the granular material in a self-averaging sit- 100
uation.
0
0 0.5 1 1.5 2 2.5 3 3.5 4
shear strain
Macroscopic Quantities
300
γ˙ =102
WedefinethestresstensorΣαβ viaCauchy’sequation: 2−200
˙γ
p
100
d
φ vα +φ vβ ∂ vα = ∂ Σαβ +fα , (9)
dth i h i βh i − β ext
0
0 0.5 1 1.5 2 2.5 3 3.5 4
where vα istheaveragevelocityintheαdirection(aver- shear strain
h i
agedoverallgrains),φisthemassdensity,andf isthe
ext
external force per volume. This equation simply states FIG. 4: Raw data of the normalized pressure as a function
oftotalshearfortwofrictionlessgranularmaterialwithpack-
thatthetotaltimederivativeoftheaveragevelocity(left
ing fraction 0.8. Data from simulations with different shear
hand side) is proportionalto the divergence of the stress rates γ˙ are shown. The top plot corresponds to γ˙ = 10−2
tensor,plusanyexternalforces. Cauchy’sequationgives and the bottom to γ˙ = 102. The pressure is normalized by
adefinition ofthe stresstensorfromwhichthereexists a γ˙2 which collapses the two data sets on to one master curve
proceduretoderivethe functionalformofthe stressten- (i.e. with this rescaling the top and bottom traces appear
sor. This is called the Irving-Kirkwood derivation [94] essentially identical, up to numerical noise), as predicted by
and it yields theinvariance for hard spheresystems.
ΣαβV = m v˜αv˜β + (R +R )nˆαFβ, (10)
i i i i j ij ij IV. TEST OF LOCAL RHEOLOGY
i i>j
X X
where v˜ is the fluctuating velocity of graini determined Wefirststudytheaveragerelationbetweenpressurep,
i
by v˜α = vα vα , and V is the volume of the granular shearstress σ, shear rate γ˙, andgranulartemperatureT
i i −h i i
material (or area in two dimensions). This symmetric insimpleshear,usingtheperiodicandtranslationallyin-
stress tensor can be written in terms of three variables: variantLEcell. Nextwecomparetheseresultswithdata
the pressure p, shear stress σ, and first normal stress obtained for the granular flow down an inclined plane.
difference N defined as
1
Σαβ = p(1+N1) −σ , (11) Simple Shear
σ p(1 N )
1
(cid:18) − − (cid:19)
Preliminary test
where the signs are chosen so that shear stress and pres-
sure are positive in our conventions.
In all of the simple shear simulations presented here
The granular temperature T is measured as
we have simulated 2500grains in a square primitive cell,
2 2 although we have conducted a limited number of simu-
1 1
T = v2 v + R2ω2 R ω lations with up to 10000 grains to ensure the accuracy
i − i! 2 i i − 2 i i! of our observations. Because the contact dynamics al-
i i i i
X X X X
(12) gorithm induces some amount of numerical noise, the
which is proportional to the kinetic temperature from motion of a collection of grains driven at different shear
equation(5), withthe factorof1/2comingfromthe mo- ratesisnotexpectedtoreproduceexactlythesamephase
ment of inertia calculated for constant density discs. spacetrajectory. InFigure4weshowrawdataofthenor-
Lastly, for all of the numerical data that will be pre- malized pressure pγ˙ 2 as function of shear strain (strain
−
sented,wequantifythedensityofthesystembyitspack- ratemultipliedbytime),atapackingfractionof0.8with
ingfractionν. Thepackingfractionisdefinedasthearea no friction and at two different values of the shear rate.
occupied by grains divided by the total area of the sys- According to the invariance in Newton’s equations
tem. Inour simulations,packing fractionis proportional pγ˙ 2 should be independent of γ˙, and this behavior is
−
to both the mass density φ (πφ = 4ν) and the number confirmed by the measurements in Figure 4. Although
density n (nπ R2 =ν). the shear rates in the two plots differ by a factor of 104,
h i
8
1 103
p
/
σ0.5
0 102
0 0.2 0.4 0.6 0.8 1
6
T
√4
˙γ/ 101
2
0
0 0.2 0.4 0.6 0.8 1
0.1 100
p
/
T0.05
0 10−1
0 0.2 0.4 0.6 0.8 1 0.5 0.55 0.6 0.65 0.7 0.75 0.8
shear strain packing fraction
FIG. 5: Invariant quantities σ/p (top), γ˙/√T (middle), and FIG. 6: Steady state values of pγ˙−2 (squares), σγ˙−2 (dia-
T/p (bottom) as a function of shear strain for a frictionless monds),Tγ˙−2 (circles),andσ/p(stars)asafunctionofpack-
granular material at packing fraction 0.8. ing fraction for a frictionless granular material.
Liquid-solid transition
thenormalizedpressureisvirtuallyidenticalforbothsys- InFigure6wepresentthesteadystatevaluesofpγ˙ 2,
−
tems. Interestingly, not only do the steady state values σγ˙ 2, Tγ˙ 2, and σ/p in simple shear for a range of high
− −
show no shear rate dependence, but the initial transient packing fraction systems that we have studied, at zero
is virtually identical for both values of γ˙. The invari- friction. Althoughthereisrelativelylittlechangeinthese
ance in Newton’s equations also predicts that σγ˙ 2 and quantities for small packing fraction, for packing frac-
−
Tγ˙ 2 areindependent ofγ˙. Forallofthe simulationswe tionslargerthat0.75thereisalargeincreaseinthevalues
−
havecarriedoutthesepredictionsfromtheinvarianceare of the stresses and granular temperature. Additionally,
upheld–althoughnumericalnoiseoftendisruptsthe per- the functional form of the stresses and granular temper-
fect invariance for large values of shear strain, we see no aturechangesfromapproximatelyexponentialtoafunc-
changeinthe steadystate valuesofnormalizedpressure, tion that grows faster than an exponential at ν 0.75.
≈
shearstress,orgranulartemperature as the shearrate is Alam and Luding [95] have measured the steady state
varied at constant density. These tests ensure that the valuesofthefirstnormalstressdifferenceN ,definedvia
1
simulations we conduct respect the invariance in New- equation(11), in simple shearflow using a monodisperse
ton’s equations for hard spheres. collectionofgrains. Indiluteflows,anon-vanishingvalue
forN resultsfromcollisionaltermsandtheanisotropyin
1
Foragranularmaterialcharacterizedbyitspressurep, the distributionofvelocitiesatBurnettorder[96]. Alam
shear stress σ, temperature T, and strain rate γ˙, we can and Luding reported that N1 becomes negative at the
construct three independent invariant quantities: σ/p, onsetof crystallization. We have thus measuredthe first
γ˙/√T, and T/p. In Figure 5 we show values of these normal stress difference in our system in order to ensure
three independent invariant quantities as a function of the absence of crystallization and as a signature of the
shearstrainforafrictionlessgranularmaterialatpacking breakdownofkineticeffects. TheobservationinFigure7
fraction of 0.8. For all quantities, steady flow is reached that N1 > 0 for all packing fractions in our system is
byashearofapproximately0.5,andwewillsubsequently consistent with the observation that our system remains
provide stationary data by time-averaging our measure- amorphous. ThedecayofN1 overallpackingfractionsis
ments between strains of 2 and 10. The values of σ/p consistent with the idea that kinetic effects become less
and T/p fluctuate much more than γ˙/√T. This is due important as the packing fraction increases, especially
to the fact that σ and p depend on the forces between after ν 0.75 when N1 begins to quickly decay.
≈
grains, which are highly fluctuating in the hard sphere Inour numericalsimulations,we expectthe collisional
limit. In our simulations with 10000 grains we observe contributions to the stress tensor to be negligible. A de-
thatthefluctuationsdecreasewhilethe averagevaluere- tailedstudyofthecrossoverbetweencollisionalandnon-
mains constant. This suggests that in the limit of large collisionalregimeswillbe the topic ofafuture work[97].
system size, the fluctuations would disappear. For now we rely on the observation that the kinetic ef-
9
1
ν0.8
0.6
10−1 0 10 20 30 40 50
1
N1 T0.5
0
0 10 20 30 40 50
0.4
T
√
10−2 /
˙γ0.3
0.2
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0 10 20 30 40 50
packing fraction height (y)
FIG. 7: First normal stress difference N as a function of FIG. 8: Profiles of packing fraction (top), granular tem-
1
packing fraction for a frictionless granular material. N1 > 0 perature (middle), and γ˙/√T (bottom) as a function of the
for all packing fractions suggests that the granular materials height(y)inthepile,measuredingraindiameters,foranon-
are amorphous, and thedecay of N1 as a function of packing frictional granular material at a 12◦ incline.
fraction suggests that non-kinetic effects dominate at high
packing fraction.
p, σ, γ˙, and T– and these lead to three independent in-
variant quantities γ˙/√T, σ/p, and T/p. We observe in
fects, as probed by the first normal stress difference, de-
oursimulationsthatalloftheseinvariantsareconstantin
cayclosetojamming. Wealsorefertopriorworkswhich
the bulk of the incline flow. Therefore it is legitimate to
support the same conclusion [77, 78].
compare these constant values with the constant values
obtained from simple shear simulations. In Figure 9 we
present how the constant values of σ/p, T/p, and γ˙/√T
Flow down an inclined plane
in the bulk of the flow depend on packing fraction and
comparewithourresultsfromthesimpleshearcell. The
Inthesimpleshearcellthedensityandstrainratewere
fact that data from different shear flows fall onthe same
specified, while stresses and granular temperature were
curvesisremarkableandsuggeststhatonetheoryshould
measured. We now focus on flow down aninclined plane
beabletodescribesimpleshearandbulkinclinegranular
which provides a complementary situation. For incline
flow.
flow, stresses are specified by the choice of an angle of
Interestingly,thedatafromdifferentflowsdonotover-
inclination and by the gravitationalfield. Then the pro-
lap over a large interval of packing fraction: the flow
files of velocity and velocity fluctuations are measured
down an inclined plane provides values at higher values
and grant access to profiles of strain rate and granular
ofpackingfractionthanthesimpleshearcell. Thisisdue
temperature.
to (i) steady flows down an incline plane are more eas-
We report in Figure 8 the packing fraction, granular
ily reachedfor lower inclinations, hence higher densities;
temperature, and γ˙/√T as a function of height for the
and (ii) the simple shear deformation is more difficult to
steadyflowofanon-frictionalgranularmaterialatanan-
integratenumericallyathigherdensities,becausethepe-
gleof12 ,withatotalheightofapproximately50grains.
◦
riodiccellinducesadditionalconstraintsthattheContact
We have conducted simulations with heights rangingbe-
Dynamics algorithm manages with difficulty. Neverthe-
tween 25 and 100 grain diameters to ensure that our re-
less, our use of two different configurations grantsaccess
sults do not depend on the size of the system.
to a broadrangeof γ˙/√T, σ/p, andT/p;andthe sets of
Weobservethat,inthebulkcentralregionoftheflow,
data are consistent with the existence of a unique, local
the packingfractionprofileis uniform,the granulartem-
relation between them as apparent in Figure 9.
perature is linear, and γ˙/√T is constant. These obser-
vations hold in our simulations for all angles where the
granular flow reaches a steady state. We only use data
from these steady state flows in this paper. These ob- Origin of a Local Rheology
servations are consistent with previous observations by
Silbert et al [41, 42]. The excellent agreement between the two sets of data
There are four quantities of interest for incline flows– challenges the belief that non-local effects arise in dense
10
0.7
0.4
p
/
σ
0.2 0.6
0
0.75 0.8 0.85 0.5
T
√0.6 p0.4
˙γ/0.4 σ/
0.2 0.3
0
0.75 0.8 0.85
0.2
p 0.1
/
T 0.1
0.05
0 0
0.75 0.8 0.85 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
packing fraction γ˙/√T
FIG. 9: The values of σ/p (top), γ˙/√T (middle), and T/p FIG. 10: σ/p plotted against γ˙/√T in steady state simple
(bottom)plottedasafunctionofpackingfraction. Datafrom shear flow (circles) and incline flow (squares). Filled sym-
simpleshearflow(circles) andflowdownanincline(squares) bols correspond to frictional granular materials with coeffi-
match on the same curves. This suggests that there is a lo- cient of friction µ = 0.4 and open symbols correspond to
cal rheology that is independent of the particular shearing non-frictional granular materials with µ=0.
geometry.
dense granular material ℓ/ R << 1, the validity of the
h i
granularflow. The datain Figure9supports averycon- hydrodynamic limit should hold over moderate heights.
servative opinion: the motion of the grains decorrelates As we will see, our theory provides constitutive equa-
beyond some finite length scale, in accord with the fast tions which relate the static friction coefficient σ/p to
decay of the pair correlation function (see Figure 3). the ratio γ˙/√T. Anticipating the following sections, we
present in Figure 10 a plot of σ/p versus γ˙/√T, with
As a consequence, in the bulk of the flow down an in-
data from both the simple shear and incline flow geome-
cline, it is possible to view layers of granular materials
tries. As we will argue, such a plot is expected to be the
as effective simple shear cells. Such a layer of granular
granular counterpart of a stress versus strain rate plot
materialatheighty respondsessentiallyas ifit wascon-
for a glassy material [3], with the shear stress normal-
finedinasimpleshearcell,intheabsenceofbodyforces,
ized by the pressure and the strain rate normalized by
withsustainedexternalstressesσ(y)andp(y). Ofcourse
the granular temperature. We see here that the analogy
theinvarianceinNewton’sequations,whichholdsexactly
is striking: when rescaledproperly the granularmaterial
for the simple shear cell, is slightly broken by the grav-
presentstypicalfeaturesofnormalyieldstressfluids[98].
itational force field. However, deep in the bulk of the
Forlargevaluesofnormalizedstrainrate,thenormalized
flow, large confining stresses eventually dominate over
stress is proportional to the normalized strain rate. For
the gravitationalfield. This approximate invariance suf-
smallvaluesofthenormalizedstrainrate,thelinearrela-
fices to predict that Bagnold’s scaling must hold for the
tionshipnolongerholdsandthereisayield(normalized)
bulk regions of incline flows and explains the numerical
stress at zero (normalized) strain rate.
data of Silbert and coworkers[41, 42].
Hydrodynamic equations are expected to arise when
thegradientsofmacroscopicvariablesbecomesmallcom-
V. CONSTITUTIVE EQUATIONS
pared to the macroscopic variables themselves. What is
surprising,however,isthatthelocalityoftherheologyis
observable for the moderate heights that we can access Severalinterestingresultsstemfromtheprecedingob-
inournumericalsimulations. Webelieve thatthe reason servations: (i) for a given density, only invariant quanti-
why the hydrodynamic regime seems to be observed at ties are relevant to describe the state of a granular flow;
accessible scales is the following: (i) the relevance of hy- (ii)wheninvariantquantitiesareconsidered,thegranular
drodynamics limit depends on how a macroscopic quan- materialdisplaysa“normal”rheologyofayieldstressliq-
tity Ψ compareswith ℓ Ψ, whereℓ is a length-scale;(ii) uid; and (iii) the rheology measured in the Lees-Edward
∇
ℓisnotthesizeofthegrainoraclusterofgrains–itisthe cellextrapolatestotherheologymeasuredfrombulkdata
sizeoftheregionsampledbyagivengrainperstrainunit, of incline flow.
which is on the order of the mean-free path. Since for a These observations foster our hope to construct local