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.