Computer simulation of random loose packings of micro-particles in presence of adhesion and friction
CComputer simulation of random loose packings of micro-particles in presence ofadhesion and friction
Wenwei Liu, Shuiqing Li ∗ , Sheng Chen Key Laboratory for Thermal Science and Power Engineering of Ministry of Education,Department of Thermal Engineering, Tsinghua University, Beijing 100084, China (Dated: October 9, 2018)With a novel 3D discrete-element method specially developed with adhesive contact mechanics,random loose packings of uniform spherical micron-sized particles are fully investigated. The resultsshow that large velocity, large size or weak adhesion can produce a relatively dense packing whenother parameters are fixed, and these combined effects can be characterized by a dimensionlessadhesion parameter ( Ad = ω/ ρ p U R ). Four regimes are identified based on the value of Ad : RCPregime with Ad < ∼ .
01; RLP regime with ∼ . < Ad <
1; adhesion regime with 1 < Ad <
Ad >
20. Force distribution of these adhesive loose packings follows P ( f ) ∼ f θ for small forces and P ( f ) ∼ exp − βf for big forces, respectively, which shares a similar formwith that in packings without adhesion but results in distinct exponents of θ = 0 . β = 0 . Jammed particle packings have been studied to under-stand the microstructure and bulk properties of liquids,glasses and crystals [1, 2] and frictional granular mate-rials [3, 4]. Two packing limits have been identified fordisordered uniform spheres: the random close packing(RCP) and random loose packing (RLP) limits [1, 5–11]. The upper RCP limit is reproduced for frictionlessspheres at volume fractions φ ≈ .
64 and has been associ-ated with a freezing point of a first order phase transition[12–15], among other interpretations [2, 16]. Only in thepresence of friction, packings reach lower volume fractionup to the RLP limit φ RLP ≈ .
55 for mechanically sta-ble packings [6, 8, 11]. However, most packings of drysmall micrometer-sized particles in nature are not onlysubject to friction, but also adhesion forces as well. Forinstance, van der Waals forces generally dominate inter-actions between particles with diameters of around 10 m or smaller. In this case, the adhesive forces begin to over-come the gravitational and elastic contact forces actingon the particles and change macroscopic structural prop-erties [17, 18].Despite the ubiquity of adhesive particle packings inalmost all areas of engineering, biology, agriculture andphysical sciences [18–21], there have been few systemati-cally investigations of these packings [22–27]. The multi-coupling of adhesion, elastic contact forces and frictionin the short-range particle-particle interaction zone andtheir further couplings with fluid forces (e.g. buoyancy,drag and lubrication) across long-range scales make ithighly difficult to single out the effect of the adhesionforces, let alone to investigate the packing properties ∗ [email protected] experimentally. With the progress of computer simula-tion techniques, discrete element simulation has becomean efficient and accurate method to study the packingproblems of micron-sized particles [17, 22, 27], which arerather difficult to achieve in the experiment condition.In a dynamic packing process, particles will be arrestedand eventually form a packing when all the kinetic energyis dissipated after a series of collision. To better under-stand the dynamic behavior during the packing process,various dynamic models were proposed to describe theimpact process with a combination of quasi-static con-tact theories and dissipation mechanisms [17, 18, 28–30].While the dissipation mechanisms account for the dy-namic effects such as viscoelasticity, the quasi-static con-tact theories fundamentally describe the contact betweenparticles, of which Hertz contact theory is the most suc-cessful to characterize the relation of applied forces andcontact area of non-adhesive particles. However, as par-ticle sizes go down to micron scales, van der Waals ad-hesion becomes the most important interaction. Severalmesoscale models are put forward to describe the effectsof van der Waals adhesion on the elastic forces betweenstatic contacts of particles, among which the JKR (John-son, Kendall and Roberts), DMT (Derjaguin, Mullerand Toporov) and M-D (Maugis-Dugdale) models arewidely accepted ones [31–33]. Recently, Li and Marshall[34], and Marshall [35] developed a three-dimensional,mesoscopic discrete-element method (DEM) for adhesivemicron-sized particles based on the JKR model, whichhas been successfully applied to dynamic simulations ofmicron-particle deposition on both flat and cylindricalsurfaces with experimental validations [17, 36]. Withthis approach of adhesive DEM, both macroscopic andmicroscopic parameters during the dynamic process canbe predicted.Previous studies using a discrete element method(DEM) have found that the packing fraction decreases a r X i v : . [ c ond - m a t . s o f t ] N ov for adhesive micron-sized particles in a range of φ =0 . − .
622 with smaller sizes [22]. A similar resultof φ ≈ . − .
55 was found for 4 − µm particlesboth in simulation and experiment [27]. As for otherexperimental investigations, a random ballistic deposi-tion and fluidized bed technique was used respectivelyto produce the packing fractions of φ = 0 . − .
33 forboth uncompressed and compressed samples [24] and of φ = 0 . − .
41 with particle diameter of 7 . − . µm [23]. From these studies, we can conclude that formicron-sized particles, packing density has a positive cor-relation with particle size. Nevertheless, the effects ofother physical parameters, such as velocities, strength ofadhesion, frictions etc., on the packing process have littlebeen discussed or analyzed.In this paper, a prototypical packing system is intro-duced for the simulation of random loose packings of soft-sphere, non-Brownian , uniform adhesive particles with adiscrete element method. We explore the very low den-sity regime of small particles with van der Waals adhe-sion interaction by changing physical parameters. Theeffects of particle velocity, size and work of adhesion onthe packing properties are investigated via a dimension-less adhesion parameter. The paper is organized as fol-lows: the adhesive DEM simulation approach based onthe JKR theory is described in detail in Section 2.1; thesimulation conditions and parameters used in this workare given Section 2.2; and the results and discussionsabout the effects of the physical parameters on the pack-ing properties are included in Section 3. Section 4 givesthe conclusions. In a novel DEM framework specifically developed foradhesive grains [17, 18], both the transitional and rota-tional motions of each particle in the system are consid-ered on the basis of Newton’s second law. m d v dt = F F + F A ,I d Ω dt = M F + M A , (1)where v and Ω are, respectively, velocity and rotationrate of an individual particle, m is the particle mass, and I = (2 / mr p is the moment of inertia. In a vacuumsystem that may apply for problems of interstellar dusts,the fluid forces and torques acting on the particle, F F and M F , are ignored. F A denotes the collision and thevan der Waals adhesion forces in the equation for transla-tional motion. Meanwhile, in the equation for rotationalmotion, M A denotes the sum of the collision and van der Waals adhesion torques on the particle. They include F A = F n n + F s t s , M A = r p F s ( n × t s ) + M r ( t s × n ) + M t n , (2)where F n is the normal force including adhesivelyelastic contact force and damping force, F s is thetangential force due to the sliding friction, M r is therolling resistance and M t is the twisting resistance. r p is the particle radius. n , t s and t r are the normal,tangential and rolling direction unit vectors, respectively. Normal Elastic and Adhesive Forces
The normal force acts in the direction of the unit vector n which points parallel to the line connecting the centersof the two particles, denoted by i and j , such that n =( x j − x i ) / | x j − x i | . We consider two particles with radii r i and r j , elastic moduli E i and E j , and Poisson ratios σ i and σ j . An effective particle radius R and effectiveelastic moduli E and are defined by1 R = 1 r i + 1 r j , E = 1 − σ i E i + 1 − σ j E j . (3)The particle normal overlap δ N is defined by δ N = r i + r j − | x i − x j | , where x i and x j denote the parti-cle centroid positions. The particle normal adhesive andelastic forces are written together as F A = F n n , where a Lennard-Jones like formula for F A between two particlesis proposed (see more details in [18, 19]) F n F C = 4( aa ) − aa ) / ,δ N δ C = 6 / [2( aa ) −
43 ( aa ) / ] , (4)which are in turn based on the JKR theory, as a = R [ F ne + 3 πRω + (cid:112) πRωF + (3 πRω ) ] /E , here a ( t ) isthe contact region radius and ω is work of adhesion dueto van der Waals interactions (namely twice the surfaceenergy ω = 2 γ ) of two contacting particles [31]. The JKRmodel assumes that the adhesive force is only acting in-side the contact radius a , and results in a larger contactarea than the classic Hertzian contact, a > a h . Basi-cally, in contrast to the DMT model (Derjaguin-Mueller-Topprov), the JKR model is appropriate for compliant,adhesive particles for which the particle’s Tabor parame-ter is large than unity, implying the length scale of elasticdeformation is larger compared to the length scale of theadhesive force [18]. Then, the critical force and overlap, F C and δ C , and the equilibrium contact radius a aregiven by F C = 32 πωR, δ C = a / R , a = ( 9 πωR E ) / . (5)Typical values of ω are about 10 − mJ/m , from eithermeasurements [37] or Lifshitz theory’s predictions. Rolling Resistance
The rolling resistance exerts a torque on the particlein the M r t r × n direction, where t r is the direction of the“rolling” velocity. An expression for the rolling displace-ment of arbitrary-shaped particles is derived by Bagi andKuhn [38]. Taking the rate of this expression and special-izing to spherical particles of equal size yields an equationfor the “rolling velocity” v L of particle i as v L = − R (Ω i − Ω j ) × n . (6)An expression for the rolling resistance torque M r is pos-tulated in the form M r = − k r ξ · t r , (7)where the direction of rolling is t r = v L / | v L | and therolling displacement is ξ = (cid:82) tt v L ( τ ) dτ . Rolling involvesan upward motion of the particle surfaces in one part ofthe contact region and a downward motion in the otherpart. The presence of an adhesive force between thetwo contacting surfaces thus introduces a torque resistingrolling of the particles. An expression for the rolling re-sistance in presence of adhesion was derived by Dominikand Tielens [19, 39], which yields the coefficient k r as k r = 4 F C ( a/a ) / . (8)The critical resistance occurs when the rolling displace-ment magnitude, ξ = | ξ | , achieves a critical value,corresponding to a critical rolling angle θ crit = ξ crit /R .For ξ > ξ crit , the rolling displacement ξ in Eq. 6 isreplaced by ξ crit t r . It is noted, according to the mea-surement by atomic force spectroscopy, θ crit is around(0 . − . Sliding and Twisting Resistance
As aforementioned, both sliding and twisting are rel-atively rare for small adhesive particles rolling is gen-erally the preferred deformation mode for agglomeratesof adhesive particles [18, 39]. It is therefore desirableto introduce relatively simple expressions for sliding andtwisting resistance in the DEM framework. The stan-dard sliding model for the case without adhesion is thespring-dashpot model proposed by [41], for which thesliding force F s is given by a linear spring-dashpot, F s = − k t δ t · t s − η t v s · t s ( t s is tangential direction),when | F s | < F crit and by the Amonton friction expres-sion F s = − F crit when | F s | ≥ F crit . The twisting modelis given similarly as M t = − k Q (cid:82) tt Ω T ( τ ) dτ − η Q Ω T withΩ T = ( Ω i − Ω j ) · n the relative twisting rate [18]. Allthe model parameters ( k t , η t , k Q , η Q ) are chosen based on [36]. Here, a simple model proposed by Thornton [42]and Thornton and Yin [43], agreeing reasonably well withexperimental data, is introduced. In this model, the onlyinfluence of van der Waals adhesion on sliding force isto modify the critical force F c rit at which sliding occurs,which is given by F crit = µ f | F n e + 2 F C | , (9)where F C is the critical normal force given in Eq. 5 and µ f is a friction coefficient that is normally about 0.3. Whenparticles are being pulled apart, the normal force ap-proaches − F C at the point of separation, at which pointthe critical sliding force in Eq. 8 approaches µ f F C . Allthe values or ranges of µ f , θ crit and F C are selected ac-cording to the data from atomic force microscopy mea-surements [37, 40, 44].The same model with twisting resistance can be usedin the presence of adhesion, with the critical force F crit used to obtain M t = π aF crit . For twisting momentswith magnitude greater than M t,crit , the torsionalresistance is given by M t = − M t,crit . The generation of the packing starts with the succes-sively random free falling of 1000 uniform spheres withan initial velocity U at a height H under gravity. Thehorizontal deposition plane has two equal edges of length L along with periodic boundary conditions on both direc-tions. A stable packing structure is achieved when all theparticles are settled after a time long enough. Here, thefluid effect is filtered out by assuming packing under avacuum condition. More importantly, the gravitationaleffect with respect to particle inertia can be neglectedwhen the system satisfies F r = U / √ gH >
1, where
F r is the
Froude number (ratio of inertia to gravity). Toprecisely characterize the gravitational effect, we definethe relative velocity increment as ∆ U = ( U − U ) /U ,where U = (cid:112) U + 2 gH is the final velocity. For allruns in the numerical simulations, we ensure that ∆ U is less than 4%. Thus the gravitational acceleration dur-ing the dynamic deposition process can be ignored, in-dicating that the deposition velocity can be treated asthe same with the initial inlet velocity. After the pack-ing is formed, however, gravity still works and can’t beneglected. As a primary concern of our work, interstellardust particles always transport with a relatively large ve-locity before they form large aggregates. Therefore, theadhesive packings simply arise due to the competitionbetween the particle inertia and particle-particle interac-tions (e.g., adhesion, elasticity and frictions). Most im-portantly, the negligible gravitational effect distinguishesour system from that of [22, 27], which generate particlesrandomly in a box without touching each other and waitthem to deposit due to gravity.Before the simulation, a sensitivity analysis betweenthe cases L = 20 r p and L = 40 r p was conducted withparameters in Table I. As also shown in Table I, the dif-ference in φ between the cases L = 20 r p and L = 40 r p is negligible, indicating the L = 20 r p is large enough toreproduce bulk properties. Then, we set L = 20 r p . Thephysical and geometrical parameters used in the DEMsimulations are listed in Table II. In addition to chang-ing particle size, we applied different work of adhesionand initial velocities to explore their effects on adhesivepackings. TABLE I: Sensitivity analysis between cases of L = 20 r p and L = 40 r p . r p L U ω N Volume Coordination( µm ) ( µm ) ( m/s ) ( mJ/m ) Fraction Number1 20 0.5 30 1000 0.155 2.241 40 0.5 30 4000 0.157 2.255 100 0.5 30 1000 0.265 2.805 200 0.5 30 4000 0.270 2.7610 200 0.5 30 1000 0.372 3.1810 400 0.5 30 4000 0.375 3.19TABLE II: Parameters used in DEM simulations.Physical Parameters Values UnitsParticle Number( N ) 1000Particle Radius( r p ) 1,5,10,50 µm Particle Density( ρ p ) 2500 kg/m Work of Adhesion( ω ) 30,20,10,5,1,0.1 mJ/m Characteristic Length( L ) 20 r p µm Particle Injection Height( H ) 4 × L µm
Gravity Acceleration( g ) 9.81 m/s Deposition Velocity( U ) 0.5-10 m/s Typical packing structures and connecting networks ofloose and dense packings are shown in Fig. 1. The volumefraction is calculated from the vertically middle part ofthe packing (0 . h ≤ X ≤ . h , with h as packing height),avoiding the so-called wall effect from both bottom andtop of the packing structure. It has been well acceptedthat volume fraction increases with the increment of par-ticle size [22, 23, 27]. However, from Fig. 1 we can seethat not only the particle size ( r p ) but also the depositionvelocity ( U ) and work of adhesion ( ω ) both significantlyaffect mesoscopic packing structures. The increased U or decreased ω will result in a relatively dense packing, FIG. 1: (Colors online) Typical packing structures with vary-ing physical parameters. (a) stands for a basic loose packingwith r p = 1 µm , U = 1 m/s and ω = 30 mJ/m . The otherthree plots change only one physical parameter for each, whichare U = 2 m/s in (b), r p = 5 µm in (c) and ω = 10 mJ/m in(d) with other parameters fixed respectively. The right part ineach plot shows the connecting network of packing structures.Different color represents different coordination number Z . which is primarily due to the competition between par-ticle inertia and adhesion that will be discussed below.Also from Fig. 1, in a very loose packing the connectingnetwork turns out to be like chains with only two con-tacts for most of the particle while for a dense one thereare more small loops in the network.We then figure out the volume fraction of all the sim-ulation conditions under different sizes of particles ( r p ),deposition velocities ( U ) and work of adhesion ( ω ). Itis seen from the four sub-plots in Fig. 2 that, with in-creased deposition velocity, the volume fraction increasesand then stays at a value approaching 0.64, which isthe widely accepted as the upper limit of random closepacking ( φ RCP ) with uniform sphere particles. However,the extent of the increment is somewhat different undervariation of r p or under variation of ω . For instance,for a relatively small particles ( r p = 1 µm ), it is inter-esting that the very loose packing ( φ = 0 . U = 0 . m/s ) with strongadhesion ( ω = 30 mJ/m ). When the deposition velocitygrows to 10 m/s , the volume fraction reaches 0.582 withan increment of ∼ . ω goes downto 0 . mJ/m (nearly non-adhesive), the volume fractionhardly changes with U . On the other extreme, for muchbigger particles ( r p = 50 µm ), despite the changes of ei-ther U or ω , the volume fraction seems to converge to ahorizontal plane, which means the effects of both adhe-sion and particle inertia on the volume fraction is suffi-ciently small for big grains.A dimensionless adhesion parameter Ad = ω/ ρ p U R ,which interprets the balance between the interparticleadhesion and the particle inertia, has been successfullyapplied to physically lump together the effects of parti-cle size, velocity and work of adhesion [34, 45]. Fig. 3 FIG. 2: (Colors online) 3D plot of volume fraction with dif-ferent velocities and work of adhesion under the conditionsof different sizes of particles. Left top ( r p = 1 µm ), righttop ( r p = 5 µm ), left bottom ( r p = 10 µm ), right bottom( r p = 50 µm ). further extends the packing properties to both higherand lower Ad conditions than those in [45]. In the caseof Ad <
1, the volume fraction varies smoothly from ∼ .
50 (close to RLP) to ∼ .
64 (RCP), with parti-cle inertia dominating over adhesion. However, when1 < Ad <
20, an adhesion-controlled regime is observed,where the volume fraction decreases exponentially withincreasing Ad . Furthermore, with very high Ad ( > ω = 50 − mJ/m while the velocityis decreased to 0 . m/s with fixed particle size r p = 1 µm in the simulations in order to achieve packings with veryhigh Ad (= 80 − Ad >
100 and the lowest volumefraction is 0.128, which is very close to the well-knownlower bound of saturated sphere packings in d dimensions φ = 1 / d [45, 46].Regarding the densification with the decrease of Ad ,two potential ways of compaction may be the causes:gravitational (or other external forces) compaction andparticle inertial compaction. Following the work of [24],the pressure caused by these two compactions can be de-fined as p g = ρ p gr p and p i ≈ ρ p U . With the param-eters given in Table II, we get p g = (0 . ∼ . P a and p i ≈ (156 . ∼ . × ) P a . Thus, the com-paction caused by gravity is negligible and particle in-ertial compression is the most important way of densifi-cation, which further approves our primary assumptionthat gravitational effect can be negligible. Furthermore,from the definition of Ad we can see that the compressionresulted from particle inertia has already been consid- FIG. 3: (Colors online) Packing volume fraction as a functionof adhesion parameter. The black open squares are the caseswith strong adhesion. The red line is a fitting of φ versus Ad . ered and characterized with respect to adhesion. Whenparticles are being packed, the van der Waals adhesionforces tend to attract particles and make them stick to-gether. On the other hand, the particle inertia which hasa quadratic correlation with particle velocity will urgethem to move and impact with other particles. If adhe-sion is stronger than particle inertia, particles are morelikely to be caught at the incipient impacts and hardlymove or roll. Thus a loose packing structure is easierto form. With the increase of particle size or velocity,particle inertia will become much stronger than adhe-sion. As a consequence, more collisions will take placeand particles tend to explore more phase space and forma denser packing structure. Note that in some previousstudies the volume fraction is also related to the forceratio, or the so-called Bond number, which is defined as Bo = f c /f e [9, 47]. f c is the attractive cohesive forceswhich may come from van der Waals interactions, liq-uid bridge cohesion and so on while f e represents theexternal driven forces that could be gravity [9] or aver-age force due to external compression [47]. However, asdiscussed above, the balance of particle inertia and adhe-sion dominates our packings instead of gravity, meaningthat we can still achieve different packing structures bychanging velocities while Bo is unchanged. Therefore, forour inertia-adhesion controlled system, Ad is more appro-priate to characterize the mechanism than Bo , which isoften applied in gravity-cohesion controlled system.As mentioned above, the volume fraction seems to de-crease exponentially with the increasing Ad under condi-tions of Ad <
20, such that an approximately liner fit canbe well obtained by using the log-plot of φ − φ ∼ Ad , ex-cept those points with a certain fluctuation under extralower Ad parameter (e.g Ad < . φ is an asymptoticpacking fraction value without adhesion effect. It may beclose to the RCP limit of φ RCP = 0 .
64 since the inertialcompression exceeds the frictional force in current par-ticulate system. Therefore, a power-low relationship isproposed as φ RCP − φ = αAd λ , (10)The parameters α and λ can be obtained as α = 0 . ± .
005 and λ = 0 . ± .
016 by a best fitting. Thereforethe relationship is φ RCP − φ = 0 . Ad . , which isalso plotted in Fig. 3 as a red line. The inset of Fig. 3shows the log-log plot of this relation, from where we cansee the linear part of ( φ − φ ) ∼ Ad along with the tran-sition to the lower limit when Ad >
20. Basically, thisrelationship connects the lowest volume fraction and thewell-known random close packing limit smoothly via asimple parameter Ad , though the intrinsically underlyingphysics behind this relationship still need to be furtherexplored.Here more discussion is essential with respect to thescaling of Ad as it indicates a universal prediction thatthe limit of small Ad will lead to only RCP instead ofRLP. From the definition we know that Ad is inverselyproportional to the square of particle velocity whichwill result in the compression caused by particle inertia.Thus, the very small values of Ad are mostly obtainedwith large velocities instead of large particle sizes orweak work of adhesion and the compression will makethe packings denser and approach RCP. However, RLPis indeed reached in our results, which is only in asmall range of Ad (around 0.1). This contrasts with theprevious friction-dependent RLP results since friction isfinite and fixed in this work. With the tuning of Ad ,transition from RLP to RCP can also be realized withoutchanging friction. As the effect of friction on adhesivepackings is not a main concern of this study, we leave itto the future work. On the other hand, all the resultswith Ad > Ad = 1 that distinguishes adhesive packingsfrom adhesive-less ones. Therefore we can define fourregimes from the results in Fig. 3: RCP regime with Ad < ∼ .
01, RLP regime with ∼ . < Ad <
1, adhe-sion regime with 1 < Ad <
20 and an asymptotic regimewith
Ad >
20, where the packings begin to approachthe lower limit φ = 1 / that is not determined by anyphysical parameters like adhesion, inertia or friction. Coordination number ( Z ), defined as the number ofneighbours touching a given particle, is introduced tocharacterize the packing structures. We calculate Z byjudging whether the distance of the centers of two par-ticles is smaller than the sum of their radius. Thus this Z is the geometrical coordination number since it is de-fined with a criterion of geometric relation that doesn’tinclude any forces. The mean coordination numbers ofall the simulation cases are plotted in Fig. 4 as a function FIG. 4: (Colors online) Coordination number as a function ofadhesion parameter. The black dash lines represent Z = 4, Z = 6 and Ad = 1 respectively. The red line is a fitting of Z versus Ad . of Ad . We can see that for Ad < d + 1 ≤ Z ≤ d where d is the spatial dimension. As mentioned in thecomputational method section, all of our packings remainfrictional with fixed friction coefficient µ f = 0 . Z does not reach the isostaticlimit Z = 6 (of frictionless particles) at Ad <
1. Formicron-sized particles, friction is often coupled with ad-hesion and is complicated to single out. Thus we don’tpay much attention to the effect of friction on adhesivepackings for the present work. However, for
Ad > φ dependence. We figure out a fitting line of Z ver-sus Ad as Z = 3 . − . log ( Ad ). Similarly, there alsoseems to be an asymptotic value of coordination numberwhen Ad >
20, where the lowest Z reached in our sim-ulations is Z = 2 .
09 with
Ad > Z ∼ Ad curve, we also investigate thedistribution of coordination number. Typical cases arepicked among all the simulations to display the commonfeatures. As shown in Fig. 5, in the packing with the FIG. 5: (Colors online) Distribution of coordination numberfor different volume fraction and Ad . The solid lines are nor-mal distribution fit. lowest density, most of the particles have coordinationnumber in the range of 1 ∼ Z shifts to the right (larger values) and becomeswider. It’s also noted that the distribution of coordi-nation number is symmetric and shapes like normaldistribution. Hence we made a normal distribution fitand the results are listed in Table III. As we can see, thenormal distribution fits well with the data, indicating aubiquity of normal distribution of coordination numberin random packings of uniform spherical particles. TABLE III: Average packing properties of different Ad alongwith a normal distribution fit of Z .Case a b c d e f g Ad
48 9.6 4.8 2.4 0.96 2.5e-4 8.9e-5 φ Z Z Radial distribution function g ( r ), which describesthe probability of finding a particle at a distance r from a given one, has been widely used to depict theinterparticle correlations and characterize the packingdensity as well as the coordination number. Fig. 6 showsthe radial distribution function g ( r ) of typical packingslisted in Table III with various Ad . For loose packingswith φ < .
4, we apply dr = 0 . d p to calculate the g ( r ) while dr = 0 . d p for dense packings with φ > . g ( r ) with respect to the increasing of volume fraction. (i) Firstly, the g ( r ) of very loose packings (about φ < .
3) drops below g ( r ) = 1 instantly after the first peak and then gradually rises to the second peak whichis inconspicuous here but implies the second contactshell of typical structural features of random packingsof granular spheres. The feature that the lowest valleyof g ( r ) is close to r = d p indicates a very low probabilityto find a particle just outside the first contact shell.The gradual increase from the valley to the secondpeak directly describes the greatest possibility to finda particle at the distance r = 2 d p , accounting for thechainlike structure of very loose packings. After thesecond peak the g ( r ) becomes almost flat, indicating auniform probability to find a particle outside the distanceof r > d p . (ii) When the volume fraction grows larger(about 0 . < φ < φ RLP ), no valleys with g ( r ) < r > d p as well. Thischaracteristic reflects the uniform probability of findingparticles between r = d p and r = 2 d p for typical randomloose packings besides the strong peak at r = 2 d p whichgrows up with the increasing of packing density. (iii) Further, as packing becomes denser, approaching therandom close packing, valleys of g ( r ) < r = √ d p , r = √ d p and r = √ d p . The valley with the peakat r = √ d p depicts the adjacent local ordering thatfour particles form two edge-sharing in-plane equilateraltriangles, which is part of the densest packing structure.Despite the ordering at farther distance represented byother peaks, the not-sharp and wide peaks still can’tmake the packings reach the densest structure like facecenter cubic (FCC). In fact, it’s always improbable toobtain FCC with random packings since the possibilityof achieving both the local and global ordering struc-tures of FCC at the same time is almost zero. Only bysufficient compression or vibration can one get a ratherdense packing which might be beyond the random closepacking. However, no matter how dense the packingis, the g ( r ) should resemble that of stage 3 since g ( r )describes the spatial correlations which are definedgeometrically. Despite that the features of g ( r ) foundin our simulations are in consistent with results from[22], nonetheless, it should be noted that g ( r ) only hasrelation with volume fraction which can be tuned vianot only particle size [22] but also Ad . A broad rangeof volume fraction and various packing structures canbe achieved by tuning Ad . The three stages of g ( r )describe the evolution of packing density and structure ofadhesive micron particles, which bridges the loosest anddensest random packings through a simple parameter Ad . For loose packings of adhesive particles, especially forthe very loose packings with mean coordination num-ber 2 ∼
4, it is of great interest how the particles getmechanical equilibrium, or what the force distribution
FIG. 6: (Colors online) Radial distribution function for pack-ings with different Ad . is like. To investigate these properties, we measure theforce of every contact. Fig. 7 displays the signed nor-mal force network and distribution of four typical pack-ings we have achieved, which are also among the caseslisted in Table III. Here the signed forces means theycould be both attractive and repulsive forces betweentwo contact particles. Forces signed negative are repul-sive while the positive forces are attractive. It can bedirectly inferred from Fig. 7 that both attractive and re-pulsive forces appear for all the packing densities. Thisproperty distinguishes the packing of adhesive particlesfrom that of adhesive-less granular matters where all theforces are repulsive. Furthermore, for low packing den-sities ( φ = 0 . , . , . φ = 0 . φ on Ad . From Fig. 3 weknow φ increases with the decreasing of Ad , which rep-resents the strength of adhesion. A relative decreasing of Ad will then lead to the reduction of attractive forces.Another important property of the force network isthat, no matter what the packing density is, most of theforces are around zero which identifies the small forcesbehavior. In order to study this behavior, we measurethe unsigned force, namely the magnitude of the attrac-tive and repulsive forces. Fig. 8 shows the distribution ofboth unsigned normal and tangential forces of the typicaladhesive packings that listed in Table III. The plateau atthe biggest forces and the deviation of the smallest forcesare believed to be caused by the finite size effect, wherethe probability of finding the biggest or smallest forceshould be around 1 /N Z ( 10 − × − P ( f ) ∼ f θ for small forces as well FIG. 7: (Colors online) Signed normal force network anddistribution of typical packings of adhesive particles. Neg-ative values are repulsive while positive are attractive. Allthe forces are normalized with the mean value of the magni-tude of normal forces. Case (a)(c)(e)(g) are consistent withthe cases listed in Table III. as the power law P ( f ) ∼ exp − βf for big forces [48–50].However, the exponents are θ = 0 . β = 0 .
839 fornormal forces while θ = 0 . β = 0 .
888 for tangentialforces, respectively. These exponents are distinct fromthe reported values of θ ≈ . ∼ . β ≈ . ∼ . q ij of load oneach particle [51]. With the sum of the fraction equalingto one, the sum of all the force components in the verticaldirection is conserved. Nevertheless, this q model is stilltoo simple to account for the force behaviour of adhesivepackings for two reasons. One is that only the verticalforces are considered. The other is that the existence ofadhesion weakens the load that the underneath particlescarry from the upper ones. Thus it might lead the frac-tion q ij to be negative, where the forces between two con-tact particles are attractive. But we still speculate thatthe adhesion will not make much effect on the distribu-tion of big forces. On the other hand, for small forces, itis reported that they have a close relation with the stabil-ity of packing. Contacts carrying small forces which arerepulsive for granular matters are easy to break. How-ever, for adhesive particles, the process of necking willtake place and a critical pull-off force F C needs to bereached to separate two contact particles. Thus contactswith small forces in adhesive packings are relatively morestable, no matter attractive or repulsive. Regarding theexponents, we conjecture that the very small and verybig repulsive forces in granular packing will not likely ex-ist in the presence of adhesion. The very small repulsiveforces could be enhanced while the very big ones will beweakened by attractive adhesion force, which will resultin the shrink of the force distribution. Therefore the ex-ponent θ of small forces will increase a little and β of bigforces will decrease. However, this is just a qualitativeanalysis where there have been few data or theory to sup- FIG. 8: (Colors online) Unsigned normal (a) and tangential(b) force distribution of typical packings of adhesive particles.Case (a) (g) in the legend are consistent with the cases listedin Table III.FIG. 9: (Colors online) Schematic diagrams illustrating thestress distributions and force balance of two contact particlesfor Hertz (a) and JKR (b) contact models. port. More insight from both experiment and theoreticalefforts is left to future work.For granular matters without adhesion, at least threeparticles (in 3D) are usually needed to support one andmake it mechanical equilibrium. However, for the veryloose packings of adhesive particles, many particles onlyhave one or two neighbours which are not able to achievemechanical equilibrium in terms of granular matters. In-deed, these particles can still be mechanically stabilizedfor two reasons. The first is the adhesive attractive forcecaused by van der Waals interactions that are especiallyprominent for micron-sized particles. The other one isfriction as well as the negligible gravity compared withvan der Waals force. Next we will conduct the mechan-ical equilibrium analysis and explain qualitatively whyadhesive particles can be stabilized more easily.For convenience and better illustration, two contactparticles p p θ is the angle betweenthe contact surface and the horizontal and a is the radiusof the contact surface. Fig. 9 also illustrates the stress distributions of Hertz and JKR contact models, which areused for adhesive-less granular matter and our adhesivepackings, respectively. For both models, the stress dis-tribution is symmetric. However, stress in Hertz modelis totally repulsive while that in JKR model is repul-sive in the central part and attractive on the rim of thecontact surface. The unknown forces that act on p f t and two normal forces f n , f n simplified by the stress distribution. When the angle θ and the external force F ext are known, two equations offorce balance and one equation of torque balance can besolved to figure out the three unknown forces, indicatingan isostatic condition. Here F ext could be gravity, elec-tric force or other kinds of external forces that act on theparticle. Usually, the particle p p θ and F ext reach some critical values, whichis the cause of rearrangements in a packing. However,due to the presence of adhesive forces, the critical slidingforce is greater than that of non-adhesive particles (seeEq. 9), meaning that particles with adhesion are moredifficult to start sliding. For rolling, on the other hand,when two contact particles start to roll or have the ten-dency of rolling, the front side of the contact surface iscompressed and the rear side will still be in touch untilthe critical pull-off force F C is reached, instead of de-tach immediately when there is no external force. As aconsequence, the simplified normal forces f n , f n willconsist of a repulsive force on the front side of rollingand an attractive force on the rear side, while the twonormal forces are both repulsive for adhesive-less par-ticles (see the normal forces in Fig. 9 indicated by thered arrows). The two normal forces of adhesive parti-cles will generate torque of the same direction, resultingin greater rolling resistance. In this case, particles withadhesion can support larger external forces as well as bestabilized within a wider range of θ with constant exter-nal forces. Based on the above analysis, an equilibriumdiagram is plotted with the parameters used in our workto demonstrate the equilibrium region in terms of angle θ and external force F ext . Fig. 10 shows two series ofequilibrium lines, the rolling equilibrium (re.) line andthe sliding equilibrium (se.) line, under which the areaindicates the equilibrium region that the particle will notroll or slide over another. We can see that with verylow external force, the adhesive particles can even bestabilized with θ = 90 ◦ . Then when the external forcegrows larger, the angle θ decreases fast and reaches about θ ≈ ◦ when the force further grows, implying that par-ticles can only be stabilized when they are placed almostvertically. Also we found that the rolling equilibrium linelies under the sliding equilibrium line, shifting to left withincreasing particle size, which means that rolling equi-librium breaks ahead of sliding and the larger size theparticle is, the more easily it will roll. This is in agree-ment with the statement that sliding is relatively rarefor small adhesive particles while rolling is generally the0 FIG. 10: (Colors online) Equilibrium diagram of two contactparticles with JKR contact model. The solid line stands forthe sliding equilibrium (se.) line. The dash, dot and shortdash lines are the rolling equilibrium lines (re.) for the par-ticles in our simulations with r p = 1 µm , r p = 10 µm and r p = 50 µm respectively. The short dot line represents theequilibrium line of Hertz contact model. The shady area in-dicates the equilibrium region of r p = 50 µm . preferred deformation mode for agglomerates of adhesiveparticles [18]. For comparison, the equilibrium line ofHertz model is also plotted in Fig. 10 as indicated by theshort dot line. The particle without adhesion can onlyreach equilibrium in a tiny range of θ ( < . ◦ ). This isbecause in Hertz model the normal force f n on the rearside of the contact surface must be repulsive and thuscannot provide additional rolling resistance like that inJKR model, leading to probable rolling motion. After all,from the analysis above, we conclude that particles areable to get mechanical equilibrium with fewer neighborsthan adhesive-less granular matter due to the existenceof adhesion. CONCLUSIONS
In this paper, random loose packings of uniform spheri-cal micro-particles are investigated by using particle-levelDEM simulations on the basis of adhesive contact me-chanics. The packing structures arise from the compe-tition between particle inertia and adhesion. A charac-teristic adhesion parameter ( Ad ) with respect to particleinertia is used to account for the combined effect of par-ticle velocity, size and adhesion. When 1 < Ad <
20, thevolume fraction and coordination number are uniquelydependent on Ad , resulting in a new regime of adhesiveloose packing. As Ad grows above 20, the packing prop-erties approach φ = 0 . Z = 2 and no longer changeswith Ad , confirming the asymptotic adhesive loose pack-ing limit. On the other hand, with Ad <
1, indicatinga weak adhesion, the packing properties go back to the range between RLP and RCP. A simple form of equa-tion φ RCP − φ = αAd λ is figured out to interpret therelationship between volume fraction and adhesion pa-rameter, where the lowest volume fraction is achievedwith Ad >
20. Furthermore, the force distribution of ad-hesive packings resembles that of adhesive-less ones with P ( f ) ∼ f θ for small forces and P ( f ) ∼ exp − βf for bigforces. The distinct exponents of θ = 0 . β = 0 . ACKNOWLEDGEMENTS
This work has been funded by the National NaturalScience Funds of China (Nos. 50976058 and 51390491)and the National Key Basic Research Program of China(No. 2013CB228506). S. Q. Li is grateful to Prof. JeffMarshall at University of Vermont, Prof. Hernan Makseat City College of New York and Dr. Guanqing Liu andDr. Mengmeng Yang at Tsinghua University for helpfuldiscussions. W. Liu acknowledges Shaojun Luo at CityCollege of New York for his support of the simulationwork. [1] J. D. Bernal, Nature , 141 (1959), URL http://dx.doi.org/10.1038/183141a0 .[2] G. Parisi and F. Zamponi, Rev. Mod. Phys. ,789 (2010), URL http://link.aps.org/doi/10.1103/RevModPhys.82.789 .[3] A. Coniglio, A. Fierro, H. J. Herrmann, andM. Nicodemi, Unifying Concepts in Granular Media andGlasses: From the Statistical Mechanics of Granular Me-dia to the Theory of Jamming (Elsevier, 2004).[4] B. Andreotti, Y. Forterre, and O. Pouliquen,
Granularmedia: between fluid and solid (Cambridge UniversityPress, 2013).[5] G. D. Scott, Nature , 908 (1960).[6] G. Y. Onoda and E. G. Liniger, Phys. Rev. Lett. ,2727 (1990), URL http://link.aps.org/doi/10.1103/PhysRevLett.64.2727 .[7] M. P. Ciamarra and A. Coniglio, Phys. Rev. Lett. , 128001 (2008), URL http://link.aps.org/doi/10.1103/PhysRevLett.101.128001 .[8] M. Jerkins, M. Schr¨oter, H. L. Swinney, T. J. Senden,M. Saadatfar, and T. Aste, Phys. Rev. Lett. ,018301 (2008), URL http://link.aps.org/doi/10.1103/PhysRevLett.101.018301 . [9] K. J. Dong, R. Y. Yang, R. P. Zou, and A. B. Yu, Phys.Rev. Lett. , 145505 (2006), URL http://link.aps.org/doi/10.1103/PhysRevLett.96.145505 .[10] C. Song, P. Wang, and H. A. Makse, Nature , 629(2008), URL .[11] G. R. Farrell, K. M. Martini, and N. Menon, Soft Matter , 2925 (2010).[12] A. V. Anikeenko and N. N. Medvedev, Phys. Rev. Lett. , 235504 (2007), URL http://link.aps.org/doi/10.1103/PhysRevLett.98.235504 .[13] B. A. Klumov, S. A. Khrapak, and G. E. Morfill, Phys.Rev. B , 184105 (2011), URL http://link.aps.org/doi/10.1103/PhysRevB.83.184105 .[14] Y. Jin and H. A. Makse, Physica A: Statistical Mechan-ics and its Applications , 5362 (2010), ISSN 0378-4371, URL .[15] A. Panaitescu, K. A. Reddy, and A. Kudrolli, Phys. Rev.Lett. , 108001 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.108.108001 .[16] S. Torquato and F. H. Stillinger, Rev. Mod. Phys. ,2633 (2010), URL http://link.aps.org/doi/10.1103/RevModPhys.82.2633 .[17] S. Li, J. S. Marshall, G. Liu, and Q. Yao,Progress in Energy and Combustion Science , 633(2011), URL .[18] J. S. Marshall and S. Li, Adhesive Particle Flow (Cam-bridge University Press, 2014).[19] C. Dominik and A. G. G. M. Tielens, The AstrophysicalJournal , 647 (1997), URL http://stacks.iop.org/0004-637X/480/i=2/a=647 .[20] J. Blum and R. Schr¨apler, Phys. Rev. Lett. ,115503 (2004), URL http://link.aps.org/doi/10.1103/PhysRevLett.93.115503 .[21] K. M. Kinch, J. Sohl-Dickstein, J. F. Bell, J. R. John-son, W. Goetz, and G. A. Landis, Journal of GeophysicalResearch: Planets , E06S03 (2007), ISSN 2156-2202,URL http://dx.doi.org/10.1029/2006JE002807 .[22] R. Y. Yang, R. P. Zou, and A. B. Yu, Phys. Rev. E ,3900 (2000), URL http://link.aps.org/doi/10.1103/PhysRevE.62.3900 .[23] J. M. Valverde, M. A. S. Quintanilla, and A. Castellanos,Phys. Rev. Lett. , 258303 (2004), URL http://link.aps.org/doi/10.1103/PhysRevLett.92.258303 .[24] J. Blum, R. Schr¨apler, B. J. R. Davidsson, and J. M.Trigo-Rodr´ıguez, The Astrophysical Journal , 1768(2006), URL http://stacks.iop.org/0004-637X/652/i=2/a=1768 .[25] C. L. Martin and R. K. Bordia, Phys. Rev. E , 031307 (2008), URL http://link.aps.org/doi/10.1103/PhysRevE.77.031307 .[26] G. Lois, J. Blawzdziewicz, and C. S. O’Hern, Phys. Rev.Lett. , 028001 (2008), URL http://link.aps.org/doi/10.1103/PhysRevLett.100.028001 .[27] E. J. R. Parteli, J. Schmidt, C. Blumel, K.-E. Wirth,W. Peukert, and T. Poschel, Sci. Rep. , 06227 (2014),URL http://dx.doi.org/10.1038/srep06227 .[28] N. V. Brilliantov, F. Spahn, J.-M. Hertzsch, andT. P¨oschel, Phys. Rev. E , 5382 (1996), URL http://link.aps.org/doi/10.1103/PhysRevE.53.5382 .[29] N. V. Brilliantov, N. Albers, F. Spahn, and T. P¨oschel,Phys. Rev. E , 051302 (2007), URL http://link.aps. org/doi/10.1103/PhysRevE.76.051302 .[30] S. Luding and H. J. Herrmann, Zur Beschreibung kom-plexen Materialverhaltens, Institut f¨ur Mechanik, S.Diebels (Ed.), Stuttgart pp. 121–133 (2001).[31] K. L. Johnson, K. Kendall, and A. D. Roberts, Proceed-ings of the Royal Society of London A: Mathematical,Physical and Engineering Sciences , 301 (1971), ISSN0080-4630.[32] B. Derjaguin, V. Muller, and Y. Toporov, Progressin Surface Science , 131 (1994), ISSN 0079-6816, URL .[33] G. Liu, S. Li, and Q. Yao, Frontiers of Energyand Power Engineering in China , 280 (2010),ISSN 1673-7393, URL http://dx.doi.org/10.1007/s11708-009-0062-5 .[34] S.-Q. Li and J. Marshall, Journal of Aerosol Sci-ence , 1031 (2007), ISSN 0021-8502, URL .[35] J. Marshall, Journal of Computational Physics , 1541 (2009), ISSN 0021-9991, URL .[36] M. Yang, S. Li, and Q. Yao, Powder Technology ,44 (2013), ISSN 0032-5910, discrete Element Mod-elling, URL .[37] L.-O. Heim, J. Blum, M. Preuss, and H.-J. Butt, Phys.Rev. Lett. , 3328 (1999), URL http://link.aps.org/doi/10.1103/PhysRevLett.83.3328 .[38] K. Bagi and M. R. Kuhn, Journal of applied mechanics , 493 (2004).[39] C.Dominik and A. G. G. M. Tielens,Philosophical Magazine A , 783 (1995),http://dx.doi.org/10.1080/01418619508243800, URL http://dx.doi.org/10.1080/01418619508243800 .[40] B. S¨umer and M. Sitti, Journal of Adhe-sion Science and Technology , 481 (2008),http://dx.doi.org/10.1163/156856108X295527, URL http://dx.doi.org/10.1163/156856108X295527 .[41] P. A. Cundall and O. D. Strack, Geotechnique , 47(1979).[42] C. Thornton, Journal of Physics D: Applied Physics ,1942 (1991), URL http://stacks.iop.org/0022-3727/24/i=11/a=007 .[43] C. Thornton and K. Yin, Powder Technol-ogy , 153 (1991), ISSN 0032-5910, URL .[44] R. Jones, H. M. Pollock, D. Geldart, and A. Verlinden-Luts, Ultramicroscopy , 59 (2004), ISSN 0304-3991, URL .[45] W. Liu, S. Li, A. Baule, and H. A. Makse, Soft Mat-ter , 6492 (2015), URL http://dx.doi.org/10.1039/C5SM01169H .[46] S. Torquato and F. H. Stillinger, Ex-perimental Mathematics , 307 (2006),http://dx.doi.org/10.1080/10586458.2006.10128964,URL http://dx.doi.org/10.1080/10586458.2006.10128964 .[47] A. Singh, V. Magnanimo, K. Saitoh, and S. Luding, Phys.Rev. E , 022202 (2014), URL http://link.aps.org/ doi/10.1103/PhysRevE.90.022202 .[48] D. M. Mueth, H. M. Jaeger, and S. R. Nagel, Phys. Rev.E , 3164 (1998), URL http://link.aps.org/doi/10.1103/PhysRevE.57.3164 .[49] E. Lerner, G. During, and M. Wyart, Soft Matter , 8252(2013), URL http://dx.doi.org/10.1039/C3SM50515D .[50] L. Bo, R. Mari, C. Song, and H. A. Makse, Soft matter , 7379 (2014).[51] C. h. Liu, S. R. Nagel, D. A. Schecter,S. N. Coppersmith, S. Majumdar, O. Narayan,and T. A. Witten, Science