Energy decay in three-dimensional freely cooling granular gas
aa r X i v : . [ c ond - m a t . s o f t ] O c t Energy decay in three-dimensional freely cooling granular gas
Sudhir N. Pathak, ∗ Zahera Jabeen, † Dibyendu Das, ‡ and R. Rajesh § The Institute of Mathematical Sciences, CIT Campus, Taramani, Chennai-600013, India Department of Physics, University of Michigan, Ann Arbor, MI 48109-1040, USA Department of Physics, Indian Institute of Technology, Bombay, Powai, Mumbai-400 076, India (Dated: November 12, 2018)The kinetic energy of a freely cooling granular gas decreases as a power law t − θ at large times t .Two theoretical conjectures exist for the exponent θ . One based on ballistic aggregation of compactspherical aggregates predicts θ = 2 d/ ( d + 2) in d dimensions. The other based on Burgers equationdescribing anisotropic, extended clusters predicts θ = d/ ≤ d ≤
4. We do extensivesimulations in three dimensions to find that while θ is as predicted by ballistic aggregation, thecluster statistics and velocity distribution differ from it. Thus, the freely cooling granular gas fitsto neither the ballistic aggregation or a Burgers equation description. PACS numbers: 45.70.Mg, 47.57.Gc, 05.40.-a, 81.05.Rm
The freely cooling granular gas, a collection of ballis-tically moving inelastic particles with no external sourceof energy, has been used to describe dynamics of granu-lar materials [1–3], large scale structure formation in theuniverse [4] and geophysical flows [5]. It is also of inter-est as a system far from equilibrium, limiting cases beingamenable to exact analysis [6, 7], has close connectionto the well studied Burgers equation [6, 8–11], and is anexample of an ordering system showing non-trivial coars-ening behavior [12–15]. Of primary interest is clusteringof particles due to inelastic collisions and the temporalevolution of the kinetic energy E ( t ) at large times.At initial times, particles remain homogeneously dis-tributed and kinetic theory predicts that E ( t ) decreasesas (1 + t/t ) − (Haff’s law) where the time scale t ∝ (1 − r ) − for constant coefficient of restitution r [16].At later times, this regime is destabilized by long wave-length fluctuations into an inhomogeneous cooling regimedominated by clustering of particles [17–19]. In this lat-ter regime, E ( t ) no longer obeys Haff’s law but decreasesas a power law t − θ , where θ depends only on dimension d [20, 21]. Direct experiments on inelastic particles underlevitation [22] or in microgravity [23, 24] confirm Haff’slaw. However, being limited by small number of particlesand short times, they do not probe the inhomogeneousregime giving no information about θ .Different theories predict different values of θ . The ex-tension of kinetic theory into the inhomogeneous coolingregime using mode coupling methods leads to E ( τ ) ∼ τ − d/ , where the relation between the average number ofcollisions per particle τ and time t is unclear [25]. This re-sult agrees with simulations for near-elastic ( r ≈
1) gases,but fails for large times and strongly inelastic ( r ≪ r = 1 is unlikely to succeed since extensivesimulations in one [20] and two [21] dimensions show thatfor any r <
1, the system is akin to a sticky gas ( r → θ mfBA = 2 d/ ( d + 2) andthe presence of a growing length scale L t ∼ t /z mfBA with z mfBA = ( d + 2) / θ BA = θ mfBA [6, 8]. However, in two dimen-sions and for dilute systems, it has been shown that θ mfBA is smaller than the numerically obtained θ BA by 17%because of strong velocity correlations between collidingaggregates [28, 29].The sticky limit has also been conjectured [20, 21] tobe describable by Burgers-like equation (BE) [30]. Thismapping is exact in one dimension [10] and heuristic intwo and higher dimensions [21], and leads to θ BE = 2 / d = 1, θ BE = d/ ≤ d ≤
4, and θ BE = 2 for d > θ mfBA and θ BE coincide with each otherin one and two dimensions and also with numerical es-timates of θ for the freely cooling granular gas in thesedimensions [20, 21]. In three dimensions, they differ with θ mfBA = 6 / θ BE = 3 /
2. However, simulations thatmeasure θ in three dimensions have been inconclusive,being limited by small system sizes and times, and themeasured value of θ ranges from θ = 1 . − . θ ∼ θ ≈ θ mfBA , conclusivelyruling out θ BE as a possible solution. Comparing withthe results of three dimensional BA, we find that θ mfBA de-scribes the energy decay in BA only when densities arehigh and multi particle collisions are dominant. We alsofind that the cluster size and the velocity distributions ofthe particles in the granular gas and BA are strikinglydifferent from each other.Consider N identical hard-sphere particles distributeduniformly within a periodic three-dimensional box of lin-ear length L and with initial velocities chosen from anormal distribution. The mass and diameter of the par-ticles are set equal to 1. All lengths, masses and timesare measured in units of particle diameter, particle mass,and initial mean collision time. The system evolves intime without any external input of energy. All particlesmove ballistically until they undergo momentum conserv-ing, deterministic collisions with other particles: if thevelocities before and after collision are u , u , and v , v respectively, then v , = u , − r n . ( u , − u , )] n , (1)where 0 < r < n is the unit vector directed from the center of particle1 to the center of particle 2. Equation (1) leaves thetangential component of the relative velocity unchanged,and reduces the magnitude of the longitudinal componentby a factor r .The above system is studied using large scale event-driven molecular dynamics simulations [36] for systemsizes up to N = 8 × . For constant coefficient ofrestitution, infinite collisions occur in finite time [37]. Anefficient scheme of avoiding this computational difficultyis to make the collisions elastic ( r = 1) when the relativevelocity is less than a cutoff velocity δ , and r = r < r = 0 .
10 and volumefraction φ = 0 . t − θ for times larger than a crossovertime that increases with system size L . We assume that E ( t ) obeys the finite size scaling form E ( t ) ≃ L − zθ f (cid:18) tL z (cid:19) , t, L → ∞ , (2)where z is the dynamical exponent, and the scaling func-tion f ( x ) ∼ x − θ for x = tL − z ≪
1. The simulationdata for different L collapse onto a single curve (seeFig. 1) when E ( t ) and t are scaled as in Eq. (2) with θ = θ mfBA = 6 / z = z mfBA = 5 /
2. The power law x − / extends over nearly 5 decades, confirming that theenergy decay in the freely cooling granular gas in three di-mensions has the exponents that are numerically indistin-guishable from the mean-field BA. The data conclusivelyrules out θ BE = 3 / f ( x ) ∼ x − η for x ≫ η ≈ . t ≫ L z , E ( t ) ∼ L . t − . .We now show that θ measured from the data in Fig. 1is independent of the volume fraction φ , coefficient ofrestitution r , and δ . The systems with varying φ -4 -2 -7 -6 -5 -4 -3 -2 -1 f ( x ) xx -6/5 x -1.83 L = 62L = 98L = 155L = 272
FIG. 1: (Color online) The data for kinetic energy E ( t ) fordifferent system sizes L collapse onto a single curve when t and E ( t ) are scaled as in Eq. (2) with θ = θ mfBA = 6 / z = z mfBA = 5 /
2. The power law fits are shown by straightlines. The data are for φ = 0 . r = 0 .
1, and δ = 10 − . are prepared by fixing L = 272 and varying N from2 × ( φ = 0 . × ( φ = 0 . φ , we find that the crossover from homogeneous( E ( t ) ∼ t − ) to inhomogeneous regime ( E ( t ) ∼ t − / )occurs at earlier times [see Fig. 2(a)]. In the inhomoge-neous regime, the curves are indistinguishable from eachother. Thus, we see that the exponent θ = 6 / φ →
0. Similarly, with increasing r , thoughthe inhomogeneous regime sets in at later times, it never-theless exists with the same power law t − θ [see Fig. 2(b)].Similar behavior has been observed in one and two di-mensions [20, 21]. We also find no discernible dependenceof the data on the parameter δ [see Fig. 2(c)]. Finally,we check that using a more realistic velocity dependentcoefficient of restitution does not change the value of theexponent θ (see supplementary material [38]).We note that θ mfBA need not be equal to the actual BAexponent θ BA [28, 29]. We study this discrepancy inthree dimensions by simulating BA directly. Two col-liding particles are replaced with a single particle whosevolume is the sum of the volumes of the colliding par-ticles. The newly formed aggregate may overlap withother particles leading to a chain of aggregation events.These multi-particle collisions result in the exponent θ BA being dependent on the volume fraction φ . We find thatas φ increases from 0 .
005 to 0 . θ BA decreases from1 . ± .
005 to 1 . ± .
005 and appears to converge tothe θ mfBA = 1 . φ . Thus, it is remarkablethat the mean field result describes well only the systemswith φ > ∼ .
2, while its derivation [27] assumes the limit φ → -8 -6 -4 -2 E ( t ) t(a) φ = 0.052 φ = 0.079 φ = 0.132 φ = 0.208 10 -8 -6 -4 -2 (b) r = 0.10r = 0.50r = 0.80r = 0.9910 -6 -4 -2 t(c) δ = 5x10 -4 δ = 1x10 -4 δ = 5x10 -5 FIG. 2: (Color online) The dependence of kinetic energy E ( t )on (a)volume fraction φ , (b) the coefficient of restitution r ,and (c) the parameter δ . The solid lines are power laws t − / .The data is for φ = 0 . r = 0 . δ = 10 − unless it is thevarying parameter. -12 -8 -4 -1 N ( m , t ) M a vg2 m/M avg t = 10000t = 20000t = 40000 10 -2 -1 -4 -3 -2 -1 N ( m , t ) M a vg2 m/M avg t = 200t = 400t = 800t = 1600 FIG. 3: (Color online) Snapshots of granular gas (upper left)and BA (upper right) in the inhomogeneous regime. Thelower panel shows the scaled mass distribution for the gran-ular gas (left) and BA (right). M avg is the mean clustersize. The solid line is a power law m − . . The data arefor φ = 0 . r = 0 . inhomogeneous regime. Snapshots of granular gas andBA (see Fig. 3) show that clusters in granular gas areextended as opposed to compact spherical clusters (byconstruction) in BA. The spatial distribution of parti-cles is partially quantified by measuring the cluster sizedistribution N ( m, t ). For the granular gas, the simula- -5 -4 -3 -2 -1 M m a x / N t granular gas, φ = 0.208granular gas, φ = 0.132BA, φ = 0.208BA, φ = 0.132BA, φ = 0.005 FIG. 4: (Color online) The largest mass M max as a functionof time. For the granular gas, r = 0 . δ = 10 − . Straightlines are power laws t . , t . , t . (bottom to top). tion box is divided into boxes of side equal to diameterof a particle. A box is said to be occupied if it containsthe center of a particle. Two occupied boxes belong tothe same cluster if connected by nearest neighbor occu-pied boxes. N ( m, t ) for the granular gas and BA, shownin the lower panel of Fig. 3, are significantly differentfrom one another. For the granular gas, N ( m, t ) con-sists of two parts: a power law ( ∼ m − . ) and a peak atlarge cluster sizes. The power law describes all clustersother than the largest cluster that accounts for the peak.The largest cluster contains about 75% of the particles.For BA, N ( m, t ) is a power law for small cluster sizes( ∼ m − . ) and exponential for cluster sizes larger thanthe mean cluster size. Both of these distributions aredifferent from the mean field result for N ( m, t ) obtainedfrom the Smoluchowski equation describing the temporalevolution of N ( m, t ):˙ N ( m, t ) = m − X m =1 N ( m , t ) N ( m − m , t ) K ( m , m − m ) − ∞ X m =1 N ( m , t ) N ( m, t ) K ( m , m ) m = 1 , , . . . , (3)where K ( m , m ) ∝ ( m − / + m − / )( m / + m / ) (4)is the collision kernel [26, 39]. For this kernel, it is knownthat that N ( m, t ) ∼ exp( − const × m − / ) for small m and N ( m, t ) ∼ exp( − const × m ) for large m [26]. Whilethe simulation results for BA matches for large m , it isdifferent (being a power law) for small m .Also, for the kernel in Eq. (4), it is expected that thelargest cluster size increases with time t as a power law t / [26, 39], the mean field answer. We compare thisprediction with the simulations for the granular gas andBA. For the granular gas, rather than a power law growthas in one and two dimensions and in mean field, there isa rapid increase in M max (see upper two curves of Fig. 4)at a time that coincides with the onset of the inhomoge-neous cooling regime. This rapid growth is similar to thegelation transition where a gel containing a fraction of thetotal number of particles is formed in finite time. How-ever the kernel for BA is non-gelling with mass dimension1 /
6, whereas the gelation transition requires mass dimen-sion to be larger than 1 [26, 39]. For BA, M max increasesas a power law (see bottom three curves of Fig. 4), withan exponent that increases with φ , and possibly convergesto the mean field value 6 /
5. Similar behavior is seen forthe growth of average cluster size of BA which growsas a power law with an exponent ranging from 1 .
06 for φ = 0 .
005 to 1 .
19 for φ = 0 . P ( v, t ),where v is any velocity component, of the granular gaswith that of BA. P ( v, t ) has the scaling form P ( v, t ) = v − rms Φ( v/v rms ), where v rms is the time dependent rootmean square velocity. The scaling function Φ( y ) is shownin Fig. 5 for different times. For the granular gas, atshort times when the system is homogeneous ( t = 5 , y ) is an exponential e − αy as predicted bykinetic theory. We find α = 2 .
65, in good agreementwith the kinetic theory value 2.60 [40]. For larger times( t = 2000 – 8000 in Fig. 5), Φ( y ) is clearly non-Gaussianand has a tail that is overpopulated compared to theGaussian (see comparison with Gaussian in Fig. 5). Aquantitative measure of the deviation from the Gaus-sian is the kurtosis, κ = h v i / h v i − /
3, shown inthe upper inset of Fig. 5. The kurtosis after an initialincrease, decreases and saturates to a non-zero value,showing quantitatively that Φ( y ) is non-Gaussian. Thelarge y behavior of Φ( y ) is shown in the bottom insetof Fig. 5. It has been argued, based on the probabil-ity that a particle never undergoes a collision up to time t , that − ln[Φ( y )] ∼ y /θ , y ≫ − ln[Φ( y )] ∼ y / , consistent with θ = 6 /
5. However, for BA we find − ln[Φ( y )] ∼ y . .This is surprising because the argument that leads to − ln[Φ( y )] ∼ y /θ [21] is quite general and does not de-pend on the detailed dynamics. Thus, in addition to hav-ing different velocity distributions, the argument basedon survival probability fails for BA.To summarize, we showed that the energy E ( t ) of athree dimensional freely cooling granular gas decreasesas t − θ , with θ ≈ /
5, indistinguishable from the meanfield result for dilute ballistic aggregation. This rules outBurgers like equations as a description of the granulargas at large times. We also showed that the relation toballistic aggregation appears coincidental with the en-ergy of the dilute ballistic gas decaying with a differentexponent. In addition, the cluster size distribution aswell as the velocity distribution of ballistic aggregationare strikingly different from that of the granular gas. We -9 -7 -5 -3 -1
0 1 2 3 4 5 6 7 Φ ( y ) y
22 exp(-2.65 y)(2 π ) -1/2 exp(-y /2)10 -2 -1 - l n Φ ( y ) y y y κ t FIG. 5: (Color online) The scaled velocity distribution func-tion Φ( y ) for the granular gas at times t = 5 ,
10 (upper col-lapsed data) and t = 2000 , , , φ = 0 . r = 0 .
10. Upper inset: The kurtosis κ as afunction of time t . Lower inset: − ln Φ( y ) as a function of y for the granular gas (lower data) and BA (upper data). ForBA, the times are t = 400 , , φ = 0 . hope that this study will prompt research into findingthe correct continuum equations for the granular gas aswell as in the design of experiments to probe the inho-mogeneous cooling regime. While frictionless freely cool-ing experiments have been limited to the homogeneousregime, inhomogeneous clustering has been observed ingranular cooling experiments where only one particle orlocation is excited [41–43]. For these systems, scaling ar-guments based on the sticky gas explain the experimentalresults [44–46]. Multiple localized excitations may resultin a crossover to the freely cooling system, making suchexperiments more suitable to probing the inhomogeneousregime.The simulations were carried out on the supercomput-ing machine Annapurna at The Institute of MathematicalSciences. ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected][1] I. S. Aranson and L. S. Tsimring, Rev. Mod. Phys. ,641 (2006).[2] T. P¨oschel and S. Luding, eds., Granular Gases (Springer, Berlin, 2001).[3] T. P¨oschel and N. V. Brilliantov, eds.,
Granular Gas Dy-namics (Springer, Berlin, 2003).[4] S. F. Shandarin and Y. B. Zeldovich, Rev. Mod. Phys. , 185 (1989).[5] C. S. Campbell, Annual Review of Fluid Mechanics ,57 (1990). [6] L. Frachebourg, Phys. Rev. Lett. , 1502 (1999).[7] S. N. Majumdar, K. Mallick, and S. Sabhapandit, Phys.Rev. E , 021109 (2009).[8] L. Frachebourg, P. A. Martin, and J. Piasecki, PhysicaA , 69 (2000).[9] R. Tribe and O. Zaboronski, Comm. Math. Phys. ,415 (2000).[10] S. Kida, J. Fluid Mech. , 337 (1979).[11] S. Dey, D. Das, and R. Rajesh, Eur. Phys. Lett. , 44001(2011).[12] S. K. Das and S. Puri, Phys. Rev. E , 011302 (2003).[13] M. Shinde, D. Das, and R. Rajesh, Phys. Rev. Lett. ,234505 (2007).[14] M. Shinde, D. Das, and R. Rajesh, Phys. Rev. E ,021303 (2009).[15] M. Shinde, D. Das, and R. Rajesh, Phys. Rev. E ,031310 (2011).[16] P. K. Haff, J. Fluid Mech. , 401 (1983).[17] I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. , 1619(1993).[18] S. McNamara and W. R. Young, Phys. Rev. E , 5089(1996).[19] E. Efrati, E. Livne, and B. Meerson, Phys. Rev. Lett. ,088001 (2005).[20] E. Ben-Naim, S. Y. Chen, G. D. Doolen, and S. Redner,Phys. Rev. Lett. , 4069 (1999).[21] X. Nie, E. Ben-Naim, and S. Chen, Phys. Rev. Lett. ,204301 (2002).[22] C. C. Maaß, N. Isert, G. Maret, and C. M. Aegerter,Phys. Rev. Lett. , 248001 (2008).[23] S. Tatsumi, Y. Murayama, H. Hayakawa, and M. Sano,J. Fluid. Mech. , 521 (2009).[24] Y. Grasselli, G. Bossis, and G. Goutallier, Euro. Phys.Lett. , 60007 (2009).[25] R. Brito and M. H. Ernst, Europhys. Lett. , 497 (1998).[26] F. Leyvraz, Physics Reports , 95 (2003).[27] G. F. Carnevale, Y. Pomeau, and W. R. Young, Phys.Rev. Lett. , 2913 (1990). [28] E. Trizac and J.-P. Hansen, Phys. Rev. Lett. , 4114(1995).[29] E. Trizac and P. L. Krapivsky, Phys. Rev. Lett. ,218302 (2003).[30] J. M. Burgers, The Non-Linear Diffusion Equation:Asymptotic Solutions and Statistical Problems (Reidel,Boston, 1974).[31] S. E. Esipov and T. J. Newman, Phys. Rev. E , 1046(1993).[32] S. E. Esipov, Phys. Rev. E , 2070 (1994).[33] S. Chen, Y. Deng, X. Nie, and Y. Tu, Phys. Lett. A ,218 (2000).[34] S. Luding, Pramana , 893 (2005).[35] S. Miller and S. Luding, Phys. Rev. E , 031305 (2004).[36] D. C. Rapaport, The art of molecular dynamics simula-tions (Cambridge University Press, Cambridge, 2004).[37] S. McNamara and W. R. Young, Phys. Fluids. A , 496(1992).[38] see Supplementary Material[39] C. Connaughton, R. Rajesh, and O. Zaboronski, in Hand-book of Nanophysics: Clusters and Fullerenes , edited byK. D. Sattler (Taylor and Francis, 2010).[40] T. P. C. van Noije and M. H. Ernst, Granular Matter ,57 (1998).[41] J. F. Boudet, J. Cassagne, and H. Kellay, Phys. Rev.Lett. , 224501 (2009).[42] X. Cheng, L. Xu, A. Patterson, H. M. Jaeger, and S. R.Nagel, Nat Phys , 234 (2008).[43] O. Johnsen, R. Toussaint, K. J. M˚aløy, and E. G.Flekkøy, Phys. Rev. E , 011301 (2006).[44] Z. Jabeen, R. Rajesh, and P. Ray, Eur. Phys. Lett. ,34001 (2010).[45] S. N. Pathak, Z. Jabeen, R. Rajesh, and P. Ray, AIPConf. Proc. , 193 (2012).[46] S. N. Pathak, Z. Jabeen, P. Ray, and R. Rajesh, Phys.Rev. E85