Role of defects in the onset of wall-induced granular convection
RRole of defects in the onset of wall-induced granular convection
Andrea Fortini , ∗ and Kai Huang Theoretische Physik II, Physikalisches Institut, Universit¨at Bayreuth,Universit¨atsstraße 30, D-95447 Bayreuth, Germany Department of Physics, University of Surrey, Guildford GU2 7XH, United Kingdom and Experimentalphysik V, Physikalisches Institut, Universit¨at Bayreuth,Universit¨atsstraße 30, D-95447 Bayreuth, Germany
We investigate the onset of wall-induced convection in vertically vibrated granular matter bymeans of experiments and two-dimensional computer simulations. In both simulations and ex-periments we find that the wall-induced convection occurs inside the bouncing bed region of theparameter space in which the granular bed behaves like a bouncing ball. A good agreement betweenexperiments and simulations is found for the peak vibration acceleration at which convection starts.By comparing the results of simulations initialised with and without defects, we find that the onsetof convection occurs at lower vibration strengths in the presence of defects. Furthermore, we findthat the convection of granular particles initialised in a perfect hexagonal lattice is related to thenucleation of defects and the process is described by an Arrhenius law.
PACS numbers: 45.70.-n, 61.72.Bb, 47.55.P-
I. INTRODUCTION
Mixing and demixing of vibrated granular matter [1, 2]are of importance in nature [3] as well as in many indus-trial processes. For example, they are used in the phar-maceutical, construction [4] or waste reprocessing [5] in-dustries. The term Brazil Nut Effect (BNE) [6], whichoriginally referred to the rise of a large particle to thetop of a container filled with smaller grains, is now usedto indicate the more general demixing of differently sizedparticles under vertical oscillations. Schr¨oter et al. [7]reviewed and identified seven possible mechanisms thatlead to the BNE and to the Reverse Brazil Nut Effect(RBNE) [8]. Convective cells induced by the walls of thecontainer were found to be a major contributing mecha-nism to the occurrence of the BNE [9–13].The complexity of the BNE characterisation is in partdue to an underlying dynamical behaviour which evenfor one component systems is very rich. At driving accel-erations smaller than the gravitational acceleration, thegranular bed comoves with the bottom wall. Upon in-creasing the acceleration the granular bed behaves like abouncing ball [14], and above this bouncing bed regioncollective undulations (also known as arches) appear [15–18]. At still higher accelerations the granular
Leidenfrost effect occurs, in which a dense granular fluid hovers overa granular gas. Recently, Eshuis et al. [19] systematicallydrew phase diagrams for all these phenomena, and foundat very high accelerations a convective regime in whichthe sample is completely fluidized [18, 19].However, this convection regime is distinct from thewall-induced convection which occurs at low accelera-tions and is one of the driving mechanisms for the BNE.The wall-induced convection is caused by the shear forces ∗ [email protected] between particles and walls. During the upward acceler-ation the mixture gets compacted and shear forces in-duced by the side walls propagate efficiently through thewhole sample. During the downward motion the mixtureis more expanded and consequently those particles ad-jacent to the walls experience stronger downward shearforces than those in the centre of the container. Thecombination of the two type of motions gives rise to con-vection [7]. This cycle of expansion and compression ofthe granular bed was studied by Sun et al. [20] and foundto be strongly dependent on wall friction. Even thoughtheir study was done in relation to BNE, no connection tothe convective motion is made. The onset of convectionhas been studied extensively in both experiments [21, 22]and with numerical simulations [23–26].In this article, we study with both computer simu-lations and experiments the dynamical phase diagramof two-dimensional vertically oscillated granular matterand analyse the mechanism behind the onset of the wall-induced convection to clarify the role of topological de-fects. II. METHODS
For the theoretical investigation we carry out Molecu-lar Dynamic(MD) [27] simulations of a two-dimensionalsystem in a box of size L x × L z delimited by flat hardwalls and gravity g = − g e z pointing in the negative z direction. The particles have two translational degrees offreedom in the x − and z − directions and one rotationaldegree of freedom about the perpendicular y -axes. Thegranular beads are described as soft disks of diameter σ ,mass m , and moment of inertia I = mσ /
8, which inter-act via a linear contact model with viscoelastic dampingbetween the disks and via static friction [28]. This modeland its parameter values (See Tab. I) have been cho-sen because they reproduce the contact properties [29] of a r X i v : . [ c ond - m a t . s o f t ] M a r granular beads, and give results for the dynamics of anintruder in a vertically oscillated granular bed, that arein good agreement with experiments [20].The simulation box is driven sinusoidally, i.e., the bot-tom of the container is moved in time according to z b = − L z + A sin ( ωt ) , (1)where z b is the height coordinate of the bottom of thecontainer, A is the amplitude of the oscillation, ω is thefrequency and t is the time.We traced the dynamical phase diagram for a fixedoscillation frequency ω = 1 . t − and lateral wall sepa-ration L x /σ = 20. Reduced units are used throughoutthe article: the particle mass m , the particle diameter σ and the gravitational acceleration g are our fundamen-tal units. Consequently, the derived units are the time t = (cid:112) σ/g , velocity v = √ gσ , force f = mg , elasticconstant k = mg/σ and damping coefficient γ = (cid:112) g/σ .Further details of the model are given in appendix A.The external driving force is characterised by a di-mensionless acceleration Γ = Aω /g , corresponding tothe maximum acceleration due to Eq. (1) divided by thegravitational acceleration g . Alternatively, we use the di-mensionless energy parameter [30] K m = A ω σg , i.e., themaximum kinetic energy per particle ’injected’ in the sys-tem every period of oscillation [31].The experiment is conducted with a monolayer ofspherical polished opaque glass beads (SiLiBeads P) witha diameter of 2 ± .
02 mm. A rectangular cell made upof two glass plates 40 mm width by 200 mm height sep-arated by a distance of 2 . N = 200 , , , , N l = N/L x . This number gives an approximate value for thenumber of particle layers and hence, the height of thegranular bed. In reality the height depends on the localstructure of the granular bed, namely, the orientation of the hexagonally-ordered particles, and the type andamount of defects. III. INITIALISATION
In the experiment the system is initialised with astrong agitation to create a completely fluidized state,followed by a slow ramping down of the vertical accel-eration. In the simulation the initial configuration isprepared with two distinct procedures. In the crystalinitialisation, the particles are placed in a perfect hexag-onal lattice resting at the bottom of the container. In therandom initialisation, the particles are placed randomlyin the box. Via molecular dynamics we evolve the sys-tem until all particles have fallen under gravity and havereached a rest position.
Figure 1. (colour online) Simulation snapshots of the typicalinitial configurations for N l = 60. The colour indicate thelocal order of the particles according to an analysis with the q order parameter [33]. Green (light grey) indicates hexago-nal order, red (dark grey) indicates a square local order whileother colours indicate a disordered configuration. From leftto right: the crystalline and random initialisations of the sim-ulation, and the experimental initial configuration. The resulting configurations are shown in Fig. 1. Wenote that in both simulations initialised randomly andthe experiments a large amount of particles have localhexagonal order with many topological defects, such asdislocations and grain boundaries. For a two dimensionalsystem there are two favoured orientations of the hexag-onal lattice, when in contact with a flat wall. In Fig. 2a)the Orientation A with the [111] direction parallel to thewall is shown. In Fig. 2b) the orientation B has the [010]direction parallel to the wall. In the crystalline initialisa-tion we chose to place the particles according to orienta-tion A. After the random initialisation processes, we findthat in the experiment almost all particles have orienta-
Figure 2. (colour online) Sketch of the possible orientationsof the hexagonal crystal in contact with the lateral wall of thecontainer box. tion B, while in simulation the orientation A seems to befavoured, but grain boundaries between the different ori-entations are visible. Figure 1 shows the configurationsobtained in the three cases.
IV. DYNAMICAL PHASES
We have systematically investigated the onset of con-vection at different dimensionless accelerations/energiesand heights of the granular bed. Figure 3 shows the lo-cation of the dynamical phases as a function of the di-mensionless acceleration Γ and the energy parameter K m for different values of the linear density N l . The param-eter N l σ is also a measure of the number of layers in thegranular bed, and therefore of its height.The transition from the comoving to bouncing bedphase occurs, as expected, at Γ slightly larger thanone, with the critical value of the acceleration increasingslightly with increasing number of layers. We find goodagreement between simulations (green triangles/blue cir-cles) and experiments (connected squares). In experi-ments, the determination of this boundary is sensitiveto the accuracy of the particle position determination,which is set by the resolution of the camera.Inside the bouncing bed regime, we observe a transi-tion from a bouncing bed phase without convection to abouncing bed phase with convection. We find reasonableagreement between experiments and the simulations ini-tialised randomly. For N l σ =60,80 the model underesti-mates the amount of energy necessary for the convectionto start. The discrepancy is probably related to vari-ability in the concentration of defects, as well as to thepresence of the front and back wall in the experiments.For simulations that start without topological defects weconsistently detect the onset of convection (blue circles)at higher values of the acceleration with respect to theexperiments and to the simulations started with defects.The enhancement of the convection due to the presence of defects explains the observation of P¨oschel and Her-rmann [11] that an intruder can initiate convection. Thepresence of an intruder induces topological defects [34]that initiate the convection and segregation. N l σ Γ K m Co-movingBouncing bed: No ConvectionBouncing bed: With ConvectionSimulation: CrystalSimulation: RandomExperiment
Figure 3. (colour online) Dynamical phases at different accel-erations Γ and linear density/number of layers N l σ . On theright y -axis the energy parameter K m is also reported. Sym-bols are for simulation results for either the system initialisedrandomly (triangles) or in a crystalline configuration (circles).The connected squares are the experimental results. V. ONSET OF CONVECTION
In crystalline materials collective particle movementsoccur via crystalline plane slips [10]. The same occurs inour granular systems, but we noted some differences inslip behaviour between simulation and experiments. Insimulation the slip occur primarily along the oblique di-rections, while in the experiment it occurs mainly alongthe vertical direction. The difference can be explainedby the different orientation of the hexagonal crystal inthe two cases. The granular particles experience a shearstress due to the wall of the container along the verticaldirection. In the experiments the particles are orientedlike in Fig. 2b), and slips occur preferentially in the ver-tical direction. On the other hand, in simulations we findmore particles with an orientation like in Fig. 2a) and thepyramidal slip planes will be activated first. As long asplane slips are concerned, the worst case scenario is rep-resented by a perfect crystal with orientation A. Sinceno defects are initially present, slips along the pyrami-dal planes can only occur if defects are nucleated first.Another difference is that in experiments we often ob-serve single convective rolls similar to the observation ina two-dimensional rotating cell [35]. This type of convec-tive motion is not detected in simulation and the reasonfor the discrepancy is likely due to a small tilt of the sidewalls [9].In order to clarify the role of defects in the onset of con-vection, we analyse in computer simulations the granulartemperature T g in relation to the average input energy K m The granular temperature is defined as T g = mN (cid:104) N (cid:88) i =1
12 ( v i ( t ) − v cm ( t )) (cid:105) (2)where v cm is the centre of mass velocity and v i ( t ) is thevelocity of particle i at time t and (cid:104)(cid:105) indicates a timeaverage. Figures 4a-b show the behaviour of the reducedgranular temperature T ∗ g = T g / ( mgσ ) as a function ofthe dimensionless energy parameter K m for systems ini-tialised with and without defects, respectively. In bothcases the temperature increases monotonically over manyorders of magnitude. The curve gives an indication ofhow much input energy K m is converted into the kineticenergy of the granular particles.We divide the energy K m parameter space in regionswith different slopes of the temperature curves. For therandom initialisation (Fig. 4a) we find two regions: A and B . The transition occurs at K m (cid:39) . B the region where convection is detected. The region A , without convection, is now divided in two subregions A and A .In order to clarify the origin of the different regionsof temperature behaviour, we investigate the slip prob-ability of a particle in the bulk of the granular bed, i.e.away from the top free interface [36]. The slip of a par-ticle is detected if a displacement larger than 0 . σ ismeasured for at least two nearest neighbours after oneoscillation. We define the slip probability p s as the frac-tion of particles which undergo a slip in one oscillationperiod. The value of p s is averaged over 200 periods ofoscillation. The slip probability as a function of 1 /K m is shown in Figs. 5a-b for the random and crystal ini-tialisation, respectively. For the random case we do notobserve a change in the slip behaviour between regions A and B . On the other hand a quite dramatic change ofbehaviour is observed for the crystal case. In particular,in region A the number of detected slip is exactly zero,while region A is characterised by a quite steep increaseof slip events, which slows down considerably upon onsetof convection (region B ).In all cases the curves can be fitted to an Arrheniuslaw, p s ∝ exp( − E b /K m ), indicating the presence of acti-vating mechanisms. The energy parameter E b is a mea-sure of the barrier height and the dimensionless energy K m assumes the role of a reservoir temperature. For therandom case just a single curve can be fitted to the entirerange of inverse energies. From the fit we find an averagedimensionless barrier height E b = 7 ±
1. For the case ini-tialised with the crystal, we find a barrier E b = 130 ± A , and a barrier E b = 7 ± B .Since the barrier height is the same for the system withinitial defects, independent of the presence of convection,we can speculate that the activation mechanism indicatedby the Arrhenius equation is the activation of slip eventsin a crystal with topological defects. Interestingly, thesame barrier is measured for the system without initialdefects in region B, suggesting the same activation mech-anism. But in order for this to occur defects must firstnucleate inside a perfect crystal, and we speculate thatthe very high barrier in region A is due to the nucleationof defects in a perfect hexagonal crystal. K m -3 -2 -1 T ∗ g A B a) Random N l =40 N l =60 K m -4 -3 -2 -1 T ∗ g A BA b) Crystal Figure 4. (colour online) Granular temperature T g versus di-mensionless energy parameter K m . The dashed lines separateregions with different slopes of the temperature curves. a) Forthe system initialised randomly, i.e. with defects. b) For thesystem initialised in a perfect hexagonal lattice, i.e. withoutinitial defects. VI. CONCLUSIONS
In conclusion, using two-dimensional computer simula-tions and experiments we locate the wall-induced convec-tion in the bouncing bed region in the dynamical phasediagram of vertically vibrated granular matter.For the onset of convection, we find a reasonable agree-ment between the experimental results and those fromsimulations initialised with defects. We believe the rea-son behind the discrepancy is the presence of a front anda back wall in the experimental box, as well as some vari- /K m -5 -4 -3 -2 -1 p s a) Random B A /K m -5 -4 -3 -2 -1 p s b) Crystal A B A N l =40 N l =80 Figure 5. (colour online) Slip probability p s versus the inversedimensionless energy parameter 1 /K m . The continuos linesrepresent exponential fit functions. The dashed lines separatethe regions defined in Fig. 4. a) For the system initialisedrandomly, i.e. with defects. b) For the system initialised in aperfect hexagonal lattice, i.e. without initial defects. ability, due to the presence of defects. Other phenomeno-logical differences between simulations and experimentsare observed for the favoured crystal orientation and slipevents. We think that these differences are due to a smalltilt of the experimental box and small variations of theexperimental box width L x .For a system initialised in a perfect crystal, i.e., with-out initial defects, we consistently observe the wall-induced convective motion to occur at higher values ofthe dimensionless energy parameter (shaking strength)with respect to the system initialised randomly. Fromthe analysis of the granular temperature we distinguishdifferent regions based on the slope of the temperature asa function of the dimensionless energy K m . For systemsinitialised randomly, i.e. with many initial defects, wedistinguish two regions ( A and B ) of temperature. Thetransition between the two regions occurs at the onset ofconvection, but we do not observe any change in the slipprobability in the transition between regions A and B .On the other hand, for systems initialised with a per-fect hexagonal lattice, i.e. without initial defects, we dis-tinguish three regions in the temperature curve. More-over, in this case the slip probability has very differentbehaviour in the three regions. In particular, in region A the number of detected slips is exactly zero, while re-gion A is characterised by a quite steep increase of slipevents with the input energy, which slows down consid-erably upon onset of convection (region B ).The slip probability follows an Arrhenius law, indicat-ing the presence of an activating mechanisms. The very high barrier in the region A of the system without de-fects is due to the nucleation of defects in the perfecthexagonal crystal. We noted that the barrier in region B is the same for the system initialised with and withoutinitial defects. In this region the behaviour of the slipprobability is related to the activation of slip events.The results of our work provide an explanation forthe onset of the convective regime in vertically oscillatedgranular systems and show that defects enhance onsetof convection, i.e. systems with topological defects showconvection at lower oscillation strengths, with respect tosystems without defects. More work is needed to quan-tify the degree of variability for the onset of convectiondue to its sensitivity to the concentration and possiblytypes of defects. Interestingly, since defects can diffuseout of the system when they reach the top of the granularbed, regions of transient convective motion are possible,provided that the time scale for the defects diffusion islarger than the time scale for the nucleation of defects.This conjecture was not studied in this work, but repre-sents an interesting avenue of future research. Further-more, we plan to study how the size polydispersity of thegranular particles changes the onset of the wall-inducedconvection as well as the influence of the box size onthe dynamical behaviour. Furthermore, we would like toexplore the effect of shock waves [37] on the dynamicalbehaviour of the system. Appendix A: Details of the Model
We carry out Molecular Dynamics simulations at fixedtime step dt , for two translational degrees of freedom andone rotational degree of freedom for a system of soft disks.An illustration of the model is shown in Fig. 6. Twoparticles at positions r i and r j with velocities v i and v j and angular velocities ω i and ω j define a system with aneffective mass m eff = m i m j / ( m i + m j ) and a normal unitvector n ij = r i − r j | r i − r j | . The tangential direction is definedas t ij = v tij | v tij | where v t ij = v ij − v n ij −
12 ( σ i ω i + σ j ω j ) × n ij , (A1)with v n ij = ( v ij · n ij ) n ij .We define the displacement in the two directions δ n ij = d − | r i − r j | , with d = 1 / σ i + σ j ), and δ t ij t ij = v t ij dt .For the disk-disk interactions we use a linear model withforces F n ij = ( κ n δ n ij n ij − γ n m eff v n ij ) (A2) F t ij = ( − κ t δ t ij t ij − γ t m eff v t ij ) , (A3)in the normal and shear tangential directions, respec-tively. The parameters κ n and κ t are the stiffness co-efficients in the normal and tangential direction, respec-tively. The energy dissipated during the contact is reg-ulated by the damping coefficients γ n and γ t . In ad-dition we model the static friction by keeping track ofthe elastic shear displacement δ t ij over the contact life-time and truncate it such that the Coulomb condition | F t ij | < | µF n ij | is satisfied, where µ is the static frictioncoefficient. The same kind of interaction is used between Figure 6. Illustration of the contact model of two granulardisks at distance r ij = r i − r j and relative velocity v ij = v i − v j , which overlap by a distance δ n ij . the particles and the container wall. We also considerthe gravitational force g = − g e z , where e z is the unitvector pointing in z direction.Once the forces on all particles are known the totalforce and torque τ i on a particle i is determined by F i = m g + (cid:88) j ( F n ij + F t ij ) (A4) τ i = − (cid:88) j σ j n ij × F t ij . Coefficient Particle-Particle Wall-ParticleNormal stiffness k n ( k ) 10 ( k )Tangential stiffness k t ( k ) 10 ( k )Static friction µ γ n
100 ( γ ) 100 ( γ )Tangential damping γ t
100 ( γ ) 100 ( γ )Table I. Numerical values of the simulation parameters. The typical numerical values of simulation parametersare shown in Table I. In this linear model it is possibleto calculate the contact duration [29] t c = π (cid:18) k n m eff − γ n (cid:19) − . . (A5)In order to obtain an accurate integration of the equationof motion during contact, the time step of the simulationis chosen to be δt ≈ t c /
50 [38].
Appendix B: Detection of the dynamical phases
The detection of the convective motion, in both com-puter simulations and experiments, is performed by using N t tracer particles, which are initially positioned at thecentre of the oscillating box in a straight horizontal line.In order to determine the threshold for convection, weanalyse the deviation of the imaginary line connectingthe tracer particles from the initial straight horizontalconfiguration. This is carried out by calculating, at thebeginning of each cycle, the variance the tracers heightfrom their average height s j = 1 N t N t (cid:88) i =1 (cid:0) z ij − (cid:104) z j (cid:105) (cid:1) , (B1)where (cid:104) z j (cid:105) is the average vertical position of the tracersat frame j and z ij is the vertical position of tracer i atthe beginning of cycle j .An average variance over all N f frames is calculated s = 1 N f N f (cid:88) j =1 s j , (B2)and the start of the convection is identified via the con-dition s c > σ , where σ is the diameter of the grains.The detection of the bouncing bed dynamical phase isperformed by detecting the detachment of the granularbed from the bottom plate, which can occur at any phaseof the oscillating cycle. Therefore, we calculate the av-erage height of the tracers with respect to the oscillatingplate for each phase ks ( k ) = 1 N f N f (cid:88) j =1 ( (cid:104) z j ( k ) (cid:105) i − z b ( k )) , (B3)where (cid:104) z j ( k ) (cid:105) i = N t (cid:80) N t i =1 z ij ( k ) is the average height ofall N t tracer particles at phase k and cycle j , z b ( k ) isthe height of the plate at phase k , and N f is the totalnumber of frames at a certain phase.Consequently, the variance over all N k phases is calcu-lated s = 1 N k N k (cid:88) k =1 s ( k ) , (B4)and the bouncing bed phase is identified with the con-dition s b > σ/
30. The value σ/
30 corresponds to half apixel of the experimental images. For the sake of com-parison, the same value is used in the analysis of thesimulation trajectories.
ACKNOWLEDGMENTS
The authors thank Ingo Rehberg and MatthiasSchmidt for discussions and acknowledge Philipp Ram-ming for the help in image processing. Maximilian vonTeuffenbach, Andreas Fischer and Philip Krinninger areacknowledged for helping with the initial development of the granular simulation code as a part of their final bach-elor projects. K.H. is supported by the DFG throughGrant No. HU1939/2-1. [1] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev.Mod. Phys. , 1259 (1996).[2] S. Luding, Nonlinearity , R101 (2009).[3] H. Miyamoto, H. Yano, D. J. Scheeres, S. Abe,O. Barnouin-Jha, A. F. Cheng, H. Demura, R. W.Gaskell, N. Hirata, M. Ishiguro, T. Michikami, A. M.Nakamura, R. Nakamura, J. Saito, and S. Sasaki, Sci-ence , 1011 (2007).[4] J. Duran., Sands, powders, and grains: An introductionto the physics of granular materials. (Springer Verlag,2000).[5] N. Mohabuth and N. Miles, Resources, Conservation andRecycling , 60 (2005).[6] A. Rosato, K. Strandburg, F. Prinz, and R. Swendsen,Phys. Rev. Lett. , 1038 (1987).[7] M. Schr¨oter, S. Ulrich, J. Kreft, J. Swift, and H. Swinney,Phys. Rev. E , 011307 (2006).[8] D. Hong, P. Quinn, and S. Luding, Phys. Rev. Lett. ,3423 (2001).[9] J. Knight, H. Jaeger, and S. Nagel, Phys. Rev. Lett. ,3728 (1993).[10] W. Cooke, S. Warr, J. Huntley, and R. Ball, Phys. Rev.E , 2812 (1996).[11] T. P¨oschel and H. J. Herrmann, Europhysics Letters ,123 (1995).[12] A. Kudrolli, Rep. Prog. Phys. , 209 (2004).[13] M. Majid and P. Walzel, Powder Technology , 311(2009).[14] A. Mehta and J. Luck, Phys. Rev. Lett. , 393 (1990).[15] S. Douady, S. Fauve, and C. Laroche, Europhysics Let-ters , 621 (1989).[16] A. Ugawa and O. Sano, Journal of the Physical Societyof Japan , 1390 (2003).[17] O. Sano, Phys. Rev. E , 051302 (2005).[18] P. Eshuis, R. Bos, D. Lohse, D. Van der Meer, andK. Van der Weele, Phys. of Fluid , 123301 (2007).[19] P. Eshuis, D. van der Meer, M. Alam, H. J. van Gerner,K. van der Weele, and D. Lohse, Phys. Rev. Lett. ,038001 (2010). [20] J. Sun, F. Battaglia, and S. Subramaniam, Phys. Rev.E , 061307 (2006).[21] E. Clem´ent, J. Duran, and J. Rajchenbach, Phys. Rev.Lett. , 1189 (1992).[22] J. B. Knight, E. E. Ehrichs, V. Y. Kuperman, J. K. Flint,H. M. Jaeger, and S. R. Nagel, Phys. Rev. E , 5726(1996).[23] Y.-h. Taguchi, Phys. Rev. Lett. , 1367 (1992).[24] S. Luding, E. Cl´ement, A. Blumen, J. Rajchenbach, andJ. Duran, Phys. Rev. E , R1762 (1994).[25] M. Bourzutschky and J. Miller, Phys. Rev. Lett. , 2216(1995).[26] D. Risso, R. Soto, S. Godoy, and P. Cordero, Phys. Rev.E , 011305 (2005).[27] D. Frenkel and B. Smit, Understanding Molecular Simu-lation (Academic Press, London, 2002).[28] P. A. Cundall and O. D. L. Strack, G´eotechnique , 47(1979).[29] J. Sch¨afer, S. Dippel, and D. Wolf, J. Phys. I France ,5 (1996).[30] H. Pak and R. Behringer, Physical Review Letters ,1832 (1993).[31] This parameter is also called dimensionless shakingstrength [18].[32] C. Kimme, D. H. Ballard, and J. Sklansky, Comm. As-soc. Comp. Mach. , 120 (1975).[33] P. Steinhardt, D. Nelson, and M. Ronchetti, Phys. Rev.B , 784 (1983).[34] V. W. de Villeneuve, R. P. Dullens, D. G. Aarts,E. Groeneveld, J. H. Scherff, W. K. Kegel, and H. N.Lekkerkerker, Science , 1231 (2005).[35] F. Rietz and R. Stannarius, Phys. Rev. Lett. , 118001(2012).[36] We have chosen to exclude particles closer than 10 σ tothe free interface.[37] K. Huang, G. Miao, P. Zhang, Y. Yun, and R. Wei,Phys. Rev. E , 041302 (2006).[38] J. Lee, J. Phys. A27