Chemical Physics Letters 370 (2003) 700–705 www.elsevier.com/locate/cplett
Avalanches and self-organized criticality in simulations of particle piles Christopher Rives, Daniel J. Lacks
*
Department of Chemical Engineering, Tulane University, New Orleans, LA 70118, USA Received 17 September 2002; in final form 7 January 2003
Abstract Simulations are carried out for an athermal pile of particles in the presence of a gravity-like force. As the force increases, the pile structure rearranges through discrete avalanche events. The sizes of the avalanche events are characterized by the change in energy, and it is found that the avalanche size follows a power-law probability distribution with a power-law exponent of )0.96. These results suggest that the particle pile is in a self-organized critical state. Ó 2003 Elsevier Science B.V. All rights reserved.
Jamming occurs in wide range of systems, including liquids, granular materials, and foams [1]. Applying a force can cause systems to unjam, with dynamics proceeding by discrete avalanche events. The nature of such avalanche events is not totally clear: For example, sandpile models suggest that avalanches follow power-law scaling due to self-organized criticality [2,3], while sandpile experiments can show power-law scaling [4] or a characteristic size scale [5], depending on the material and the type of experiment. The present Letter addresses avalanche events in particle-level simulations of a model system. Previous particle-level simulations have addressed the jamming transition in granular piles [6], but the avalanche-size distribution has not been investigated. The particle-level simulation of avalanches
*
Corresponding author. Fax: +1-504-865-6744. E-mail address:
[email protected] (D.J. Lacks).
in granular piles is difficult because interparticle friction is the important interaction that stabilizes the pile. To circumvent this difficulty, we address avalanches in a somewhat different system, in which attractive interparticle forces are used to stabilize the pile. This methodology allows the simulation to be carried out unambiguously, but the results are not necessarily applicable to granular piles. Simulations are carried out for a two-dimensional system with 9370 particles and a random distribution of particle sizes. A two dimensional system is used, so that finite size effects are smaller for a given number of particles (e.g., a system with a length of 100 particles along each dimension requires 106 particles in three dimensions but only 104 particles in two dimensions). A distribution of particle sizes is used to make the system disordered, since two-dimensional athermal systems with identical particles do not form fully disordered structures.
0009-2614/03/$ - see front matter Ó 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0009-2614(03)00168-4
C. Rives, D.J. Lacks / Chemical Physics Letters 370 (2003) 700–705
The particle i has the size ri , where the ri are randomly distributed between 0:8r0 and 1:2r0 . A force f acts on the particles in the direction of decreasing y, and there is a ÔfloorÕ that prevents the particles from moving below y ¼ 0. The total energy of the system is obtained as E ¼ Epforce þ Epfloor þ Epp
ð1Þ
with Epforce ¼
N X
ð2Þ
fyi ;
i¼1
Epfloor ¼
N X
4e
i¼1
Epp ¼
N N X X
ri0 yi
12 ;
V ðrij ; rij Þ;
ð3Þ
ð4Þ
i¼1 j¼iþ1
where ri0 ¼ ð1=2Þðri þ r0 Þ and rij ¼ ð1=2Þ ðri þ rj Þ. The Weber and Stillinger [7] potential is used for the particle–particle interactions, h r p r q i V ðr; rÞ ¼ eA e½1=ðar=rÞ r r ðr=r < aÞ; ð5aÞ ¼ 0 ðr=r P aÞ; ð5bÞ where p ¼ 12, q ¼ 0, a ¼ 1:652194, A ¼ 8:805997. The Weber and Stillinger potential is used because it contains a finite interaction distance while maintaining continuous derivatives at the interaction cutoff distance; since we are interested in tracking the normal mode eigenvalues with applied force, it is essential that the derivatives are continuous up to the second derivative. The following units are used throughout the Letter: e for energy, 21=6 r0 for length, and e=ð21=6 r0 Þ for the force f. The system is initially setup in a triangular pile, with the particles at the lattice positions of a hexagonal lattice with the lattice constant of 1. The force f is initially at a very small value ðf ¼ 106 Þ, because a nonzero value of f is needed to keep the pile on the floor. The particle positions are varied to minimize the energy, using a conjugate gradient method; the resulting structure is disordered due to the distribution of particle sizes (but the overall triangular shape is maintained).
701
The simulations determine how the system responds to a slowly increasing force f. Since the system is athermal and f is slowly increasing, the particles are always located at a minimum energy configuration (i.e., the time for the system to adjust to reach an energy minimum is much faster than the timescale of the change in the force f). The simulations are carried out by repeatedly increasing the force f in very small increments ðDf ¼ 105 Þ that approximate continuous change; the energy is minimized with respect to the particle positions after each increment. The results for the energy of the 9370-particle system as a function of the force f are shown in Fig. 1, and the changes in the structure of the pile with increasing f are shown in Fig. 2. These results arise from 30 813 energy minimizations at successive increments of f. As expected, the pile flattensout as the force f is increased. Due to the attractive forces between particles, the structure of the particle pile differs from the structure of a real sandpile (where the particles do not have attractive forces) – for example, the right corner of the pile is raised from the floor in Fig. 2a, and the sides of the pile are more curved than in a real sandpile. Fig. 1 shows that the energy usually increases continuously with increasing f, but these energy
Fig. 1. Energy as a function of the force f, for the 9370-particle system. The structures corresponding to points a, b, c, d and e are shown in Fig. 2. Note that the energy is on a per particle basis (i.e., the total energy divided by 9370).
702
C. Rives, D.J. Lacks / Chemical Physics Letters 370 (2003) 700–705
Fig. 2. Structure of the system at increasing force f, for the 9370-particle system. The values of the force f for the structures a, b, c, d and e are shown in Fig. 1. The dashed boxes are guides that show the dimensions of the initial structure. The percentages indicate the factors by which pictures are reduced (so as to fit on the page).
increases are punctuated by discontinuous energy drops. The discontinuous energy drops occur only after an initial period of elastic loading ðf < 0:02Þ. Similar discontinuous property changes have been observed in simulations of strained glasses [8–14]. The discontinuous energy drops correspond to avalanches in the particle pile. As an example, the energy drop and structural change associated with one of the avalanches are shown in Fig. 3. The pictures of the structure before and after the avalanche show that this particular avalanche affects
essentially the entire particle pile – i.e., this avalanche is not localized to a small number of particles. This avalanche is one of the largest that occur in the simulation (as quantified by the magnitude of the energy drop, as discussed below). The size of an avalanche can be quantified by the energy released, DE [15]. The magnitude of DE is calculated with a correction for the elastic energy involved in the force increment, " # oE DE ¼ EðfÞ Eðf DfÞ þ Df ; ð6Þ of fDf
C. Rives, D.J. Lacks / Chemical Physics Letters 370 (2003) 700–705
703
Fig. 3. Example of a discontinuous drop in energy and structural change associated with an avalanche, for the 9370-particle system. The dashed boxes are guides that show the dimensions of the structure before the avalanche. Note that the energy is on a per particle basis (i.e., the total energy divided by 9370).
where dE=df is evaluated numerically. A small cutoff value ðDEcut ¼ 0:1Þ is used to eliminate spurious energy drops that are due to noise. A total of 169 avalanches occur in the simulation of the 9370particle system (where an avalanche is defined by DE < DEcut ), and the probability distribution for avalanches releasing energy DE is shown in Fig. 4. Small avalanches have the highest probability of occurrence, and the avalanche size follows a powerlaw probability distribution with an exponent of )0.96 (i.e., probðDEÞ DE0:96 ); this power-law distribution spans four orders of magnitude. In regard to the calculation of the distribution of DE, we note two points: (1) The correction for the elastic energy (Eq. (6)) only significantly affects
the distribution for DE < 1; the distribution for DE > 1 (i.e., 10 of the 14 points in Fig. 4) is essentially unchanged if the correction is omitted. (2) The size distribution of the avalanches does not depend on the shape of the pile, which changes during the course of the simulation – as shown in Fig. 4, the results are very similar for the first and second half of the set of avalanches. As shown above, the avalanches in the particle pile follow a power-law distribution in size, with a power-law exponent near )1. Two mechanisms have been shown to give rise to a power-law distribution of avalanche events: Self-organized criticality [2] and the sweeping of an instability [16,17]. In the self-organized criticality mechanism, the
704
C. Rives, D.J. Lacks / Chemical Physics Letters 370 (2003) 700–705
Fig. 4. Probability distribution for avalanche events releasing energy DE, for the 9370-particle system. The values of DE are for the total energy (i.e., not per particle as in Figs. 1 and 3). Circles represent results for the full set of avalanches; XÕs represent results for the first half of the set of avalanches ðf < 0:081Þ; +Õs represent results for the second half of the set of avalanches ðf > 0:081Þ.
system spontaneously evolves toward a critical state characterized by fluctuations of all length scales. In contrast, the sweeping of an instability mechanism involves the slow change of a control parameter to a global instability, and the large fluctuations occur only as the global instability is approached [16,17]. In the present results, the power-law behavior including large fluctuations occurs in each of the first and second halves of the set of avalanches (see Fig. 4), which rules out the possibility that the large fluctuations are just occurring as a global instability is approached. The power-law distribution in the present results thus appears to be due to self-organized critical behavior of the particle pile. A previous study of ours investigated avalanches in a two-dimensional glass under tensile strain [18], using a very similar simulation method, the same interparticle potential, and a similar system size (9952 particles) as in the present investigation. Fig. 5 compares the present results and the strained glass results for the avalanche size distribution – the avalanche size distribution is the same for these two sets of results (to within the noise), except at the very large avalanche sizes where the strained glass results deviate from power-law behavior. This deviation from power-
Fig. 5. Probability distribution for avalanche events releasing energy DE, for the present 9370-particle system (d), compared with results for avalanches in a glassy material under tensile strain from [18] (). See text for more details of the results from [18]. Note that Figure 7 of [14] shows the results in terms of DE=E0 , where E0 ¼ 2:68 (thus the values both the x and y coordinates of the points in Figure 7 of [14] are smaller by a factor of 2.68 than in the present figure).
law behavior at very large scales for the strained glass results may be due to finite system-size effects (similar finite size effects are seen in other simulations [2]). Differences in the strained glass simulations, as compared to the present simulations, are that strain is increased rather than a force, and periodic boundary conditions are used. These differences do not seem to affect the nature of the avalanche events, except that the periodic boundary conditions apparently cause the finite systemsize effects to be more pronounced. A systematic study of system-size effects will be carried out for these two systems to further elucidate the difference in the avalanche size distribution at large avalanche sizes.
Acknowledgements Funding for this project was provided by the NSF (Grant No. DMR-0080191).
References [1] A.J. Liu, S.R. Nagel, Nature 396 (1998) 21. [2] P. Bak, C. Tang, K. Wiesenfeld, Phys. Rev. Lett. 59 (1987) 381.
C. Rives, D.J. Lacks / Chemical Physics Letters 370 (2003) 700–705 [3] For a recent review, see D.L. Turcotte, Rep. Prog. Phys. 62 (1999) 1377. [4] V. Frette, K. Christensen, A. Malthe-Sorenssen, J. Feder, T. Jossang, P. Meakin, Nature 379 (1996) 49. [5] S.R. Nagel, Rev. Mod. Phys. 64 (1992) 321. [6] L.E. Silbert, D. Ertas, G.S. Grest, T.C. Halsey, D. Levine, Phys. Rev. E 65 (2002) 051307. [7] T.A. Weber, F.H. Stillinger, J. Chem. Phys. 80 (1984) 2742. [8] D. Deng, A.S. Argon, S. Yip, Philos. Trans. R. Soc. London, Ser. A 329 (1989) 613. [9] P.H. Mott, A.S. Argon, U.W. Suter, Philos. Mag. A 67 (1993) 931. [10] C.F. Fan, Macromolecules 28 (1995) 5215.
705
[11] D.J. Lacks, Phys. Rev. Lett. 80 (1998) 5385. [12] D.L. Malandro, D.J. Lacks, Phys. Rev. Lett. 81 (1998) 5576. [13] D.J. Lacks, Phys. Rev. Lett. 87 (2001) 225502. [14] M. Utz, P.G. Debenedetti, F.H. Stillinger, Phys. Rev. Lett. 84 (2000) 1471. [15] S. Tewari, D. Schiemann, D.J. Durian, C.M. Knobler, S.A. Langer, A.J. Liu, Phys. Rev. E 60 (1999) 4385. [16] D. Sornette, J. Phys. I (France) 4 (1994) 209. [17] D. Sornette, Phys. Rev. Lett. 72 (1994) 2306. [18] G. Gagnon, J. Patton, D.J. Lacks, Phys. Rev. E 64 (2001) 051508.