Density profiles of loose and collapsed cohesive granular structures generated by ballistic deposition Dirk Kadau and Hans J. Herrmann IfB, HIF E12, ETH Ho¨nggerberg, 8093 Zu¨rich, Switzerland (Dated: January 7, 2011) Loosegranularstructuresstabilizedagainst gravitybyaneffectivecohesiveforceareinvestigated onamicroscopicbasisusingcontactdynamics. WestudytheinfluenceofthegranularBondnumber onthedensityprofilesandthegenerationprocessofpackings,generatedbyballisticdepositionunder gravity. Theinternalcompactionoccursdiscontinuouslyinsmallavalanchesandwestudytheirsize distribution. We also develop a model explaining the final density profiles based on insight about 1 thecollapse of a packingunderchanges of theBond number. 1 0 PACSnumbers: 47.57.-s,45.70.Mg,83.80.Hj 2 n a I. INTRODUCTION the density profiles depend on the granular Bond num- J ber, i.e. the ratio of cohesive force by gravity, and what 6 theinfluenceofthedynamicsofdeposition/collapse. As Loose granular packings, metastable granular struc- discussedaboveloosestructuresaregeneratedinnature, ] tures and fragile granular networks play an important industrial application, experiments or simulation by dif- t f role in a wide range of scientific disciplines, such as ferent processes. Here, we focus on ballistic deposition. o collapsing soils [1–4], fine powders [5] or complex flu- However, we expect the findings of this paper to be of s ids [6, 7]. In collapsing soils without any doubt there . relevance to all systems involving compactiondue to the t is a metastable or fragile granular network involved [1– a particles’ own weight. m 3, 8, 9]. A similar failure behavior can be found in col- Afteradescriptionofthesimulationmodelandabrief loidal gels [10] and snow [11, 12]. But also powders have - discussion of possible experimental realizations in sec- d inmostcasesaneffectivecohesiveforce,e.g.duetoacap- tion II, we first study the resulting density profiles when n illarybridgebetweentheparticlesorvanderWaalsforces gravity acts during deposition, in particular the influ- o (important when going to very small grains, e.g. nano- c ence of the granular Bond number (sec. III). To under- particles) leading to the formation of loose and fragile [ stand the shape of the density profiles we study in the granularpackings [5, 8, 9, 13]. In many complex fluids a following (sec. IV) the role of the dynamics of the col- 2 fragile/metastable network of colloids/grains is believed lapse occurring in small avalanches. We study the aver- v to be the essential ingredient for the occurrence of shear 9 age “avalanche profile” defined here as the average dis- thickening [6] or yield stress behavior [7]. 3 tance a particle moves downwards after being deposited 3 Thegeneralfeatureofsuchfragilenetworksisthatthey depending on its height. We observe characteristic pro- 5 cancollapse/compactunder the effect of anapplied load files which can be used to relate the final density pro- . 0 [9, 11, 13]. This load can be an external load or exerted file to the deposition density, given by the number of 1 internally by a force acting on all particles within the deposited particles per unit volume (sec. V). To under- 0 structure. This “internal collapse” is important in dif- stand this phenomenological profile we study a simpler 1 ferent applications like cake formation of filter deposits system where first all particles are deposited followedby : v [14–16],wherethe compactionforceinmostsituations is the collapse of the whole structure leading to an even i the drag force exerted on the grains by the flow which simpler profile (sec.VI). Our calculationsyield that this X istypicallyporositydependent [17]. Thestructure’sown linear profile is obtained in all processes where a homo- ar weight leads to compaction of snow after deposition [18] geneousinitialconfigurationiscollapsed/compactedtoa and during aging [19, 20], or to sediment compaction homogeneous final state. In Sec.VII we show that the [21–23]. In all cases, typically a depth dependent poros- phenomenologicalobtained avalancheprofile obtained in ity is observedandquantified by continuumdescriptions sec. III can be derived from the linear avalanche profiles [15,16,18,21–23]. Inmostcasesthedetailsoftheporos- of the homogeneous collapse. ityprofileareinfluencedbyacombinationofdifferentme- chanicalandchemicalprocesses[15,21]. Itiswellknown that the porosity of a structure is of major importance II. DESCRIPTION OF SIMULATION MODEL for its mechanicalproperties [11, 21, 24, 25], in filtration processes [26] and its chemical properties like catalytic The dynamical behavior of the system during gener- activity[27]. Theaimofthispaperistostudythemicro- ation is modeled with a particle based method. Here scopicprocesses,i.e.onthegrainscale,fortheseinternal we use a two dimensional variant of contact dynamics, compaction processes. For this, we will investigate the originally developed to model compact and dry systems compactionduetogravityinasimplifiedmodelsystemof with lasting contacts [28–31]. The absence of cohesion grains held together by cohesive bonds. We analyze how between particles can only be justified in dry systems 2 on scales where the cohesive force is weak compared to fect of gravity as described in the previous section. It is the gravitational force on the particle, i.e. for dry sand expected that the density and the characteristics of the and coarser materials, which can lead to densities close densityprofilesaremainlydeterminedbytheratioofthe to that of random dense packings. However, an attrac- cohesive force F to gravity F , typically defined as the c g tive force plays an important role in the stabilization of granular Bond number Bo =F /F [46, 47]. Obviously g c g large voids [32], leading to highly porous systems as e.g. the case of Bo =0 corresponds to the cohesionless case g infinecohesivepowders,inparticularwhengoingtovery whereas for Bo → ∞ gravity is negligible. A similar g small grain diameters. Also for contact dynamics a few dimensionless quantity had been identified as most im- simple models for cohesive particles are established [32– portant parameter in previous studies on compaction of 35]. Here we consider the bonding between two particles cohesive powders [32, 36, 48]. in terms of a cohesion model with a constant attractive In the following, we use monodisperse systems with a force Fc acting within a finite range dc, so that for the friction coefficient µ = 0.3 and a rolling friction coeffi- opening of a contact a finite energy barrier Fcdc must cientofµr =0.1(inunits ofparticleradii). Theeffectof be overcome. In addition, we implement Coulomb and varying these parameters is also studied exemplary and rolling friction between two particles in contact, so that will be discussed later. Typically the values of the den- large pores can be stable [32, 36–39]. sity can depend on these parameters as shown in Ref. To generate the loose structure we use ballistic depo- [37] whereas the qualitative behavior does not change. sition where each deposited particle, chosen at random Figure 1 shows the final structures obtained for differ- horizontalposition, is attachedto the structure atmaxi- ent values of granular Bond number ranging from 0 to mal possible height with zero velocity. At the same time 106. Also the limit of infinite Bond number is shown, weallow forallparticlestomovewhichcanleadtoapar- leading to pure ballistic deposition [42] well studied al- tialcollapseofthestructuresduetogravity[8,9,13,40]. ready in the past. For small Bond numbers, here repre- The structure is deposited ona flatsurface,i.e. a wallat sented by Bo = 0, the system typically reaches a ran- g the bottom. We use periodic boundaries in horizontal domclosepackingwhichalsohasbeenstudiedintensively directiontoavoideffectsofsidewalls,likeJansseneffect. in the past. Note that our case of monodisperse parti- During this process the time interval between successive cles typically leads in dense packings to crystallization depositions crucially determines the structure and den- effects which could be avoided by using a small poly- sityprofilesofthefinalconfigurations. Herewewillfocus dispersity. As our focus in this paper is on the looser onthe twoextreme casesofverylargetime intervals,i.e. structureswherethiseffectisnotveryimportantwepre- thesystemcanfullyrelaxaftereachdepositionofasingle ferthe monodispersesystemtokeepthe modelassimple grain,andvanishingtime interval,i.e.the collapseofthe as possible. In the intermediate range of Bond numbers systemshappensafterthedepositionprocessiscomplete. the density varies between the two limiting values. Inthe firstcasethe intervalis chosenlargeenoughto let Plotting the density profile depending on the vertical the system compactify and relax due to the additional positiony (Fig.2)providesamorequantitativeanalysis. weightofthe depositedgrain. This is verifiedonthe one Itcanbeseenthatinthetwolimitingcases(Bo =0and g handbycheckingthatthefinaldensityisindependenton Bo →∞)the densityis constant. Forthe infinite Bond g the time interval and on the other hand by monitoring number this can be explained easily as no collapse at all the dynamics of the process. Having no time between occurs and the density profile is that of a ballistic depo- depositions in practice means that first we perform pure sition and thus constant [41, 42]. For the non-cohesive ballistic deposition [41, 42], and then switching on the case a close packing is expected, also leading to a con- full particle dynamicsleading to a collapseofthe system stantdensity. Thiswillbediscussedagaininmoredetail due to gravity. Experimentally, the two cases can be re- laterinthispaper. Intheintermediaterangethedensity alized in a Hele-Shaw cell [43–45] which can be tilted to decreases with increasing height. This is a result of the effectively change gravityIn the slowdeposition process, generationprocesswherethefragilestructureispartially simply the cell is slowly filled in an upright position so collapsed due to the weight of the added particles which that full gravity acts on the grains. In the other case happensdiscontinuouslyinrelativelysmallavalanchesas the Hele-Shaw cell will be almost horizontal, so that the will be discussed in more detail in the next section (sec. grains can be filled in with nearly vanishing gravity,and IV). then the cell is tilted so that gravitycan fully act on the Knowingthat the density depends onverticalposition grains, leading to an abrupt collapse of the structure. a general dependence of the total density on the Bond number cannot easily be defined. Instead, for a given system size as in Fig. 2 the density at a fixed position III. DENSITY PROFILES WHEN GRAVITY can be measured. In Fig. 3 the averaged density in the ACTS DURING DEPOSITION lower half excluding the region very close to the bottom is shown versus the granular Bond number. The den- In this section we analyze the density profiles for the sity varies between the two limiting cases Bo = 0 and g case of large enough time intervals between successive Bo → ∞. Note that the Bond number is plotted in a g depositions to allow the systems to relax under the ef- logarithmic scale, i.e. to see substantial changes of vol- 3 Bo g =0 Bo g =1e2 Bo g =1e3 Bo g =1e4 Bo g =1e5 Bo g =1e6 Bo g infinity FIG. 1: (Color online) Final structures achieved by the deposition/collapse process for different granular Bond numbers Bog. Inaddition totheparticles compressiveforces areillustrated byred (darkgray) lines connecting thecenterof masses between theparticles. In thecase of Bog →∞ no forces are present as it is realized in the simulations byswitching off gravity. 1 Bo=0 Bo=0 0.8 g g 0.8 Bo=102 g ν0.6 BBoogg==110034 wer half00..67 BBoogg=→0∞ 0.4 Bog=105 νlo Bo=106 0.5 g 0.2 Bog→∞ 0.4 0 0.3 0 100 200 300 400 10 100 1000 10000 1e+05 1e+06 y Bo g FIG. 2: (Color online) Density profiles for different granular FIG. 3: (Color online) Average volume fraction νlowerhalf Bond numbers Bog (cf. fig. 1). Here, the volume fraction ν depending on granular Bond number Bog. The density is is plotted. In this case the volume fraction is measured in averaged in thelower half of thesystem excludingtheregion thinslicesofgivenwidth(here: 3.97particleradii)atvarying very close to the bottom to avoid border effects (here we height y. For Bog =0 no cohesion is active and the random excluded the region below the height of 50 particle radii so close packing is reached. In the limit Bog → ∞ the system that clearly boundary effects are removed for all curves, cf. does not collapse at all, and the simple ballistic deposition fig. 2). The volume fractions vary between the two limits case [42] is obtained. given by random close packing (Bog = 0) and pure ballistic deposition (Bog →∞). ume fractionthe cohesiveforceor the gravitationalforce have to be changed by orders of magnitude. Particles has been found [49, 50]: varying the friction coefficient withsimilargravityandcohesiveforcewillshowthesame on a logarithmic scale leads to a variation between the typical behavior. As typically both forces depend on the values 0.84 for the packing fraction of a random close sizeoftheparticlesitappearstobenaturaltocharacter- packingandthe value0.77forinfinitely largefrictionco- ize the behavior of granular matter and powders by the efficient(intwodimensions,inthreedimensionsbetween grainsize. For non-cohesivematerial recent experimen- 0.64 and 0.55). In the cohesive case as discussed here tal, numerical and theoretical studies [49–52] investigate this range of accessible volume fractions is much higher theinfluenceofthefrictioncoefficienton,e.g.thevolume andlimitedbythepreparationprotocol,i.e.inthispaper fraction. A similar behavior as found here for the co- by the ballistic deposition. This limit of course can be hesivematerialwhenvaryingthegranularbondnumber, changedwhenchangingthe preparationprotocol,e.g.by 4 introducing a capture radius (cf. sec. VI). 0.8 H=397 H=797 H=147 0.7 ν 0.6 0.5 0.4 FIG. 5: Particle trajectories of particles of the small system (H = 147, cf. Fig. 4) for the whole deposition/collapse. For 0 200 400 600 800 better visibility only each 5th particle’s trajectory is shown, H-y i.e. the trajectories of 640 particles (instead of all 3200 par- ticles). Viewing the total system (and on the right more de- tailed when zooming in) illustrates that parts of the system FIG.4: (Color online)Illustration oftheeffect ofsystemsize movecollectivelydownwardsaccompaniedbyasidewardsmo- for intermediate density range (Bond number Bog = 103). tionorrotation. Whenzoomingintheindividualtrajectories PlottingthedepthH−y measuredfromthesurfaceH ofthe can be identified which are composed of the sum of paths final packings smaller systems show thesame profile as large duringall thesmall avalanches experienced by theparticle. systems. cles are plotted for a relatively small system of height For allresults presentedabovethe totalsystem height H = 147 consisting of about 3200 particles (for better H was fixed, i.e. the deposition process stops when no visibility only each 5th trajectory is shown, i.e. the tra- more particles can be deposited below a specified value jectories of 640 particles, instead of all 3200 particles). H. Whencomparingdensityprofilesfor differentsystem The avalanches are a collective motion of parts of the heights H plots depending onthe verticalposition y will system. This mainly downwards motion is accompanied show different densities. A scaling canbe achievedwhen by a sidewards motion or rotation. When zooming in plottingthedensityversusthedepthH−y asillustrated individual trajectories can be identified. These trajec- in Fig. 4. This means the upper partof the large system tories represent the motion of each particle during de- isdepositingandcollapsinginthe samewayasthe small position/collapsing. Thus, they show the paths that a systemwhile additionallyleadingto a further collapseof particle experiences in all avalanches at different times. the structure deposited previously below, accompanied Neighboring particles can have very similar trajectories, by a downwards motion of the whole upper part. Obvi- i.e.theybelongto the samesetofavalanchesatdifferent ously the slow depositionprocessguaranteesthat inertia times. is not important(cf. sec. VI). In Figure 6 we show the size of avalanches depending The specific behavior of the density profiles shown in on initial and final vertical position. This size is mea- this section results from a deposition process combined sured by ∆y, the total downwards displacement of the with a collapse of the current structure due to gravity. particle after its deposition, i.e. initial position minus fi- The deposition is characterized by the number of de- nalposition. Thisrepresentsforeachparticle the sumof posited particles per volume, which we call “deposition all avalanches occurring during the generation process, density”andwhichhereisnotconstant(sec.V). Thecol- resulting in as many data points as particles in the sys- lapsehappenssuccessivelyinrelativelysmallavalanches, tem. In Fig. 6 this data is averaged in bins of size two analyzed in detail in the following section. In section V particle diameters. The fluctuations within each bin are we will show that these avalanches can be used to relate shownby theverticalerrorbars. Bothcurves(for y and the final density profile to the “deposition density”. i y ) can be relatively well approximated by parabolas: e ∆y(y )=a′+b′y +c′y2, ∆y(y )=a+by +cy2 (1) IV. ANALYSIS OF THE AVALANCHES i i i e e e DURING DEPOSITION/COLLAPSING It is obvious that both curves cannot obey exactly the parabolic behavior as y and y are related by y (y ) = i e e i Typicallythecollapsingofthestructures,asmentioned y +∆y(y ). However,in the cases presented in this sec- i i earlier,happensdiscontinuouslyinsmallavalanches. As tion, obtained by slow deposition, the value of ∆y is rel- these avalanches are important also for the final den- atively small compared to y so that y (y ) is very close i e i sity profiles (see sec. III) their characteristics is stud- to a straight line, leading only to a very small horizon- ied in detail in this section. To illustrate the nature of tal shift. This behavior is typical for intermediate Bond these avalanches in Fig. 5 the trajectories of the parti- numbers whereas in the limiting cases no noticeable de- 5 0 axes by the system height H. From this scaling one can deduce the system size dependence of the pre-factor of the quadratic term in eq. (1). The scaling becomes: -20 ∆y(y,H)=H ·f(y/H)∝H ·(y/H)2 ∝1/H (2) ∆y(y) ∆y-40 i when assuming that ∆y ∝ y2 (parabolic behavior, see eq. 1). This 1/H dependence could be verified by fitting the curves in Fig. 7. Note that the parabolic shape was -60 ∆y(y) alsofoundwhenvaryingthefrictioncoefficientµandthe e rolling friction coefficient µ . r -80 0 200 400 600 800 y, y i e 1000 FIG. 6: (Color online) Size of avalanches depending on ver- cy tical position for Bog = 103. Here the size is measured by en100 u ∆y, the average total downwards motion of a particle after q e deposition(initialpositionminusfinalposition). Onthehor- r f all slides izontalaxistheinitialpositionyi (bluesquares)andfinalpo- without bottom and top sition ye (red circles) are plotted. This leads to two slightly 10 shifted curves as ye <yi. The solid lines represent parabolic fits∆y(yi)=−1.8−0.31yi+0.00035yi2 (black,full line) and ∆y(ye)=−10.2−0.31ye+0.00038ye2 (violet, dashed). Addi- 0 20 40 60 80 100 tionallyshownisafitby∆y(ye)=−aye(1−ye/H)predicted |∆y| by the considerations in sec. VII leading to a ≈ 0.39 (green, dashed-dotted). FIG.8: (Coloronline)Thehistogramofthesizeofavalanches |∆y| for Bog = 103 follows basically a Gaussian. Deviation pendence of ∆y on the vertical position could be found. from this behavior can almost fully be suppressed when re- For Bo = 0 a small constant value, below the particle moving thebottom and top part of thesystem. g diameter (around 1.5 particle radii) is observed. In the case Bo →∞ no collapse happens, i.e. all ∆y =0. Whereas the averageofthe avalanchesize ∆y as func- g tion of the vertical position shows a parabolic profile of 0 reasonable quality, there are of course large fluctuations around this value. In Fig. 8 we show the distribution of -0.02 H=797 theavalanche sizes(here|∆y|)forthe entiresystem,i.e. H=397 independentontheverticalposition. Whenremovingthe upper and lower part of the system to decrease bound- H-0.04 ary effects we obtain a Gaussian distribution, i.e. we get y / an estimate of a typical avalanche size. This typical size ∆-0.06 decreaseswith increasing Bondnumber, and in the limit of Bo →∞, where no avalanches occur, it vanishes. In g the limit of Bo =0 (no cohesion) the behavior is differ- -0.08 g ent, an exponential decay is obtained (Fig. 9). Here the boundaries have no effect, i.e. we get the same behavior -0.10 0.2 0.4 0.6 0.8 1 whenremovingtheupperandlowerpartofthesystemas y /H donepreviously. ForthisBondnumbertypicallythesur- i face of the structure during deposition grows relatively flat, so that large |∆y| are unlikely as expressedby the FIG. 7: (Color online) Collapse of the size of the avalanches exponential decay. Due to the monodispersity this sur- for two different system sizes can be obtained scaling both faceislocallyalmostaflatcrystallinesurfacewithheaps axes by the system height (here: Bog = 103). Under the which consists of a few particles only, in most cases one assumption of a parabolic profile this scaling leads to a 1/L particle. When a particle is deposited on a one particle dependenceof c (pre-factor of quadraticterm in eq.1). heapitrollsofftoresteventuallyasa“crystalline”neigh- borasidetheparticle,resulting in|∆y|betweenoneand The parabolic behavior can be reproduced also for two. Thisleadstotheverysmallrange of|∆y|withcon- othersystemheights. InFig.7twodifferentsystemsizes stantprobabilityinFig.9. Takingaslightlypolydisperse againfor Bo =103 are showncollapsedby scaling both system this region would disappear. g 6 givenheighty (L widthofthetwodimensionalsystem i,e x all slides in units of particle radii): without bottom and top 1000 yi,e ′ ′ N (y )=L dy ρ (y ) (3) y d,f i,e xZ d,f c 0 n e100 u The final position y of particles can be related to the q e e position y of deposition by the avalanche profile ∆y: r i f 10 y (y )=y −∆y(y ), or y (y )=y +∆y(y ) (4) i e e e e i i i In this notation∆y is negative as the motion of the par- 1 ticles is downwards (due to gravity). Therefore, y is 0 2 4 6 8 10 12 14 i |∆y| larger than or equal to ye. As particles are never de- stroyed the number of particles deposited up to a given height, N (y ), will stay the same, but shifted to a lower d i FIG. 9: (Color online) The histogram of avalanches|∆y|for height, Nf(ye), where yi and ye are related by eq. (4). Bog = 0 basically shows an exponential decay. Deviation Together with eq. (3) this leads to: from this behavior can be found for |∆y| between 1 and 2 particle radii where the probability is about constant. This yi(ye) ′ ′ effect cannot be suppressed when removing the bottom and Nf(ye)/Lx =Nd(yi(ye))/Lx = dy ρd(y ) (5) Z top part of thesystem (as, e.g. for Bog =103). 0 ≡G(yi(ye)) | {z } In this section we studied the collapse of the struc- ThisrelatesNf tothe depositiondensitywhereaseq.(3) tures occurring in small avalanches we analyzed statis- relates Nf to the final density. The function G here is tically. We suggest to characterize these avalanches by formally introduced for the integral as abbreviation, by their“size”,showingatypicaldependenceonverticalpo- derivationofGthedensityisretrieved. Thefinaldensity sition, a parabolic shape for the specific systems investi- can be obtained by derivation of Nf/Lx using eq. (5): gatedinthis section. Inthe followingsectionwe will use d N (y ) d dG(y )dy this characteristic behavior to be able to derive the final ρ (y )= f e = G(y (y ))= i i f e i e density profile from the “deposition density”. In section dye Lx dye dyi dye VI the same concept will be shown to be applicable also dy i = ρ (y (y )) (6) for other protocols of generating loose structures. d i e dy e d∆y(y ) e = ρ (y (y )) 1− d i e (cid:18) dy (cid:19) e V. THEORETICAL ANALYSIS OF THE AVALANCHES The deposition density ρ (y (y )) in principle can be ex- d i e ′ presseddirectly byy introducingρ (y ). As usually the e d e In the previous sections we mentioned that the dy- functional behavior of both functions is not known, but namics leading to a final configuration is determined by only valuesfor specific y andy ,the transformationcan i e smallavalanchesoccurringduringthedepositionprocess. be done for each point by simply using eq. (4), i.e. re- All these compaction events contained in the function placing each y by y = y +∆y(y ). Summarizing, to i e i i ∆y(y ), which is given by the difference between the ini- calculate the final density profile one needs to know the e tialpositiony andfinalpositiony . Notethat∆ycanbe depositiondensityρ andtheavalancheprofile∆y. Note i e d plotted (e.g.fig. 6) as function ofthe final positiony or that the avalanche profile dependence on both y and y e e i alternatively as function of the position of deposition y . is needed, which can be calculated from each other for i Theaimofthissectionistorelatethefinaldensityprofile some cases as shown later. For experimental situations to the dynamic process of deposition and collapse by us- these quantities are not known. However, the relation ing∆y(y ),showinghowtheavalanchesproducethefinal between ρ and ρ (eq. 6) can be used to calculate the e f d density ρ (y ) from the deposition density ρ (y ). The deposition density from the final density in the slow de- f e d i deposition density is defined by the number of particles positionlimit,whenassumingaparabolicprofileasfound deposited within a volume. As the structure collapses in the simulations before. between the depositions the deposition density is not in- In figure 10 we use eq. (6) to calculate the final den- dependentonthecollapsing,anditispossiblethatat(al- sity from the deposition density by using the parabolic most) the same position several particles are deposited. fit for ∆y (fig. 6). In practice first the deposition den- Thus, locally within a fixed volume even more particles sity curve is shifted on the horizontal axis by y = e could be deposited than typical for a dense packing. y −(a′+b′y +c′y2),thenmultiplyingthedepositionden- i i i We firstcalculatethe number N ofparticlesupto a sity withtherighthandsideofeq.(6),1−(b+2cy ),i.e. d,f e 7 0 0.8 0.6 -5 d ν 0.5 0.7 -10 0.4 ∆y(y) νf 0 200 400 600 800 ∆y -15 i y 0.6 -20 0.5 -25 ∆y(ye) 0 200 400 600 800 -300 100 200 300 400 y y, y i e FIG. 10: (Color online) Using the parabolic approximation FIG. 11: (Color online) Size of avalanches dependingon ver- (fig. 6) for theaverage avalanches thefinaldensity (here vol- tical position for Bog = 102. Here the size is measured by ume fraction ν ) can be calculated from the deposition den- f ∆y, thetotal downwards motion of the particle after deposi- sity (here shown: volume fraction νd for Bog = 103, inset). tion (initial position minusfinalposition). Onthehorizontal Therearestrongfluctuationsinthedeposition densitywhich axis the initial position yi (blue squares) and final position are induced by the irregularity of the avalanches. To obtain ye (redcircles) are plotted. Thelines representtheparabolic awmiasttmhchoeixontpghonrceeulnarttvie0v.e1wl5ye)wuteoslelcwaalictfiuhtlattfhueentmchteeioafinsnua(rhlededreden:esniptsyiotwy(beprlaroclakfiwlleinfifoetr), fi∆tys(∆yey)(=yi)−=4.7−−1.70.−230y.e24+yi0.+0000.5090y0e256(yfui2ll(dlianseh,ebdl,avckio)l.et) and sufficiently large y (except close to thebottom) 1 0.8 0.7 usingthederivativeof∆y(ye) whichisalinearfunction. 0.9 d If the deposition density would be constant this would ν 0.6 leadto alinearprofileforthe finaldensity. However,the 0.8 0.5 deposition density is not constant, explaining the non- linear behavior for the final density. Additionally the νf 0 100 2y00 300 400 deposition density shows strong fluctuations, but by as- 0.7 suming the avalanches to follow the averaged parabolic behaviorthe correspondingfluctuations inthe avalanche 0.6 profile are not included. The calculated curve matches relativelywelltheprofilemeasuredinthesimulationsfor 0.5 sufficiently large values of the vertical position. Close to 0 100 200 300 400 y the bottom, however,the calculated curve deviates from themeasuredone. Inthisregionthedepositiondensityis verysmall,i.e.almosttheoneofpureballisticdeposition. FIG. 12: (Color online) Using the parabolic approximation Thiscanbeunderstoodasthesystemneedstogainasuf- (fig. 11) for the average avalanches the final density can be ficientamountofweightforthe collapseto start(cf. also calculated from the deposition density (here shown: volume sec.VI). This shouldcorrespondto ahigher initialslope fraction for Bog =102). There are strong fluctuations in the of∆y(y )whichisnotreflectedintheparabolicapproxi- e depositiondensitywhichareinducedbytheirregularityofthe mation(eq.1). Inthisregionahigherordertermswould avalanches. As before (for Bog = 103), to obtain a smooth be necessary for reproducing also the system bottom. curve we use a fit function (here: power law fit, with expo- The same analysis has been done also, e.g. for Bo = nent0.13) tocalculatethefinaldensity(blackline)matching g 102 as shown in Figs. 11 and 12. In this case the de- relativelywellforsufficientlylargey (exceptclosetothebot- tom). position density shows somewhat lower fluctuations as for Bo = 103 (Fig. 10). To quantify this we esti- g mated the fluctuations of the deposition density at ver- tical position y = 200 for both cases. For Bo = 102 age profile. Away from this limit, but still close enough g we obtained around 15% whereas we estimated around that the density profile is very similar to the Bo → ∞ g 20% for Bo = 103. For the case of Bo → ∞ there is case,as e.g. for Bo =104, very largefluctuations in the g g g no avalanching at all (cf. sec. IV), and trivially the fi- avalanche profile are observed. Thus, the parabolic pro- nal density equals the deposition density. This is very filecannoteasilybeidentified. Stillthetheoryworkswell similar for very large Bo , but as some avalanches occur as the final density is very close to the deposition den- g there are some relatively small fluctuations in the aver- sity so that even a very inaccurate fit for the avalanche 8 profile does not affect the calculated density profile too Similarly as we did before we analyze the size of the much. In the limit of Bo = 0 the avalanche profile is avalanches ∆y as defined in sec. IV. Figure 15 shows g a constant (cf. sec. IV), i.e. all grains are slightly shifted a linear dependence of ∆y on either y and y . The fit e i downwardsby the same amount (except boundary effect parameters of the two lines can be related to each other atthe bottom). As then the derivativevanishes the final by the relation between y and y (eq. 4). Assuming i e ′ ′ ′ density equals the deposition density. ∆y(y ) = a−by and ∆y(y ) = a −by the values a e e i i ′ Here, we showed how the parabolic avalanche profile and b can be calculated from a and b (see app. A) as: canbeusedtocalculatethefinaldensityprofilesfromthe b a deposition density in the case where gravity acts during b′ = , a′ = (7) 1+b 1+b deposition. In the next section the same concept will be used to the simpler case of collapse of the system after The vertical dependence of ∆y can be used similarly as deposition is complete. These two cases could be then before to calculate the final density from the initial den- related in section VII. sity by using eq. (6). The calculated density profile us- ing this linear dependence reproduces the obtained final density profileverywellasshowninFig.14. In this case VI. COLLAPSE AFTER DEPOSITION the agreement is better as now the initial density is not COMPLETE fluctuating very much in contrast to the cases discussed in sec. V. The density increase ∆ρ (or volume fraction increase ∆ν) can be directly calculated by the constant Intheprevioussectionsweinvestigatedthecasewhere slope of ∆y(y ): gravityacts duringdeposition,leadingto relativelycom- e plex shape of the density profiles and a parabolic char- ∆ρ ∆ν d∆y(y ) e acteristics ofthe avalanchesize. For this casewe showed = =− (8) ρ ν dy that these avalancheprofiles can be used to relate the fi- i i e naldensityprofiletothe depositiondensity. Inthissec- For the same parameters (Bo = 103) we studied the g tion we will analyze the case when the particles are first effect of the system height H on the density increase deposited,thengravityisswitchedonandthestructures while still keeping the initial density fixed (fig. 16). A collapse. This case is even simpler and can later be used logarithmic fit matches the data best. This fit certainly to understand the more complex system studied before. cannotcontinuetoinfinityasthereisalimitfortheden- In this case the initial density ρi characterizes the sys- sityρmax givenbytherandomclosepacking(seealsofig. tem(insteadofthedepositiondensityasinthepreviously 3), leading to a (∆ν/ν ) of 1.19(≃ρ /ρ −1). i max max ini discussed situation). Using the off lattice version of bal- Usinginitialcaptureradiiasdescribedabovewestudy listic deposition as presented in Refs. [41, 42] with stick- the influence ofthe initialdensity onthe relativedensity ingprobabilityone,verticallyfallingparticlesstickwhen increase ∆ν/ν (fig. 17). We could obtain the best fit i they touch an already deposited particle. This leads to when using a power law with exponent of about 1.64. a fixed initial density. Lower densities can be obtained We showed in this section that the linear avalanche by using acapture radiusrcapt i.e.particles stickto each profile is a characteristic feature of compacting from a otherwhen they arewithin acertaindistanceduring the depthindependenttoadepthindependentstructure,ob- falling of the depositing particle. More precisely: When tained here for systems generated by ballistic deposition thedistancebetweenthecenterofmassesoftwoparticles collapsing due to gravity. More complex avalanche pro- isbelow2·rcapt theparticlesstickandthefallingparticle fileswithnon-constantderivativewilltransformhomoge- is pulled along the connecting line towards the already neous structures into inhomogeneous structures. Thus, deposited particle. This capture radius is a measure for we expect the linear profile to be obtained in all cases thedistancebetweenthebranchesofthe depositandthe where a homogeneous initial system compacts to a ho- resulting density is inversely proportional to rcapt [32], mogeneousfinalsystem. Thesehomogeneouscompaction rcapt = 1 gives the original method. The resulting ini- processes are investigated in different research areas as tial structures are shown in Fig. 13. These structures e.g. discussed in Refs. [53–57]. In addition in the next obtained with different capture radius will be later used sectionwe will show that also for the more complex pro- to study the influence of the initial density. cess when gravity acts during deposition (sec. III) this First we will investigate the behavior using rcapt = 1. linear profile can be used to derive the parabolic profile Figure 14 shows the density profile before the collapse of the avalanches. which is the same that we got in the limit of Bo → ∞ g in sec. III, also independent on vertical position. After this deposition is complete gravity is “switched on” and VII. RELATION BETWEEN DEPOSITION the structure abruptly collapses. Here we choose a Bond UNDER GRAVITY AND SWITCHING ON number of Bo = 103. This leads to a final structure GRAVITY AFTER DEPOSITION g with higher density, in this case also independent on the verticalposition(Fig.14). Asnoparticlesareaddedafter For the very fast process a linear profile for ∆y de- the initial deposition the final system height is lower. pending on vertical position has been found (cf. fig. 15) 9 r =1 r =1.25 r =1.5 r =1.75 r =2 r =2.25 r =2.5 capt capt capt capt capt capt capt FIG. 13: Initial structuresgenerated by ballistic deposition with increasing captureradius rcapt. initial density 0.6 final density calc. from ∆y 0.5 0.4 ν 0.3 0.2 0.1 0 0 100 200 300 400 y FIG.14: (Coloronline)Theinitialandfinaldensityareabout constantwhendepositingfirstandthencollapsingthesystem FIG. 15: (Color online) The linear dependence of the (Bog = 103). Using the linear dependence of the avalanches avalanche sizes ∆y(yi,e) explains the homogeneous density ∆y ontheverticalposition (seeFig.15)thefinaldensitycan increase as seen in fig. 14. The linear fits are ∆y(ye) = be calculated from the initial density using eqs. (4) and (6). 2.2−0.62ye and ∆y(yi)=1.9−0.39yi. The results support the analytical considerations. profile for ∆y depending on vertical position (cf. figs. whereas the slow deposition limit shows a parabolic 6,11). In this section we will discuss how a relation be- 10 1 0.8 0.6 1 ρ0 h0 ρ1 h1 νi ν / ∆ 0.4 0.2 0 500 1000 1500 1 ρ0 hh00 ρ1 h1 H 2 ρ1 h1 ρ2 h2 FIG. 16: Dependence of volume fraction increase ∆ν/νi on system height H when first deposited and then collapsed. A logarithmic fit matches the data best (here y = −0.94 + 0la.r2g6elsnt(xp)o)s.sibTlehevalilmueitfoorf r∆aνn/dνoimofcl1o.s1e9,pawchkiicnhg wdeilfilnbees athpe- ii+1 ρρii−1 hhii−1 ρρii+1 hhii+1 proached for infinitesystem heights. 3 n ρn−1 hn−1 ρn hn 2.5 FIG. 18: Sketch illustrating the procedure of depositing the grains slice by slice. The first slice deposited is compacted by internal collapse. The same is true for each “freshly” de- νi 2 ν / posited slice. The slices below are compacted by the added ∆ weight of the slices above. Periodic boundary conditions in 1.5 horizontaldirectionareimposed(illustratedbydashedlines). The figure also illustrates the definition of the symbols used here. The n slices are numbered from 1 to n. A slice i col- 1 lapsesfrom ρi−1 toρi whileitsheightdecreases from hi−1 to hi, where hi =ρi−1/ρihi−1. 0.1 0.15 0.2 0.25 0.3 0.35 0.4 ν i Letusfirstconsidersystemscomposedofasmallnum- FIG. 17: Dependence of volume fraction increase ∆ν/νi on ber n of slices. The case n=1 (one slice) is the same as volumefractionνioftheinitialsystem(systemfirstdeposited discussed in the previous section: The system collapses and then collapsed). Different densities could be reached by “internally”leadingtoanincreaseofthedensityfromρ 0 increasingthecaptureradiusforballisticdeposition (fig.13). toρ whiletheheightreducesfromh toh . Herewede- 1 0 1 Apowerlawfitwith exponent1.64fitsrelativelywell(power note the slice number as 1 (cf. Fig. 18). As shownin the law fit results in y=0.158x−1.64). previous section the avalanche sizes have a linear profile ∆y(y(1)) = S y(1). S is the slope in slice 1, and is the e 1 e 1 same for all freshly deposited slices when the height h tween both can be established. By this relation also the 0 is kept constant. The vertical position y(1) within slice parabolic profile is put onto a more fundamental basis e like the linear profile for the homogeneous collapse. 1 is measured from its bottom (ye(1) = 0...h1). This notation will be used in the following for each slice i: Let us imagine depositing particles slice by slice as (i) sketched in Fig. 18. The slices are thin parts of the sys- ye =0...hi. The case n=2 (two slices) means adding teminverticaldirectionspanningthefullsystemwidthin an additional slice to the case n = 1. Then the lower horizontal direction. They can be considered as systems slice(slice2)experiencesanadditionalcompactionbythe with very small initial height h . In each slice the depo- added weight expressed by the corresponding avalanche 0 sitionwillbeimmediatelyfollowedbythecollapse. How- size C y(2) assuming a linear behavior for this relatively 2 e ever,there will be not only an“internalcollapse” within fast process similar as for the internal collapse. This is the“freshly”depositedslice,butalsoacompactionofthe justified at least for the limit of small slices later consid- slices belowdue to the additionalweightofthe “freshly” eredinthissection. Theupperslice(slice1)willbecom- deposited slice. pacted internally and additionally will move downwards