History dependence of mechanical properties in granular systems Shio Inagaki,1,∗ Michio Otuski,2,† and Shin-ichi Sasa3 1Department of Physics, Kyoto University, Kyoto, 606-8502, Japan 2Yukawa Institute for Theoretical Physics, Kyoto, 606-8502, Japan 3Department of Pure and Applied Sciences, University of Tokyo, Tokyo, 153-8902, Japan (Dated: January 27, 2009) We study the history dependence of the mechanical properties of granular media by numerical simulations. We perform a compaction of frictional disk packings in a two-dimensional system by 9 controllingtheareaofthedomainwithvariousstrainrates. Wethenfindthestrainratedependence 0 of the critical packing fraction above which the pressure becomes finite. The observed behavior 0 makes a contrast with the well-studied jamming transitions for frictionless disk packings. We also 2 observe that the elastic constants of the disk packings depend on the strain rate logarithmically. This result providesa experimental test for thehistory dependenceof granular systems. n a J PACSnumbers: 45.70.-n;46.25.-y;83.10.Rs;62.20.D- 7 2 Granularmaterialshaveattractedtheattentionsofsci- by its time scale, which is precisely given by the inverse entists for a long time [1]. Despite the extensive efforts, ofthestrainrateoftheshrinkageinthepackingprocess. ] n our understanding has not yet been deepened enough to Byvaryingthistimescaleinasystematicmanner,weat- n constructa fundamentallaw. One majordifficulty arises tempt to observe the process dependence of mechanical - s from the process dependence of the macroscopic behav- properties such as macroscopic elastic constants. i d iors. It is rather surprising to know that the mechanical t. properties of even static granular media cannot be de- Here, it should be noted that our numerical experi- a termined from all the information of the constituents. ments are carried out with identical constituents. In or- m An example is given by a so-called stress dip, which is der to examine history dependence in an explicit way, - a local minimum of a vertical pressure distribution just d it is of the essence to investigate a characteristic fea- n under the apex of a sandpile. This stress dip appears tureofgranularmediawithpersistingwithidenticalcon- o by localdepositionofgrainswitha pointsourcewhereas stituents. Thismakesacontrastwithpreviousstudiesin- c it does not appear by homogeneous deposition with an vestigatinghow macroscopicelastic constantsdependon [ extended source [2]. Such process dependence is often the material properties of the particles such as the stiff- 1 called as ’memory’ which sounds mysterious. A sandpile ness, the friction coefficient, and the shape [14, 15, 16]. v is merely a collection of sand grains between which the The main discovery in the present work is an existence 0 interactionisrathersimple;theyinteracteachotheronly of two time scales which are separated over two digits. 5 when they are in contact. What attracts us is that only 1 Interestingly, between the two time scales, macroscopic 4 granulesstoreinformationonhowthey arepiledup. See elastic properties are expressed as logarithmic functions . Refs. [3, 4, 5] as other related memory effects. of the strain rate. 1 0 Granularmaterialshavealsobeenofinterestasajam- 9 mingsystem[6,7]. Asanidealizedcaseofjammingtran- Model: We consider N frictional circular disks in a 0 sitions,thepackingproblemhasbeenstudiedintensively two-dimensional square box. We assume that the di- : v withfrictionlessparticles[8,9,10,11]. Theimportantre- ameter of disk i, denoted by d , is either d or d/1.4 i i sultisthatthepressureofdiskpackingsisnotzeroabove X with an equal ratio, and that its mass mi is given by the random close packing [8]. On the other hand, with m = πd2ρ/4. The interaction between disks appears r i i a frictional particles, the transition point appears between onlywhentwodisksareincontact. Supposediskiatpo- the random loose packing and the random close pack- sitionr withvelocityv andangularvelocityω contacts i i i ing [12, 13]. It might be naturally conjectured that the withdiskj atr withv andω . Usingtherelativeposi- j j j packingfractionofthetransitiondependsonthepacking tionr =r −r andthe relativevelocityv =v −v , ij i j ij i j process. However, such a dependence has not been well the normal and tangential parts of the relative velocity studied. We aim to get insight into the history depen- ofcontactingpoints,vn andvt ,aregivenasvn =v ·n ij ij ij ij dence of jamming transitions. Furthermore, we explore andvt =v ·t−(d ω +d ω )/2,respectively,with the ij ij i i j j a possibility to describe the history dependence of me- normalunit vector n=r /|r | andthe tangentialunit ij ij chanical properties of frictional disk packings. vector t. Let k and k be normal and tangential elastic n t In this Letter, for the purpose of illuminating the na- constants, η and η be normal and tangential viscous n t ture of the history dependence, we propose a frictional coefficients,andµbe the Coulombfrictioncoefficientfor disk packing in a square box without gravity as the sim- sliding friction [17]. Then, the normal force Fn and the ij plest system. We characterize the packing process only tangentialforceFt actingondiskifromdiskj aregiven ij 2 by We investigate cases with several values of δτ with δW = 5.0×10−2 fixed. We confirmed that our claim Finj = −kn(di/2+dj/2−|rij|)+ηnvinj, (1) in this Letter is independent of the value of δW if it is Ft = min( ht ,µ Fn )sign(ht ) (2) sufficiently small. We stop moving the walls at certain ij (cid:12) ij(cid:12) (cid:12) ij(cid:12) ij (cid:12) (cid:12) (cid:12) (cid:12) time t = τ and we wait until the total kinetic energy with ht = k ut +η vt . Here, ut = t vt (t)dt is the of disks becomes smaller than 10−14 so that we can as- ij t ij t ij ij Rt0 ij sumeittobeinastaticstate. Thiscompactionprocessis tangential displacement, where t is the time when the 0 characterizedbythestrainratemeasuredfromtheinitial disks come into contact. The interaction between a disk state: and the side walls of the box is calculated in the same way with the one for the interaction between disks. W(t=0)−W(t=τ) ǫ˙= . (3) The parameters in our simulations are converted into W(t=0)τ dimensionlessquantitiessothatd=1,k =1,andρ=1. n Notethatsincethefirstcontactoccursimmediatelyafter Then,theparametervaluesweusedareasfollows. d=1, we start a compaction, we assume the initial configura- µ = 0.4, k = 0.1, and η = η = 0.35. Note that the t n t tion as a reference state to determine the strain rate. normalrestitutioncoefficientisapproximatelyestimated Below we will regard the strain rate ǫ˙ as the control pa- as 0.6. We restrict our investigation to the systems with rameter of the problem we consider. N =1000. Forareferencetolaboratoryexperiments,for Evidence of the process dependence: We first study example, in the case of disks with maximum diameter of the process dependence of the final pressure with the 10mm, stiffness of 2.5MPa, thickness of 6mm, and mass 3 packing fraction fixed. In this study, we stop the shrink- density of 1.6×g/cm , the time unit in our simulations corresponds to about 3.1×10−4s. ageof the wallswhen the packingfractionreachesa pre- scribed value. The time evolution of the position, velocity and rota- For comparison, we consider the case µ = 0.0. As tionangleofdiskiisdescribedbytheequationofmotion shown in the inset of Fig. 2, the results with four differ- with using the forces described above. In numericalsim- ent strain rates yield one pressure curve as a function of ulations, we solve a set of equations of motion with em- ployingtheAdamsmethodwithatimestepof2.2×10−2. packing fraction. Here, the packing fraction of the sys- tem is measured in the domain away from the boundary (a) at the distance of 2d in order to remove boundary ef- fects. We observe that the pressure starts increasing at a critical packing fraction φ ≃ 0.84, which corresponds c to the so-calledjamming transition point [11]. It is a re- markable property that the critical packing fraction for frictionless particles is determined uniquely irrespective of packing processes. In contrast, under the existence of the tangential fric- tion with the friction coefficient of µ=0.4, the pressure curve strongly depends on the packing process as shown FIG.1: Illustration ofthepackingprocess. (a)Initialconfig- in Fig. 2. In particular, as we pack slower, the pressure uration: Disks are placed in a box. Open (green online) and curve approaches the one in the case µ = 0.0 (displayed filled (blueonline) circles represent theparticles with diame- by the triangle symbols). In order to extract a quantita- ters of 1.0 and 1.0/1.4, respectively. (b) Schematic diagram tiverelationofthe processdependence, wemeasuredthe of the motion of walls: Width of a box, W(t), is reduced as critical packing fraction φ as a function of ǫ˙. As shown a function of time. The box is shrunken isotropically with c various strain rates. in Fig. 3, our result suggests ǫ˙ As an initial state at t = 0, disks are scattered in a φc(µ=0.4,ǫ˙)−φc(µ=0)≃−αφlog (4) ǫ˙ c squareboxwithnocontact. ThewidthoftheboxW(t= 0) is determined so that the initial packing fraction is to in the regime ǫ˙ ≤ ǫ˙ ≤ ǫ˙′ with ǫ˙ ≪ ǫ˙′, where ǫ˙ and ǫ˙′ c c c c c c be 0.5. (See Fig. 1 (a).) We then carry out the following correspondto the inverseoftwo cross-overtime scalesτ c compaction procedure with moving all sides of the box and τ′, respectively. c simultaneously toward the center of the box. Let δW We shall present a phenomenological argument to un- be a small value. The width of the box is reduced with derstand Eq. (4). Let us fix p of the final state. When δW at intervals of time δτ. The time dependence of the the shrinkage is faster than the change of contact net- width of the box W is illustrated in Fig. 1 (b). At the work and the force balance is realized in this network, moment the width of the box is reduced with the strain, the packing fraction in the final state does not seriously centers of disks are shifted with an affine transformation depend on the strain rate. This provides the first cross- with an equal amount of the strain of the box. over time scale τ′. When the shrinkage is slower than c 3 this time scale, the disks have enough time to search for 0.04 0.85 ) aurmatoiroendwenitshenφemtwiogrhkt.bTehpersoepaorrcthiionngaltitmoetτhefonruamcboenrfigo-f =0.4 0.03 =0.0) 0.84 µ µ ftoorctheebacolamnbceincaotonrfiigaulrcaotmiopnlsexΣityfoorfafogricveenbaφl.anHceerset,adteuse, φ)-(c 0.02φ(c 0.83 10-5 10ε-4. 10-3 0 Σ≃exp(aN) inthe largeN limit withaconstanta that 0. = 0.01 characterizes the number of local configurations satisfy- µ ing force balance conditions. Since it is reasonable to φ(c 0 . assume that a can be expanded in φ near the packing 10-6 10-5 ε 10-4 10-3 fractionφ′ atwhichtheshrinkageratedependencestarts c increasing, Σ depends on φ in an exponential manner. The consideration leads us to FIG. 3: Deviation of the critical packing fraction φc(µ = τ ≃τc′exp(A(φ−φ′c)), (5) 0.4) from φc(µ = 0.0) as a function of the strain rate of the shrinkage, where φc(µ = 0.0) = 0.842 is assumed. Inset: where A is a constant. This tendency ceases when the Strain rate independence of critical packing fraction φc with µ=0.0. packing fraction φ is close to φ (µ=0), which provides c c thetimescaleτ asτ =τ′exp(A(φ (µ=0)−φ′))inthe c c c c c case p ≃ 0. Since the searching time τ is related to 1/ǫ˙ These elastic constants of disk packings can be mea- linearly, we obtain Eq. (4). suredbyperformingabi-axialcompressiontest[18]with the final states of the packing process. We measure in- 0.04 0.04 crements of stress, δσij, of the system when the strain, p 0.02 µ=0.0 δǫyy =0.01is givenwhileδǫx =0.0. Then,assumingthe 0.03 linear elastic theory, we estimate the values of E and ν. p 0 Note that we carefully performed the compression test 0.02 0.8φ 0.85 so that the quasi-static condition is satisfied during the 0.01 µ=0.4 compression. The results are summarized in Fig. 4. We find that 0 the largest Young’s modulus is about 1.54 times larger 0.78 0.8 0.82φ 0.84 0.86 0.88 0.9 than the smallest one. This fact clearly indicates that the elastic properties depend on the packing process de- spitethematerialpropertiesofconstituentsareidentical. Furthermore, the shapes of curves in Fig. 4 suggests FIG. 2: Pressure p versus packing fraction φ for the system with µ=0.4 and various values of strain rate ǫ˙, 1.42×10−3 E =E0−αElog(ǫ˙), (6) (Cross symbol), 3.56×10−4 (Plus symbol), 1.78×10−4 (As- ν =ν +α log(ǫ˙), (7) terisksymbol)),and4.45×10−5 (Squaresymbol)). Forcom- 0 ν parison, the result in the case µ = 0.0 and the strain rate in the regime ǫ˙ ≤ǫ˙≤ǫ˙′. of 4.45×10−5 is superposed as the triangle symbols. Inset: c c Results in the case µ = 0.0, with various strain rates. The 0.55 symbols coincide with those in the main frame, but all the 0.6 ν symbols areon onecurve. Theaverage is takenover20 sam- 0.5 ples. E 0.5 0.45 10-6 10-5 10-4 10-3 Experimental method: In many laboratory experi- ments,the packingfractionisnoteasilymeasurable. We 0.4 thus proposeamethod bywhichthe processdependence can be detected without observing the packing fraction. 10-6 10-5 1.0-4 10-3 10-2 Thekeyideahereistofocusourattentiononmechanical ε propertiesofthepacking. Concretely,wefixthepressure p in the packing. That is, we stop the compaction when the final pressure becomes an prescribed value. In the FIG. 4: Young’s modulus E and Poisson’s ratio ν (inset) as presentstudy, we consider the case p=0.025. As is seen functions of thestrain rate ǫ˙. Young’s modulus E =1 corre- from Fig. 2, the packing fraction might depend on the spondstothedimensionlessnormalstiffnessoftheconstituent strain rate of the shrinkage, and therefore it is naturally itself. 40 samples are taken for each τ. expected that Young’s modulus E and Poisson’s ration ν of the packings also depend on the strain rate ǫ˙. Here we wish to remark on a possibility of realizing 4 laboratory experiments with our idea. For example, for that θi = 2π for each i. As shown in Fig. 5, the Phjki jk the experimental system we consider in the paragraph extent of crystallization decreases rapidly between the Model, the regime 10−6 ≤ ǫ˙ ≤ 10−4 can be realized by two cross-overtime scales. Therefore,the understanding reducing the size of a box from 46cmto 38cmwith com- of the tendency to crystallization in slower operations paction time τexp satisfying 0.3(s) ≤ τexp ≤ 30(s). This might be a heart of future problems. We will study it operation is accessible in laboratory experiments. Our from several viewpoints. simple settingwillprovideusagreatadvantageto inves- tigate history dependence of granular media. We expect that the logarithmic dependence of E and ν on ǫ˙ will be The authors thank Hisao Hayakawa and Hiroki Ohta examined by laboratory experiments. for fruitful discussions. S. I. acknowledges Ken-ichi Conclusion and discussion: To conclude, we have Yoshikawa for his helpful comments and warm encour- found the history dependence of the jamming transition agement. This work was supported by grants Nos. points. For this phenomenon, the tangential friction is 18GS0421 (S. I.), 19840027 (S. I.), and 19540394 (S. S.) found to be essential. Furthermore, we have present fromtheMinistryofEducation,Science,SportsandCul- a clear evidence for the packing process dependence of ture of Japan. M. O. thanks the Yukawa Foundation for the elasticpropertiesofdiskpackingswithidenticalcon- the financial support. stituents. 0.375 0.35 ∗ Electronic address: [email protected] Ψ † Department of Physics and Mathematics, Aoyama 6 0.325 Gakuin University,Kanagawa, 229-8558, Japan [1] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev. 0.3 Mod. Phys.68, 1259 (1996). 0.275 . [2] L. Vanel, D. Howell, D. Clark, R. P. Behringer, and E. 10-6 10-5 ε 10-4 10-3 Clement, Phys.Rev. E 60, R5040 (1999). [3] C.Josserand,A.V.Tkachenko,D.M.Mueth,andH.M. Jaeger, Phys. Rev.Lett.85, 3632 (2000). [4] A.BarratandV.Loreto,Europhys.Lett.53,297(2001). [5] A. Nakahara and Y. Matsuo, Phys. Rev. E 74, 045102 FIG. 5: Bond-orientational order parameter ψ6 as afunction (2006). of the strain rate ǫ˙. [6] I.K.Ono,C.S.O’Hern,D.J.Durian,S.A.Langer,A.J. Liu,andS.R.NagelPhys.Rev.Lett.89,095703(2002). The most important question left to be answeredis to [7] O. Dauchot, G. Marty, and G. Biroli, Phys. Rev. Lett. uncover the universal nature of the logarithmic depen- 95, 265701 (2005). dence of E and ν on ǫ˙. As discussed above, the sim- [8] C. S. O’Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, pleargumentsuggeststhe logarithmicdependence ofthe Phys. Rev.Lett. 88, 105507 (2002). [9] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, packing fraction on ǫ˙ under the condition that the pres- Phys. Rev.E 68, 100306 (2003). sure p in the final state is fixed. Then, one naively con- [10] J. H. Snoeijer, T. J. H. Vulugt, M. v. Hecke, and W. v. jectures that the state with larger packing fraction with Saarloos Phys. Rev.Lett. 92, 054302 (2004). pfixedcorrespondstoastateclosertocrystalinrandom [11] A.Donev,S.Torquato,F.H.Stillinger,andR.Connelly, packings. In order to confirm this conjecture, we mea- J. Appl.Phys. 95, 989 (2004). suretheextentofcrystallizationofthepackingswiththe [12] K. Shundyak, M. v. Hecke, and W. v. Saarloos, Phys. six-fold order parameter of interbond angle defined by Rev. E75, 010301(R) (2007). [13] E. Somfai, M. v. Hecke, W. G. Ellenbroek, K. Shundya, and W.v.Saarloos, Phys.Rev.E 75,020301(R) (2007). (cid:12) 1 (cid:12) ψ6 =(cid:12)(cid:12)(cid:12)Nbond XXexp(6iθjik)(cid:12)(cid:12)(cid:12), (8) [14] FR.evR.aLdejtati.,7M7,.2J7e4an(,19J9.6J).. Moreau and S. Roux, Phys. (cid:12) i hjki (cid:12) [15] S. Luding, Phys.Rev.E. 63, 042201 (2001). (cid:12) (cid:12) [16] J.W.Landry,G.S.Grest,L.E.SilbertandS.J.Plimp- wherehjkirepresentsapairofdiskswhichareincontact ton, Phys. Rev.E 67, 041303 (2003). with disk i next to eachother,θjik is the interbondangle [17] P. A. Cundall and O. D. I. Strack, G´eotechnique 29, 47 between bond ij and ik, and Nbond is the total number (1979). of such interbond angles. Here, bond ij means the line [18] S. Inagaki and J. T. Jenkins, J. Phys. Soc. Jpn. 73, 926 segmentconnectingthecentersofdiskianddiskj. Note (2004).