Critical packing in granular shear bands

0.19 MB·

by  S. Fazekas

Critical packing in granular shear bands
S. Fazekas1,2, J. T¨or¨ok1 and J. Kert´esz1
1Department of Theoretical Physics,
2Theoretical Solid State Research Group of the Hungarian Academy of Sciences,
Budapest University of Technology and Economics (BME),
H-1111 Budapest, Hungary
(Dated: October13, 2006)
7 In a realistic three-dimensional setup, we simulate the slow deformation of idealized granular
0 media composed of spheres undergoing an axisymmetric triaxial shear test. We follow the self-
0
organizationofthespontaneousstrainlocalizationprocessleadingtoashearbandanddemonstrate
2
theexistenceofacriticalpackingdensityinsidethisfailurezone. Theasymptoticcriticality arising

n fromthedynamicequilibriumofdilationandcompactionisfoundtoberestrictedtotheshearband,
a while the density outside of it keeps the memory of the initial packing. The critical density of the
J shearbanddependson friction (andgrain geometry) andin thelimit of infinitefriction it definesa

3 specificpackingstate, namely the dynamic random loose packing.

] PACSnumbers: 45.70.Cc,81.40.Jj
t Keywords: granularcompaction, stress-strainrelation
f
o
s
. I. INTRODUCTION While this hypothesis forms the basis of many contin-
t
a uum constitutive models of soil mechanics since decades
m
[11],ageneralmicromechanicaltheoryofshearbandfor-
Packing density of different particulate systems is of
- mation and of the involved criticality is still missing.
d main interest for scientific fields including, but not lim-
Progress, needed in order to deepen our understanding
n ited to, suspensions, metallic glasses,molecular systems,
o and granular materials. In three dimensions, for identi- ofthe criticalstateinshearbands,canbeexpectedfrom
c the remarkable development of experimental techniques
calspheres,theface-centeredcubic(FCC)packingisthe
[ (including Computer Tomography [12, 13] and measure-
maximum possible [1]. This fills the space with a vol-
ments in microgravity [14, 15]) and of simulation tech-
3 ume fractionofπ/(3√2) 0.74. Randomarrangements
v ≈ niques which become increasingly efficient as computa-
havemuchlowerdensities [2]. Differentexperiments and
3 tional power grows [16].
computersimulationsrevealedthatthelargestobtainable
7
volume fractionofa randompacking ofidentical spheres Shearing of granular materials has been investigated
4
4 is around 0.64. This is known as the random close pack- in many different geometries and specially designed lab-
0 ing (RCP) limit. A mathematical definition of this limit oratory tests (for recent results see [12, 14]). Such ex-
6 can be given through the concept of maximally random perimentalstudiesrevealedcomplexlocalizationpatterns
0 jammed state [3, 4]. andpresentedevidencefortheexistenceofacriticalpar-
/ ticle density inside the shear bands. The importance of
t ReachingtheRCPlimitneedscarefulpreparation(e.g.
a computer simulations is enhanced by the fact that they
m tapping and compression). If glass or marble beads are
makepossiblestudieswhicharedifficulttocontrolinex-
simplypouredintoacontainerthevolumefractionisusu-
- periments (e.g. friction dependence) and they facilitate
d ally only around 0.6. A random loose packing (RLP) at
themeasurementofhardlyaccessiblequantities(e.g. vol-
n its limit of mechanical stability obtained by immersing
ume fraction inside the shear bands).
o spheres in a fluid and letting them settle has a volume
c fraction of 0.555 [5]. The volume fraction of RLPs ob- Thecriticaldensity,innumericalstudies,isoftenstud-
:
v tainedwithdifferentmethods(bothexperimentalandnu- ied only in special conditions when shearing extends to
i merical) show that this packing state is less well defined the whole volume of the samples. This allowed for dis-
X
than the RCP limit. Attempts made in order to relate cussing the criticality based simply on global behavior
ar RLP to rigidity percolation [5] and to critical density at (e.g.dilatancy). Withoutreferencetoshearbands,many
jamming of an assembly of (infinitely) roughspheres [6], qualitativeeffectswerealreadypointedoutinbothmath-
are to be mentioned. ematical models [17, 18] and simulations [18, 19]. How-
ever,suchstudiesneglectedtheinvolvedlocalizationphe-
Already in 1885, Reynolds noted that dense granular
nomena inevitable in real situations and disregardedthe
samplesdilateduringslowdeformation[7]. Ontheother
self-organizing manner in which the packing state of the
hand, it is wellknownthat loosegranularmaterialsden-
shear bands is usually formed.
sifyinsuchaprocess[8,9]. Underslowshearthestrainis
usually localized to narrow domains called shear bands. A principal parameter which controls the dynamic
As it was firstsuggestedby Casagrande[10], it is tempt- equilibriumbetweendilationandcompactioninfully de-
ing to assume that in these failure zones the system self- velopedshearbandsisthefrictionbetweenthegrains. In-
organizes its packing density to a critical value indepen- tuitively,asystemoffrictionlessgrainscanbeshearedat
dent of the initial packing state of the material. a large packing density (close to the RCP limit) because
2
µ0 0.8 0.5 0.3 0.2 0.1 0.0
η0 0.555 0.562 0.578 0.599 0.621 0.641
TABLE I: Volume fraction η0 of samples prepared with dif-
ferent coefficients of friction µ0. With each µ0 we prepared 2
samples havingthesame η0 within 0.2% relative error.
The numericalvalues are chosento realize the hardest
material that we could safely simulate with the minimal
dampingpreservingthenumericalstabilityofthecalcula-
tions. With the abovestiffness anddamping coefficients,
FIG. 1: (Color online) a) Grains placed between two hori- the inverse of the average eigenfrequency of contacts, in
zontal platens and surrounded by an elastic membrane were both normal and tangential direction, is more than one
subjected to a vertical load and a lateral confining pressure. orderofmagnitude largerthanthe usedintegrationtime
b)Themembranewasmodeledwithoverlappingspheres[20]. step ∆t = 10−6s. This assured that the noise level in-
duced by numerical errors is kept low.
With relatively small samples (made up of 27000 par-
thegrains(underslowshear)caneasilyrearrangeincom-
ticles) but in a realistic geometry we have succeeded to
pact configurations. At large friction the rearrangement
reproduce shear band morphologies [20, 23] known from
of the grains is hindered by friction, consequently the
experiments [12, 14]. In order to study the criticality
packing density of the shear bands is expected to define
of these shear bands, we prepared homogeneous initial
a low density state close to the RLP limit.
configurationsof differentvolume fractions using the de-
The aim of this Paper is to study the emergence of a
position method described in [20].
critical packing state in sheared granular media and to
We used a particle distribution similar to those en-
present its relation to shear bands as well as its depen-
countered in experimental studies of idealized granular
dence on friction.
materials. Our particles are spherical, they have equal
mass density (2.5 103 kg/m3), equal friction coefficient,
·
and their diameters are set according to a narrow Gaus-
II. SIMULATIONS
sian distribution with mean d = 0.9 mm and standard
deviation of 2.77%. The prepared cylindrical samples,
We investigate numerically an axisymmetric triaxial havingdiameter D =23.3d, consistedof20000to 27000
shear test (see Fig. 1). This consists of the slow com- sphericalgrainsasrequiredbya prescribedpackingden-
pression of a cylindrical sample enclosed between two sity and the H 2.2 D geometricalconstraint,where H
end platens. The sample is surrounded by an elastic is the height of≈the samples.
membraneonwhichanexternalconfiningpressureisap- Initially the particles were placed randomly in a tall
plied. The end platens are pressed against each other in cylinder (about 3 times taller than H). They were given
a strain controlled way. The upper platen is allowed to smalldownwardsvelocitiesinsuchawaythattheyallcol-
tilt. In certain conditions a planar shear band is formed lided approximatelyat the same time. The upper platen
[12]. Using different initial packingdensities and friction was pressed on top of the packing to hold it together.
properties of the grains as well as identifying the grains This method provides an efficient way to produce a ho-
inthefailurezonesmakesitpossibletostudythecritical mogeneous random packing. The volume fraction of the
packing density inside the shear bands. prepared samples could be controlled in the full RLP to
Thesimulations,whicharegoingtobepresentedhere, RCP range (see Tab. I) by varying the coefficient of fric-
are based on a standard three-dimensional Distinct El- tion µ which was applied during this phase.
0
ement Method (DEM) [21]. We implemented the Hertz After preparation, the friction coefficient of the parti-
contact model [22] with appropriate damping combined cles was set to a new value µ independent of µ . During
0
with a frictional spring-dashpot model [16]. The normal thesimulations,similarlyto[20],wecompressedthesam-
Fn and the tangential Ft forces are calculated as plesverticallyatzerogravityand0.5kPaconfiningpres-
sure. The bottom platen was fixed. The upper platen
F = κ δ3/2 γ δ1/2v (1)
n n n − n n n moved downward with a constant velocity, inducing an
Ft = κtδt γtvt (2) axialstrainrateof20mm/s. Duringcompressiontheup-
−
perplatencouldfreelytiltalonganyhorizontalaxiswith
where κ =106N/m3/2, κ =104N/m, γ =1Ns/m3/2, rotationalinertia I=10−7kgm2.
n t n
and γ =1Ns/m are normal and tangential stiffness and The lateral membrane surrounding the sample was
t
damping coefficients, δ and δ are normal and tangen- modeled with approximately 15000 identical, overlap-
n t
tial displacements, and v and v are the normal and ping,non-rotating,frictionalspheresconnectedwithelas-
n t
tangential relative velocities. tic springs. The stiffness of the springs was set to
3
κ =0.5N/m. Thispreventedtheparticlesfromescaping
s
bypassingthroughthemembrane. Themembraneparti-
cles were initially arrangedin a triangular lattice (Fig. 1
(b)). The confining pressure was applied on the triangu-
lar facets formed by the neighboring“membrane nodes”
as described in [20] (see also [24, 25, 26]).
It is worth mentioning that Cui and O’Sullivan [26]
haveintroducedatechniquewhichspeedsupcalculations
by computing only a section of the cylindrical sample.
This allowsforlargersamplesbutrequiresthatthe sym-
metryofthesystemiskeptduringcompression,andthus
eliminatesthepossibilityofsymmetrybreakingstrainlo-
calization,whicharisesspontaneously[20]iftiltingofthe
upper platen is not suppressed (see Fig. 1 (a)).
We have executed several simulation runs. The grain-
platenandgrain-membranecontactswerecalculatedsim-
ilarly to grain-graincontacts including the friction prop-
erties. The samples of different initial volume fractions
(see Tab. I) were first compressed using the same coef-
ficient of friction µ = 0.5. Later, we compressed the
densest samples (η = 0.641) with 10 different friction
0
coefficients µ 0, 0.1,...0.9 . For each set of parame-
∈{ }
ters, two simulation runs were executed using specimens
prepared with different random seeds.
During compression,we measuredlocally the shearin-
tensity S and the volume fraction η. The regular trian-
gulation of the spherical grains [27, 28] was used to de-
fine these quantities. The local volume fraction is given
by the ratio of the volumes of a grain and its regular
Voronoicell. The localshear intensity is calculatedfrom
the macroscopic strain tensor derived from particle dis-
placements [20,29]. Usingthe eigenvaluesε ofthis ten-
k
sor, we defined the local shear intensity as
(cid:12) 1 (cid:12)
S =mkax(cid:12)(cid:12)εk− 3Xεl(cid:12)(cid:12). (3)
(cid:12) l (cid:12)
(cid:12) (cid:12)
To overcomefluctuationsdue to randompackingandre-
arrangements, we calculated spatial averages up to 3rd
order neighbors along the regular triangulation.
FIG. 2: (Color online) Example of shear intensity and vol-
III. IDENTIFICATION OF HIGH SHEAR
INTENSITY REGIONS ume fraction histogram maps (η0 = 0.641, µ = 0.5). The
shear intensity is measured in arbitrary units. On (a, b) the
occurrences(scaled to[0,1]) are encodedwith thecolor scale
Strain localization in dense and loose samples shows shown at the top. On (a) the white curve marks the shear
substantialdifferences[12]. Indensesamplesshearbands intensitythresholdSTH. On(b)itmarkstheshearbandvol-
are usually formed after a short plastic deformation and ume fraction ηSB. The dotted and dashed vertical lines on
insidethemthelocalpackingdensityislowerthaninthe (a, b) at 4.4% and 16% axial strain mark the position of the
bulk (i.e. the regions outside of the shear bands). Since histograms (c,d) and (e, f), respectively.
denser parts are more stable, the position of the shear
bands remains unchanged for the whole duration of a
shear test. Contrary, in loose samples the shear bands The general algorithmic identification of failure zones
have a slightly higher packing density than the bulk and based on geometric methods is difficult, especially, re-
hence the position of the shear bands is likely to change garding the identification of the non-persistent shear
andtomovearoundthewholesample. Thisleadstomore bands of the loose samples. Nevertheless, based on the
orlesshomogeneoussampleswithlocalpackingdensities calculated local shear intensity, individual grains can be
close to the packing density of the shear bands. categorized to be part of the failure zones or the bulk
4
providingthata goodenoughthresholdseparatingthese
two classes can be found. A histogram technique seems
to be a perfect candidate for this.
Fig. 2 presents histograms of the local shear intensity
and local volume fraction. It shows the main aspects
of the strain localization process observed in one of our
simulations executed with a sample having η = 0.641
0
andµ=0.5. Atthebeginningofthetest,uptoapproxi-
mately6%axialstrain,inshearintensityhistograms(see
Fig. 2 (a, c)), we find a single peak at medium S. This
means that almostallparticlesrearrangesimultaneously
and consequently the sample experiences a more or less
plastic deformation. This is underlined by the fact that
thelocalvolumefractionhasjustonestrongpeak(Fig.2
(d)), indicating that the sample is still homogeneous.
At higher (>6%) axial strain a shear band is formed.
This is localized to a planar failure zone of width of ap-
proximately10particlediametersandischaracterizedby
much higher S than the bulk. In the bulk the shear in-
tensityfluctuationsaresmall,whilethesefluctuationsare
largein the shear band. Consequently,inshear intensity
histograms (Fig. 2 (e)), we find a narrow peak at low S,
corresponding to the bulk, and a wide peak at high S,
corresponding to the shear band. The volume fraction
histogramdoesalsobecomemorestructuredshowingev-
idenceofanon-homogeneousmaterial. Thenarrowpeak FIG. 3: (Color online) Volume fraction measured in high
at low volume fraction corresponds to the shear band shearintensityregions(a) andglobally (b)as functionof the
while the bulk produces a much wider distribution at a axial strain. The different lines correspond to different ini-
higher volume fraction (Fig. 2 (f)). tial packing densities η0 (see Tab. I) decreasing from top to
bottom. The two panels usethesame notations.
Motivated by this separation, we computed a shear
intensity threshold S , which could be used to define
TH
two classes of shear intensity values (low and high) and
byOtsu’smethod, whichfalls onthe middle ofthe peak,
to classify the grains accordingly into shear band and
is not physically relevant. However, as the sample is ho-
bulk. ForthisweusedOtsu’sthresholdselectionmethod
mogeneous, the selected samples in high shear intensity
[30] described in Appendix A. This histogram technique
regionsstillgivethevolumefractionwhichisclosetothe
minimizes the within-class variance and maximizes the
average volume fraction of the whole sample. This can
separation of classes, and thus gives an ideal solution to
be verified on Fig. 2 (d).
our problem. We have also tested another threshold se-
lectionmethodmodelingthe histogramswiththe sumof
twoGaussianfunctions,however,numericallythisproved
to be less stable and less reliable. IV. RESULTS
The condition S > S , made it possible to identify
TH
the grains in high shear intensity regions – where shear As expected [12, 14], we found that due to strain lo-
bandsemerge–andthustheaveragevolumefractionηSB calizationatthe end ofthe sheartests the globalvolume
of these regions could be calculated. The local packing fraction of the samples is not equal to the packing den-
densityinshearbandsisfoundtohavesmallfluctuations sity of the high shear intensity regions and thus global
and to give a peak in the volume fraction histograms. measurements cannot be used to characterize the prop-
This coincides with ηSB (see Fig. 2 (b, f)), giving a self erties of failure zones. The behavior of both dense and
validation of the method. loose samples demonstrates that in the shear bands, the
Let us note that even if Fig. 2 (e) does not suggest a initial packing conditions are canceled and a criticalvol-
clean separationofthegrainsintothosewithinandthose ume fraction η is reached in a self-organizing manner
c
outside the shear bands – i.e. there is no shear intensity independently of the initial density of the tested granu-
gapbetweenthetworegimes–thisisnotcrucialforη . lar specimens (Fig. 3 (a)).
SB
Adding an artificial random noise of 10% to STH does The criticality is found to be restricted to the shear
influence the resulting ηSB only within 0.5%. bands. The global volume fraction ηg calculated from
We also mention here, that before strain localization the total volume of the samples does not converge to
takes place the shear intensity histograms have only one η (Fig. 3 (b)). This behavior is expected to be more
c
peak (see Fig. 2 (c)). In this case, the threshold given pronounced on larger systems. The dense samples are
5
bands. Wemeasuredthepackingdensityη insidethese
SB
failure zones and we found that in fully developed shear
bands η approaches a critical value η independent of
SB c
the initial density of the samples. This in agreement
with Casagrande’s [10] observation made for sandy soils
seven decades before and also with recent experiments
[12, 13, 14, 15] and numerical studies [17, 18, 19].
Rothenburg and Kruyt[18] obtained similar results in
two-dimensionalsimulations of biaxial shear tests. They
havepresentedatheoryoftheaveragecoordinationnum-
ber of sheared granular media and derived a law for its
evolution during slow deformations. Analyzing the re-
lationship between volume fraction and average coordi-
nation number, they conclude that a proper characteri-
FIG. 4: (Color online) Critical packing of shear bands as
zation of granular media undergoing shear deformation
function of friction measured at 20% axial strain. The fitted
curveshowsηc(µ)=ηc∞−(ηc∞−ηc0)exp(−µ/µ0c),whereηc0 = should be based on packing density.
0.637, ηc∞ =0.578, and µ0c =0.23. We have shown that the criticality is restricted to the
shearbandsandglobalmeasurements(suchasdilatancy)
are unsuitable for the investigation of the properties of
characterized by η > η . This demonstrates that the sheared granular materials in realistic situations, where
g c
dilatancy [7] is concentrated to the shear bands. Con- strain localization is inevitable. To our knowledge, it is
trary, for loose samples η < η , however, the specimen the first time that the critical packing density of shear
g c
isonlyslightlylooseroutsidetheshearbands. Thisgives bands was evidenced based on simulations of a realistic
a direct proof of shear induced compaction [9]. three-dimensionalsetupandspontaneousstrainlocaliza-
As we could see, at given friction, the critical pack- tion revealing the self-organizing manner in which the
ing of shear bands is independent of the initial den- packing state of the shear bands is developed.
sity. It is, however, a further question whether it de- We have further shown that η depends on the coeffi-
c
pends on friction. We studied quantitatively this effect cientoffrictionµandinthe limitµ itconvergesto
for samples with initial density η =0.641. At µ=0 the avalueη∞,whichwehavecalculated→wi∞thintheaccuracy
0 c
high shear intensity regionshave a large volume fraction of our simulations. The found η∞ defines a low density
c
η0 η (0)=0.637 0.002,whichisonly slightlysmaller dynamic random loose packing (DRLP) state, which is
c ≡ c ±
than η , showing that frictionless granular systems can characteristic to the dynamic equilibrium between dila-
0
be sheared at densities very close to the RCP limit. tionandcompactionintheshearbandsanddependsonly
In frictional systems, the volume fraction of the fully on the geometry of the grains. Based on the underlying
developedshearbandsissubstantiallylowerthantheini- mechanism, we argue that the asymptotic packing state
tial volume fraction (see Fig. 4). With increasing µ, this of shear bands differs from the static RLP limit.
decreases and converges to a limit which we estimate to This result should be also compared with the findings
ηc∞ ≡ limµ→∞ηc(µ) = 0.578±0.003 based on the expo- presentedrecentlybyZhangandMakse[6]regardingthe
nential extrapolation of our data. critical density of granular materials at jamming tran-
This limit volume fraction depends only on geometry sition. The importance of these findings lie in the fact
factors such as shape and size distribution of the grains that jamming is a basic concept which through a uni-
andischaracteristictothedynamicequilibriumbetween fying phase diagram [31, 32] connects granular matters
dilation and compaction developed in a self-organized with a variety of other systems including dense particu-
manner through strain localization. Based on the value latesuspensionsandeffectssuchasdivergingviscosityat
of η∞ and this latter aspect, the corresponding asymp- a maximum packing fraction [33, 34].
c
totic state – which we refer to as the dynamic random In quasistatic limit, Zhang and Makse [6] reported a
loosepacking (DRLP)–shouldbedistinguishedfromthe monotonous decrease of the critical packing density as
static RLP limit. function of friction. For low friction, the density of the
shear bands found in our simulations is lower than at
jamming found by Zhang and Makse [6], while at high
V. CONCLUSIONS AND DISCUSSION
frictionthesituationisreversed. Thisindicatesanatural
separation of low and high density regions with possibly
We have presented distinct element simulations of ax- different mechanisms of dissolving the jammed state.
isymmetric triaxial shear tests at zero gravity and low As a final remark let us note that our results are de-
confining pressure. Due to spontaneous strain localiza- rivedforidealizedgranularmaterialscomposedofspheres
tion, shear bands were formed. Using a histogram tech- having a narrow size distribution. It is well known that
nique, we identified the grains in high shear intensity re- for non-spherical grains and wide size distributions the
gions,whichatlargeaxialstrainscoincidewiththeshear packingefficiencyincreases[35], whichshouldbe alsore-
6
flected in the packing density of the shear bands. This The within-class variance
couldbethe reasonwhyexperimentalresultsonsandre-
vealsmallervolumefractionsinshearbands[12,14]than
σ2 (t)=q (t)σ2(t)+q (t)σ2(t) (A7)
the values found in our simulations. W 1 1 2 2
is an inverse measure of the compactness of classes. The
VI. ACKNOWLEDGMENTS
between-class variance
This research was carried out within the framework
σ2(t)=q (t)q (t) µ (t) µ (t) 2 (A8)
of the“Center for Applied Mathematics and Computa- B 1 2 (cid:0) 1 − 2 (cid:1)
isameasureoftheseparationofclasses. Itiseasytoshow
tional Physics” of the BME, and it was supported by
thatσ2 (t)+σ2(t)=σ2. Otsu[30]proposedtocalculate
OTKA F047259 and T049403, and the P´eter Pa´zm´any W B
anoptimalthresholdt=t byeitherminimizingσ2 (t)
program RET-06/2005. S.F. thanks D. Chetverikov for opt W
or maximizing σ2(t).
the introduction to Otsu’s method. B
Maximizing σ2(t) is easier. It can be seen that
B
APPENDIX A: OTSU’S THRESHOLD
SELECTION METHOD q2(t) = 1 q1(t), (A9)
−
µ q (t)µ (t)
1 1
µ (t) = − , (A10)
2
Otsu’s method [30] is a histogram technique known q (t)
2
fromDigitalImage Processing,where it is typically used
totransformgrayscaleimagesintotwocomponent(black
and thus
and white) images.
Let us consider a normalized histogram P(i), i.e. a
histogram with the property σ2(t)= q1(t) µ (t) µ 2. (A11)
B 1 q (t)(cid:0) 1 − (cid:1)
1
−
X P(i)=1, (A1)
i
For each candidate threshold t, q (t) and µ (t) can be
1 1
where i is the bin index. The mean µ and the variance calculated with the recursive formula
σ2 can be calculated as
µ = X iP(i), (A2) q1(t+1) = q1(t)+P(t+1), (A12)
i q1(t)µ1(t)+(t+1)P(t+1)
µ (t+1) = , (A13)
σ2 = X(i µ)2P(i). (A3) 1 q1(t+1)
−
i
Let us further consider a candidate threshold t and where q1(0)=P(0) and µ1(0)=0.
split the histogram in two parts I1(t) = {i|i ≤ t} and Bothq1(t)andµ1(t)areincreasingmonotonouslywith
I2(t)={i|i>t}. With k ∈{1,2} and t,consequentlythemaximumofσB2(t)iswell-defined,ex-
cept for degenerated cases which must be handled sepa-
qk(t)= X P(i), (A4) rately. Theoptimalthresholdtoptisgivenbythesmallest
i∈Ik(t) candidate threshold s which satisfies the equation
the mean µ (t) and variance σ2(t) of the two parts are
k k σ2(s)=maxσ2(t). (A14)
defined by the equations B B
t
qk(t)µk(t) = X iP(i), (A5) Becauseσ2 (t )+σ2(t )=σ2,themethodbothmin-
i∈Ik(t) imizes theWwitohpitn-claBss voaprtiance and maximizes the sep-
qk(t)σk2(t) = X (cid:0)i−µk(t)(cid:1)2P(i). (A6) aration of classes.
i∈Ik(t)
[1] D. A.Weitz, Science 303, 968 (2004). Phys. Rev.Lett.84, 2064 (2000).
[2] T. Aste, M. Saadatfar, and T. J. Senden, Phys. Rev. E [4] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel,
71, 061302 (2005). Phys. Rev.E 68, 011306 (2003).
[3] S. Torquato, T. M. Truskett, and P. G. Debenedetti, [5] G. Y. Onoda and E. G. Liniger, Phys. Rev. Lett. 64,
7
2727 (1990). [23] S.Fazekas,J.T¨or¨ok,J.Kert´esz,andD.E.Wolf,inPow-
[6] H. P. Zhang and H. A. Makse, Phys. Rev. E 72, 011301 ders and Grains 2005, edited by R. Garc´ıa-Rojo, H. J.
(2005), and references therein. Herrmann,and S. McNamara (Balkema, London, 2005),
[7] O. Reynolds,Philos. Mag. 20, 469 (1885). pp.223–226, cond-mat/0606720.
[8] J. T¨or¨ok, S. Krishnamurthy, J. Kert´esz, and S. Roux, [24] H. Tsunekawa and K. Iwashita, in Powders and Grains
Phys. Rev.Lett.84, 3851 (2000). 2001,editedbyY.Kishino(Balkema,Rotterdam,2001),
[9] T. L. Youd, Journal of the Soil Mechanics and Founda- pp.177–180.
tions Division 98, 709 (1972). [25] H. Sakaguchi and H.-B. Mu¨hlhaus, in Exploration Geo-
[10] A. Casagrande, Journal of the Boston Society of Civil dynamics Chapman Conference (American Geophysical
Engineers 23, 257 (1936). Union,Dunsborough,WesternAustralia,2001),pp.153–
[11] A.N.SchofieldandP.Wroth,CriticalStateSoilMechan- 155.
ics (McGraw-Hill, London,1968). [26] L. Cui and C. O’Sullivan, in Powders and Grains 2005,
[12] J. Desrues, in X-ray CT for Geomaterials, edited by edited by R. Garc´ıa-Rojo, H. J. Herrmann, and S. Mc-
J. Otani and Y. Obara (Balkema, Amsterdam, 2004), Namara (Balkema, London,2005), pp.301–305.
pp.15–41. [27] C. Lee, in Applied Geometry and Discrete Mathemat-
[13] J. Desrues, R. Chambon, M. Mokni, and F. Mazerolle, ics: The Victor Klee Festschrift, editedbyP.Gritzmann
G´eotechnique46, 529 (1996). andB.Sturmfels(AmericanMathematicalSociety,Prov-
[14] S. N. Batiste, K. A. Alshibli, S. Sture, and M. Lankton, idence, RI,1991), pp.443–456.
Geotechnical Testing Journal 27, 568 (2004). [28] H. Edelsbrunner and N. R. Shah, Algorithmica 15, 223
[15] K.A.Alshibli,S.Sture,N.C.Costes,M.L.Frank,M.R. (1996), and references therein.
Lankton,S.N.Batiste,andR.A.Swanson,Geotechnical [29] D. Daudon, J. Lanier, and M. Jean, in Powders and
Testing Journal 23, 274 (2000). Grains1997,editedbyR.P.BehringerandJ.T.Jenkins
[16] T. P¨oschel and T. Schwager, Computational Granular (Balkema, London,1997), pp.219–222.
Dynamics: Models and Algorithms (Springer, Berlin, [30] N.Otsu,IEEE Trans.Systems,Man and Cybernetics 9,
2005), and references therein. 62 (1979).
[17] M. Piccioni, V. Loreto, and S. Roux, Phys. Rev. E 61, [31] V. Trappe, V. Prasad, L. Cipelletti, P. N. Segre, and
2813 (2000). D. A.Weitz, Nature411, 773 (2001).
[18] L.RothenburgandN.P.Kruyt,InternationalJournalof [32] A. J. Liu and S. R.Nagel, Nature 396, 21 (1998).
Solids and Structures41, 5763 (2004). [33] J. J. Stickel and R. L. Powell, Annual Review of Fluid
[19] X.Zhuang,A.K.Didwania,andJ.D.Goddard,Journal Mechanics 37, 129 (2005).
of Computational Physics 121, 331 (1995). [34] G. Ovarlez, F. Bertrand, and S. Rodts, Journal of Rhe-
[20] S. Fazekas, J. T¨or¨ok, J. Kert´esz, and D. E. Wolf, Phys. ology 50, 259 (2006), and references therein.
Rev.E 74, 031303 (2006), cond-mat/0506661. [35] A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Still-
[21] P. A.Cundalland O.D. L.Strack,Geotechnique 29, 47 inger, R.Connelly, S. Torquato, and P. M. Chaikin, Sci-
(1979), and articles citing this one. ence 303, 990 (2004).
[22] L. D. Landau and E. M. Lifshitz, Theory of Elasticity
(Pergamon, New York,1970), chap.9, 2nd ed.