Force networks and elasticity in granular silos ∗ John F. Wambaugh , Robert R. Hartley, and Robert P. Behringer Department of Physics and Center for Nonlinear and Complex Systems, Duke University, Durham, NC 27708† (Dated: February 2, 2008) Wehavemadeexperimentalobservationsoftheforcenetworkswithinatwo-dimensionalgranular silosimilartotheclassicalsystemofJanssen. ModelslikethatofJanssenpredictthatpressurewithin a silo saturates with depth as the result of vertical forces being redirected to the walls of the silo wheretheycanthenbecarriedbyfriction. Byaveragingensemblesofexperimentally-obtainedforce networksin differentways, wecomparetheobserved behaviorwith variouspredictionsfor granular silos. We identify several differences between the mean behavior in our system and that predicted 8 by Janssen-like models: We find that the redirection parameter describing how the force network 0 transfers vertical forces to the walls varies with depth. We find that changes in the preparation 0 of the material can cause the pressure within the silo to either saturate or to continue building 2 with depth. Most strikingly, we observe a non-linear response to overloads applied to the top of n the material in the silo. For larger overloads we observe the previously reported “giant overshoot” a effect where overload pressure decaysonly after an initial increase [G. Ovarlezet al., Phys. Rev. E J 67, 060302(R) (2003)]. For smaller overloads we find that additional pressure propagates to great 2 depth. Thiseffectdependsontheparticlestiffness,asgivenforinstancebytheYoung’smodulus,E, 2 of the material from which the particles are made. Important measures include E, the unscreened hydrostatic pressure, and the applied load. These experiments suggest that when the load and ] the particle weight are comparable, particle elasticity acts to stabilize the force network, allowing t f non-linearnetwork effects to beseen in themean behavior. o s PACSnumbers: 45.70.-n,45.70.Cc,83.80.Fg,45.05.+x . t a m I. INTRODUCTION to be in these positions [1]. Determining a statisticalde- - scriptionofthesecontactsisatthecoreofunderstanding d dense granular matter. n Thephysicsofgranularmatterisofgreatinterestboth o Grain-scale statistical descriptions are further compli- because of the enormous breadth of physical systems in c cated by the presence of force networks [2, 3]. An ex- the granular regime as well as the lack of a characteri- [ traordinary feature of granular contact networks is that zation of the granular state in fundamental terms. Al- the distributionofforcesonthe contactscanbe strongly 1 though there has been some success understanding ener- v inhomogeneous [4]. Certain sequences of contacts carry getic granular gases in terms of kinematic models, dense 7 forceswithmagnitudemanytimesthemeanforlongdis- granular materials continue to pose difficult problems. 8 tances along “force chains”. As shown in Fig. 1, force 3 We consider dry granular matter at low humidity where networkscan developto support the weight of the parti- 3 there are no attractive forces between the constituent cles as well as any additional applied load. . particles and the particles do not significantly interact 1 Despitethedifficultiesinherenttounderstandinggran- with the interstitial fluid (air). Contact forces dominate 0 ular materials, phenomenologicalengineering models ex- 8 dense granular materials in this regime, with each parti- istthat describe the mean behaviorwellenoughto allow 0 cle having several persistent contacts. for the design of granular storage facilities and handling : v Coulomb’s law of friction is an inequality relating the processes. These models can provide useful tools for in- Xi magnitude of the friction force, Ff to the normal force, vestigating granular materials. In particular, Janssen’s Fn, depending upon the mobilization of the friction. If model of granular materials in vertical silos is interest- r a µs is the static friction coefficient, then the mobilization ing because it has been thoroughly tested, is realizable is given by Ff/(Fnµs). For a non-moving contact the inthe laboratory,andultimately is viewedto be a work- mobilization ranges from −1 to 1 and reflects both the ing qualitative and even quantitative description of the direction in which the force of friction acts to resist mo- mean behavior of granular matter [5, 6]. tion and how close the forces acting on the contact are to causingmotionin the opposite direction. Becausethe forceoffrictionisdependentuponhowthecontactarose, II. THE JANSSEN MODEL the granularstate is determined bothby the positionsof the particles and the history of how the particles came Janssen treats the material in a silo as a continuum and considers horizontal slices, for which the difference of stress, the weight of the material in the slice, and the force of friction at the wall must be balanced ver- ∗Currentaddress: National Center forComputational Toxicology, USEPA,ResearchTrianglePark,NC27711 tically. Here we consider a quasi-two-dimensional silo, †Electronicaddress: wambaugh@phy.duke.edu such as the one we use in the experiments reported be- 2 1 0.8 e ur 0.6 s s 0.4 e r P 0.2 0.5 1 1.5 2 2.5 3 Depth Harbitrary unitsL FIG. 2: Janssen’s model predicts that stress exponentially approachesasaturationvaluewithdepthinasilo. Overloads decay exponentially and if an overload equalto thethesatu- ration stress is placed on thesilo thestress is constant. Here we plot predictions in a system that saturates at 0.5 in arbi- traryunitsfornoload, asaturation overload (=0.5), andan abovesaturation overload. |σ |≤ µ |σ | for wall friction µ . Janssen further as- xy w xx w sumes that the vertical and horizontal normal stresses are linearly related to each other to write σ =kµ σ xy w yy where k represents both the mobilization of interparti- cle contact friction (with friction coefficient δ) and the mobilization of friction at the walls. By considering the extremes of maximum upward and downward mobilization, it has been shown that 1−δ −1 P10 ) < m c 40 >)/ ( λ P n 30 < e − −2 s P10 s n 20 (( a b J o r 10 P 0 −2 0 2 4 0 0.5 1 1.5 2 Depth from Surface (L) (P−

FIG. 12: Due to the insensitivity of the harder photoelas- ticparticles, determiningthelocalJanssen screeninglength- ) scaleλnumerically isonlypossibleforthelargest 3402 g(N) > −1 P10 overload. The observed screening length is roughly half as < long as for thesoft (E =4 MPa) particles. / ) > P < to E = 210 MPa. For the harder photoelastic material, − −2 P10 the sensitivity to applied stress was lower, with the re- ( sult that the stress due to gravity was not visible even b( atgreatdepth in the silo. Inthis case,any pressure pro- o r file is effectively an overload decay profile, since there is P no unloaded profile to subtract. As shown in Fig. 11, the additional pressure due to the overloadsdecays with −2 0 2 4 depth after overshootingfor allthree cases. For allthree (P−

overloads,the mean pressure eventually drops below the range of observable sensitivity, though large fluctuations FIG. 13: Probability distributions for the difference in pres- aboutthe meanwerestillpresent. Allthree casescanbe sure from the mean value at the center (solid) and sides fit with exponentials if the initial overshootis neglected. (open) of thesilo forsoft (top) particles in an unloaded silo Since the non-decaying, deep-penetrating overloads for (•)andwithoverloadsof20g((cid:4))and80g(N). Wealsoshow the softer particles were close to the saturationpressure, distributions for hard (bottom) particles with overloads of itisnotpossibletodeterminewhetherornotsucheffects 1837 g (•), 2267 g ((cid:4)) and 3402 g (N). The curves differ for are present in the harder particles using our method. low pressures duetothe differences in sensitivity of thehard Determining the Janssen screening length-scale λ nu- andsoft photoelastic polymers, butgenerally thereisagreat deal of overlap and we cannot distinguish the sides from the mericallyforthehardparticleoverloadswasproblematic center. since the pressure profiles eventually decayed below the sensitivity of the hard photoelastic material. In Fig. 12, we show the calculated λ(z) for the largest, 3402 g overload in the top third of the silo. What we see is qual- itatively similar to the behavior of the largest overload in Fig. 10 — λ approaches a constant value from above. maximumload onthe harderparticles, ∼3400g,this ra- Surprisingly, although the system geometry is the same, tio is ∼1.0×10−4. Thus, in appropriately scaled terms thescreeninglengthisroughlyhalfofwhatwasfoundfor the loads are comparable in both cases, and hence on the soft (E = 4 MPa) particles used to obtain the first this basis, elasticity cannot explain the difference in the two pressure profiles (Fig. 5 and Fig. 9). experimentalobservations. However,the ratioof the ap- The two sets of experiments nominally contrast soft plied load to the weight of the material is significantly and hard particles. However, there is another aspect different in the two cases, namely about 0.5 for the soft which may be at least as important as granular elastic- particles, and 17 for the hard particles. Thus, we con- ity. We note that the ratio of the applied loads to the clude thatforthe softerparticles,the loadsarenotlarge Young’s modulus, E, are comparable for the hard and enough to overcome the effects caused by gravitational the soft material. For the maximum load on the softer loading, which necessarily reflects the preparation his- particles, ∼ 100g, this ratio is ∼ 1.6×10−4 and for the tory. 9 C. Horizontal Variations and Fluctuations 12000 Itisalsointerestingtoinvestigatevariousotheraspects 10000 of the force profiles which differ from the usual Janssen m) picture. In particular, we consider horizontal variabil- c ity of the forces, and also force fluctuations. In this re- e/ 8000 n gard, we note model studies by Pitman which exhibited y d horizontal dependence on fluctuations and stress within ( 6000 e granular silos[16] due to inhomogeneous mobilization of ur friction. ss 4000 e As shown in Fig. 13, we measured the local deviation r P from the mean pressure at a given height for both the 2000 centerandsidesofthesiloforalloverloadsforbothhard and soft particles. We did not see significant horizontal 0 1 2 3 4 5 6 variation in the distributions, and for the soft particles Depth from Surface (L) the distributions were nearly identical for all cases and 7 veryreminiscentofforcedistributions observedforother x 10 2 granular systems [30]. Forces below the mean fall off quickly, while forces larger than the mean tail off slowly, and roughly follow an exponential. m) 1.5 Recent simulations for vertical granular matter with c / e periodic boundary conditions (in the horizontal direc- n y tions)haveindicatedthatthe averaginglength-scalesuf- d ( 1 ficienttodeterminethemeanpressureispossiblyassmall e as a grain diameter [31]. Though our geometry has side ur s walls, it is still worthwhile to determine the length over s which we must average for the pressure profile to con- re0.5 P verge. To do this, we determine a grid of points within the silo andaveragethe pressurewithin atunable radius toassignapressuretothatpoint. InFig.14,weplotthe 0 1 2 3 4 5 6 pressureinthe middle of the silo (onaxis)withdepth as Depth from Surface (L) determined by averaging over regions with radii of 1, 5, 10 and 15 small particle diameters. For the two larger FIG. 14: Pressure profiles on axis converge as the length R averagingregions,theprofilesconvergetothoseobserved overwhichaveragingisconductedvaries,asindicatedbylines using horizontal binning by depth in Fig. 9. This con- of increasing thickness for R=1, 5, 10 and 15 small particle vergence indicates a correlation length-scale between 10 diametersforbothsoftparticles(top)inanunloadedsilo(•) and 15 particle diameters — approximately the same as and with overloads of 20 g ((cid:4)) and 80 g (N) and for hard the Janssen screening length-scale obtained numerically. particles (bottom) with overloads of 1837 g (•), 2267 g ((cid:4)) Thisisinterestingsinceinthisinstance,themacroscopic and 3402 g (N). scale of the Janssen parameter λ is comparable to the averagingscale,eventhoughtheJanssenmodeldoesnot include particle size. As can be seen for the lower image between two particles must be equal to the sum of their in Fig. 14, the the effect of averaging length is roughly radii. Thisallowsustodeterminealocalpackingfraction the same regardless of particle stiffness. throughout the silo. Anadditionalaspectofinterestisthegeometricnature Since previous studies have shown that the packing of the contacts. Relevant quantities include the packing fraction of a pile relates to the redirection parameter, fraction, φ, the number of contacts per particle, Z, and weinvestigatedepth-dependentvariationsofthepacking clustering coefficient, defined below. fraction [8]. Due to the dilation of the material, we find With our second preparation procedure, we are able thatthe packingfractionis slightlylowerforloadedsilos to obtain to obtain all of these quantities. We do this for the depths availablein the unloadedcase. We do not by capturing both polarized images showing the force find systematic correlation between the packing fraction network within the photoelastic disks and unpolarized and the screening length-scale. images showing the locations of those particles. Particle Theparticlelocationsandcontactscompriseacontact centers can be identified using a convolution with sep- network for each pile. Analyzing these contact networks arate kernels for the two particle sizes. We then use a allowsthedeterminationofgraphandnetworkproperties numerical algorithm to reject overlapping particles and such as the average number of contacts per particle, Z. identify neighbor particles in contact with each particle We can also determine the ‘clustering coefficient’, which baseduponthegeometricconsiderationthatthedistance gives a measure of how tightly interconnected is the net- 10 Contact Network Force Network 1 5 5 4 4 0.9 φ n Z 3 3 S o o acti 0.8 cle 2 2 ft g Fr 0.7 Parti 1 1 kin er 00 2 4 6 00 2 4 6 c a P P 0.6 s ct 5 5 a 0.5 nt 4 4 0 2 4 6 o C 3 3 Depth (L) H n 1 Mea 2 2 ard 1 1 0.9 0 0 φ 0 2 4 6 0 2 4 6 n Depth (L) o cti 0.8 a Fr FIG. 16: The average number of contacts per particle is g 0.7 roughly constant for the contact network (left) but varies n ki at the top of the pile for the force network (right). This is c true for both soft particles (above) in an unloaded silo (•) a P 0.6 and with overloads of 20 g ((cid:4)) and 80 g (N) as a well as for hard particles (below) with overloads of 1837 g (•), 2267 g ((cid:4)) and 3402 g (N). 0.5 0 2 4 6 Depth (L) material, although slightly lower due to the presence of FIG. 15: The average packing fraction for soft (left) and the walls [34]. The average number of contacts in the hard (right) particles varies slightly depending upon over- force network is slightly less than to two, which is not load. Whileitisroughlyconstantwithdepthforanunloaded surprising since many particles in a chain have only two silo (•) it varies for overloads of 20 g ((cid:4)) and 80 g (N). contacts and, in our relative small silo, many particles make contact with the walls. Significantly, for the force network within both the work [32]. For each particle we consider the n neighbor- hard and soft particles, the average number of contacts ingparticleswithwhichthatparticleisincontact. There at the top of the pile varies depending upon the over- are n−1 possible neighbor-neighbor contacts for round load, before converging with depth, unlike the contact particles in two dimensions, so we define the clustering network which does not vary with overload. This indi- coefficient as the ratio of the m actual contacts present cates that the force network is quantitatively different, between those neighbors to the n−1 possible contacts, depending uponthe size ofthe appliedoverload. For the C ≡ m . deep-penetratingoverloadonthesoftparticles,themean n−1 For each contact network, we can determine a subset numberofcontactsatthe topisclosesttothesaturation ofcontactsforwhichtheparticle-scalepressureisgreater value, perhaps indicating that the structure of the force thanthemeanforthewholesilo. Thissubsetofcontacts network is made visible by the additional pressure, but canbethoughtofasroughlycomprisingaforcenetwork, is not disturbed as in the case of larger overloads. for which the same graph and network properties as the Fortheclusteringcoefficient,Fig.17showsthecontact contact network can also be calculated [33]. network to be better connected than the force network. We first examine the dependence of the mean number Aswiththenumberofcontacts,theclusteringcoefficient of contacts per particle with depth for both the over- is roughly constant for the contact network but varies all contact network and the force network. As shown in with depth for the force network. For hardparticles, the Fig.16,wefindthatforthecontactnetwork,thenumber overloads do not seem to alter the clustering coefficient of contacts is nearly independent of overload and con- of the contact or force network. stant with depth at a value close to k = 4, the typically Thedifferencesbetweenthecontactandforcenetworks observed coordination number for a disordered granular may indicate an explanation for the deep propagation

Force networks and elasticity in granular silos

Similar Force networks and elasticity in granular silos