text:
# Statistics of the contact network in frictional and frictionless granular packings

## I Introduction

The properties of granular materials, even static sandpiles, continue to perplex engineers and physicists alike . Processing of granular materials play a central role in the pharmaceuticals industry as well as engineering communities, such as ceramic component design. For example, one may wish to evenly distribute the ingredients in a tablet or pill, or reduce the likelihood of component failure. However as yet there is no clear indication of how the individual particle properties determine the final state of the system.

Although packings of frictionless, monodisperse, cohesionless, hard spheres have been well studied , little is known about the effect of including particle friction. Recent discrete element simulations of granular materials, where packings were generated for particles for different static coefficients of friction $`\mu `$ between individual particle pairs, showed that the local particle coordination of the packing varied strongly as a function of friction . From a different perspective, experimental studies of static granular assemblies have shown many interesting facets of the stress state of these systems. One method of analysis appears to dominate in describing the statistics of granular packings: computations of the distributions of normal forces are de rigeur.

Experimental studies on granular packings use a variety of techniques to measure the distribution of normal contact forces $`P(f)`$ between particles and container walls , where $`fF/\overline{F}`$ \- all measured normal forces $`F`$ are normalized with respect to the average force $`\overline{F}`$. The Chicago group  utilized carbon paper to measure $`P(f)`$ at the base and sides of a cylindrical container packed with glass spheres with a normal load applied at the top of the packing. Forces several times the average force were observed, with resolution down to the weight of a few particles. Blair et al.  measured $`P(f)`$ for amorphous and ordered granular packings for particles with different values of $`\mu `$, which varied by a factor of approximately three. These experiments demonstrated that $`P(f)`$ is indiscriminate towards the effects of particle friction and structure of the packing, and the general form of $`P(f)`$ remained robust within the resolution of the experiment.

In a different experimental set up, L$`\varphi `$voll et al.  used a pressure transducer device to measure $`P(f)`$ at the bottom of an unloaded granular packing under its own weight, on a fixed substrate of particles glued to the supporting base. This experiment was able to resolve forces down to the weight of a few grain masses and showed that the spatial distribution of contact forces were correlated over a few particle diameters. Using a novel modification of the carbon paper technique, Tsoungui et al.  actually measured $`P(f)`$ inside the bulk of a 2D packing. Despite the poorer statistics of this study, the results agreed well with Blair et al.  and L$`\varphi `$voll et al. . Experimental studies on static granular packings show that $`P(f)`$ exhibits several generic features; an approximately exponential tail at large $`f`$ and a plateau or peak near $`f1`$.

Computational studies, such as contact dynamics or molecular dynamics, of compressed packings provide good comparison with the experimental data . However, as yet, there has been no systematic study of the effects of particle friction on the force distributions within a granular assembly. Here we show how the effects of friction change the behaviour of $`P(f)`$ in the small force region but only weakly affect the large-$`f`$ region. We show that the local contact geometry of the packing is not isotropic and only weakly influenced by friction. We also discuss aspects of the force network whereby high force-bearing structures, or force chains, do not seem to be a dominant feature of these unloaded packings.

We computed $`P(f)`$ in the bulk of various packings (which is presently inaccessible in 3D experiments) that had settled onto either a rough bed or a planar base. We compared these results with $`P(f)`$ for particles in contact with the flat base (similar to experiment) of a periodic packing and with the $`P(f)`$ generated at the side walls of a cylindrical packing. We resolve the components of the contact force that are normal ($`n`$) and tangential ($`t`$) to the line of centres between two particles in contact.

In the next section we briefly describe the model, though a more thorough description of the technique is available elsewhere . In section III, we present results for the force distributions, force correlations, and the contact geometry. We also discuss some aspects of the force network with respect to a force cut-off scheme, highlighting some pros and cons of this method. In section IV we summarize and conclude this work.

## II Model and Method

We performed three dimensional (3D) molecular dynamics simulations with $`N`$ monodisperse, cohesionless, inelastic spheres that interact only on contact via a Hooke (linear) spring or a Hertz contact law and static friction . Contacting particles $`i`$ and $`j`$ positioned at $`𝐫_i`$ and $`𝐫_j`$ experience a relative normal compression $`\delta =|𝐫_{ij}d|`$, where $`𝐫_{ij}=𝐫_i𝐫_j`$, which results in a force $`𝐅_{ij}=𝐅_n+𝐅_t`$. The normal and tangential contact forces are given by,

$$𝐅_n=f(\delta /d)\left(k_n\delta 𝐧_{ij}\frac{m}{2}\gamma _n𝐯_n\right),$$
(1)
$$𝐅_t=f(\delta /d)\left(k_t\mathrm{\Delta }𝐬_t\frac{m}{2}\gamma _t𝐯_t\right),$$
(2)
where $`𝐧_{ij}=𝐫_{ij}/r_{ij}`$, with $`r_{ij}=|𝐫_{ij}|`$, $`𝐯_n`$ and $`𝐯_t`$ are the normal and tangential components of the relative surface velocity, and $`k_{n,t}`$ and $`\gamma _{n,t}`$ are elastic and viscoelastic constants respectively. $`f(x)=1`$ for Hookean springs and $`f(x)=\sqrt{x}`$ for Hertzian contacts. $`\mathrm{\Delta }𝐬_t`$ is the elastic tangential displacement between spheres, obtained by integrating surface relative velocities during elastic deformation of the contact. The magnitude of $`\mathrm{\Delta }𝐬_t`$ is truncated as necessary to satisfy a local Coulomb yield criterion $`F_t\mu F_n`$, where $`F_t|𝐅_t|`$ and $`F_n|𝐅_n|`$, and $`\mu `$ is the particle-particle friction coefficient. For the present simulations we set $`k_n=210^5mg/d`$, $`k_t=\frac{2}{7}k_n`$, $`\gamma _n=50\sqrt{g/d}`$. For Hookean springs we set $`\gamma _t=0`$ while for Hertzian springs, $`\gamma _t=\gamma _n`$. For Hookean springs the coefficient of restitution $`ϵ_{n,t}`$, is related to $`\gamma _{n,t}`$ through,

$`ϵ_{n,t}=\mathrm{exp}(\gamma _{n,t}t_{\text{c}ol}/2),`$

where the collision time $`t_{\text{c}ol}`$ is determined by the contact frequency between two particles. For the parameters chosen, $`ϵ_n=0.88`$ for Hookean springs. For Hertzian contacts $`ϵ`$ is velocity dependent. We chose a time-step $`\delta t=10^4\tau `$, where $`\tau =\sqrt{d/g}`$ and $`g`$ is the gravitational acceleration.

Amorphous packings (with packing fraction $`\varphi 0.60`$) were generated by allowing an initially dilute system to settle under gravity acting in the vertical direction. Particles settled onto a bottom wall that was either a planar base or a bumpy bed of particles frozen into a close packed random configuration. This process was run until the kinetic energy of the system was much smaller than the potential energy . The base had the same frictional and elastic properties as the particles being poured.

Most of our results are for packings that are spatially periodic in the horizontal plane, i.e. we ignored the effects of sidewalls. Because of this, the pressure in a packing does not saturate with depth. Therefore, to make a direct comparison with experiment, our definition of the average force $`f_{n,t}F_{n,t}/\overline{F}(z)_{n,t}`$, was normalized by $`\overline{F}(z)_{n,t}`$, the average contact force at a depth $`z`$ in the packing. The generation of these packings is fully discussed in Ref. . We also compared results for packings poured into a cylindrical container with ‘flat’ walls and the same properties as the particles. In this case, there is no need for depth-average normalisation, as the walls carry a significant fraction of the weight of the system . Results for the periodic packings with depth-average normalisation are consistent with the cylindrically confined packings and experiment. As a consequence, the depth-average normalisation proves to be the correct method for dealing with periodic packings.

## III Results

### A Force Distributions

Force distributions in all granular packings exhibit several general features. Measurements of the distribution of normal contact forces $`P(f_n)`$, for granular packings that are either free-standing under the influence of gravity  (as we simulate here), confined packings that have been loaded (as in experiments) , or axially compressed systems (as in previous simulation studies and experiment) , as well as a lattice model , all purport exponential tails in $`P(f_n)`$ at large forces (typically for $`f_n>1`$). Mueth et al.  used an empirical fit to their experimental data of the form,

$$P(f_n)=a(1be^{f_n^2})e^{\beta f_n},$$
(3)
and found $`a=3.0`$, $`b=0.75`$, and $`\beta =1.5\pm 0.1`$ for loaded glass spheres confined in a cylindrical container.

In Fig. 1 we show our computations of the force distributions for the normal contact force $`f_n`$ for different systems. In Fig. 1(a) we see that the form of $`P(f_n)`$ is the same for both Hookean or Hertzian contact force laws. Varying the system size has no effect (other than improving the statistics of the data). Similarly, in Fig. 1(b) computations of $`P(f_n)`$ in the bulk of a periodic or confined system, at the base of the periodic system, or at the sidewalls of the cylinder are indistinguishable. Recent 2D simulations have shown that $`P(f_n)`$ at the base can depend on the properties and geometry of the base . Computations of $`P(f_n)`$ for those particles in contact with the flat base and at the side walls also show the generic form seen in the other data although the statistics here are poor due to the number of contacts in the plane ($`10^4`$) compared with the number of particle-particle contacts in the bulk ($`10^510^6`$).

To compare with existing experimental data, we fit Eq. 3 to our data for the largest system. We show this comparison in Fig. 2. The $`P(f_n)`$ computed over all contact forces is denoted by the solid circles in Fig. 2 with best fit parameters $`a=2.55`$, $`b=0.65`$, and $`\beta =1.35`$, agreeing well with Eq. 3 up to $`f2`$, but falling off more quickly than Eq. 3 for large $`f`$. We find a better fit to Eq. 3 if we exclude the data for $`F_n<mg`$, i.e. ignoring all data in the limit $`f0`$, essentially mimicking the finite resolution in experiment. This alters the average value such that our original data set has now been ‘squeezed’ together. We denote this data as the partial set in Fig. 2. The fit to Eq. 3 with $`a=3.1`$, $`b=0.78`$, and $`\beta =1.55`$, is much better than when data for small forces is included. Our simulation data is in good qualitative and quantitative agreement with previous experimental results  and similar to Radjai et al. .

The empirical fit of Eq. 3 is poor for the total bulk $`P(f_n)`$ over a large range of the data and we only achieve agreement by renormalising our data, using the partial data set in Fig. 2. We also note that on closer inspection of existing simulation and experimental data, whether the tails of $`P(f_n)`$ are truly exponential or not is questionable and may be an indication of the averaging technique used in computational studies .

The distribution of tangential forces $`P(f_t)`$ is shown in Fig. 3. In comparison with the normal forces, $`P(f_t)`$ decays more slowly than $`P(f_n)`$. Fitting Eq. 3 to the bulk data for the largest system ($`N=128000`$), we find good agreement with $`a=2.5`$, $`b=0.7`$, and $`\beta =1.4`$.

While there is clearly some agreement on the behaviour of $`P(f_n)`$ for large $`f_n`$, the characteristic nature of the small force region of $`P(f_n)`$ remains in dispute. Experimental data show that $`P(f_n)`$ approaches a finite value as $`f_n0`$. However, some numerical works have suggested that $`P(f_n)0`$ for small $`f_n`$ . In Fig. 4 we show the small force region of $`P(f_n)`$ for packings identically prepared but with different particle friction coefficients. We do not show the full $`P(f_n)`$ curve as friction only weakly influences the behaviour of the large-$`f`$ region. However, our large system size data suggests that the exponential tail becomes slightly steeper with decreasing friction, i.e. $`\beta `$ increases as $`\mu `$ decreases. The defining feature of these packings is that for purely frictionless systems, $`P(f_n)`$ shows a well-defined peak at small forces, while for $`\mu >0`$, $`P(f_n)`$ has an upturn at very small forces. The amplitude of this upturn increases with increasing friction coefficient.

The Chicago group  studied different particle packings where $`\mu `$ varied by a factor of approximately three. Within the resolution of their experiment they did not find any systematic trend with friction. Because of the higher resolution in simulation, the following comments are relevant to such studies: the fraction of particle-particle contacts, or bonds, experiencing small forces increases with increasing $`\mu `$ even though the total number of contacts decreases with increasing $`\mu `$ . Further study shows that the fraction of particles that are close to the Coulomb yield criterion $`F_t\mu F_n`$, i.e. those particle pairs that are most likely to undergo local plastic rearrangement, increase as $`\mu 0`$. Indeed, we have previously reported  that frictionless packings are always isostatic, whereas frictional packings are hyperstatic and this may be related to the behaviour of $`P(f_n)`$ at small $`f_n`$.

For completeness we show the corresponding distributions $`P(f_t)`$, for the tangential forces in Fig. 5. In this case, we do not find any significant systematic trend with $`\mu `$. The role of $`\mu `$ in the determination of $`P(f_{n,t})`$, is subtle. In frictionless packings, $`P(f_n)`$ does not show an upturn at small $`f_n`$, therefore the generation of this upturn in frictional packings comes from the very presence of the frictional forces $`f_t`$, which influence the nature of particle contacts such that $`P(f_n)`$ itself observes an upturn at small forces.

### B Force Correlations

The spatial force-force correlation function $`(r)`$ measures spatial correlations between forces separated by a distance $`r`$. We use the same definition as in Refs. ,

$$(r)\frac{_i_{j>i}\delta (|𝐫_{ij}|r)f_if_j}{_i_{j>i}\delta (|𝐫_{ij}|r)},$$
(4)
where $`𝐫_{ij}`$ is the distance between particle contacts $`i`$ and $`j`$, and $`f_i`$ is the normalized contact force acting at contact $`i`$. In experiment, spatial force correlations can, at present, only be measured at container walls: the points of force contact coincide with particle contacts at the container surface lying in a 2D plane. The minimum separation between measurements in experiment is coincident with the particle size, $`r_{\text{m}in}d`$. In a 3D packing, contact forces transmitted at the points of particle-particle contacts are only restricted by excluded volume effects. For monodisperse spheres in 3D the minimum separation, $`r_{\text{m}in}\frac{d}{2}`$. A locally four-particle pyramid configuration would give this minimum separation. Mueth et al.  found no evidence for spatial correlations between the contact forces within the resolution of their measurements. L$`\varphi `$voll et al. , using a different measuring technique, resolved their force data showing weak force correlations at the base of their packings which extend out to approximately five particle diameters . This may only come about from the induced order of the sample at the container wall.

Because of the restricted geometry of experimental measurements, we found it instructive to compare our computations of $`(r)`$ for the normal contact forces within the bulk of amorphous packings, and see how these might depend on $`\mu `$. For comparison we also computed the correlation function between tangential contact forces for $`\mu =0.5`$. In Fig. 6 we present the spatial force correlation function for a frictionless packing ($`\mu =0`$) and a frictional packing ($`\mu =0.50`$). Within the bulk of the packing, forces are correlated, but only over short distances, extending to less than three particle diameters in the bulk, indicative of the diffuse nature of the force transmission network. However, the effect of friction on these correlations is very weak, with the frictional packing exhibiting only a very slight increase in local correlation.

Similar to Mueth et al. , in Fig. 7 we also show the radial distribution function $`g(r)`$, between contact points inside the bulk of a frictionless ($`\mu =0`$) and a frictional packing ($`\mu =0.50`$). Clearly, the frictionless packing has a higher first peak, representative of the higher coordination of the frictionless packing compared with the frictional one , and also local correlations between the positions of the contact points are stronger in the case of the zero friction packing indicating a more ordered distribution of contact points in the system.

### C Contact Geometry

We have so far shown that computations of $`P(f)`$ for various particle parameters yield essentially the same data, except for small $`f`$. It is ironic then, that although the generic features of $`P(f)`$ are a signature of the granularity of the system, it offers little distinctive information on the grain-level properties of the packing. Keeping in the spirit of particle pair information, in Fig. 8 we show the probability distributions for particle-particle contact angles defined in the local spherical coordinate system that bonds make with respect to the vertical (parallel to gravity direction). In Fig. 8 we compare packings with different $`\mu `$ (=0, 0.1, 0.5) and found that the distribution of contact angles has only a weak dependence on friction indicating that all the systems locally appear similar. In all cases, the majority of contact angles lie in the range $`45^{}<\theta <90^{}`$.

It is a simple exercise to further compute $`P(\theta )`$ only between particles that carry a large force, i.e. to identify or distinguish between “weak” and “strong” forces, as Radjai and co-workers have done for compressed systems . In Fig. 9, we compare $`P(\theta )`$ computed between all contacting particle pairs and $`P(\theta )`$ computed for the subset of particles in contact whose normal contact force $`f_n>f_{\text{c}ut}`$, where $`f_{\text{c}ut}`$ is some given threshold value. Here we set $`f_{cut}=2.0`$, i.e. all particles whose normal contact force is greater than the twice the average. Resolving the contact angle distribution according to a force cut-off as in Fig. 9 reveals that high force-bearing clusters are more directional and the anisotropy grows with increasing $`f_{cut}`$ (not shown here).

### D Contact Network

The existence of heterogeneous force networks is supported by experimental visualisation. Photoelastic particle packings  or piles subject to local force perturbations , demonstrate inhomogeneity in the magnitude of the forces propagating through granular assemblies - “force chains”. However, it is still not clear how relevant these structures are in determining the stress state of the system. The 2D simulations of Radjai et al.  suggested for compressed granular packings, a distinction can be made between the “strong” force network, those particles in contact that carry a force greater than the average normal contact force, and the “weak” force network, the network of particles that experience a force smaller than the average. In some theoretical approaches, the strong force network is assumed to support all the stress in the system, with the weak force network acting merely as a supporting framework to this which can essentially be neglected .

To investigate the relative importance of the force networks, we computed the normal force that sub-networks of force chains contribute to the bulk average contact normal force. In Fig. 10, we varied $`f_{cut}`$ and then computed the fraction of bonds remaining in the force network whose contact force was greater than $`f_{cut}`$ (‘strong’ force network), and computed the contribution that this network made to the average force. The computation of the relative force network contributions in Fig. 10 indicates only a weak distinction between the ‘strong’ force network for particle contacts with $`f_n2`$, and a weak force network with $`f_n2`$, say. Therefore it is questionable whether the so-called strong network actually does carry most, if not all, of the stress in the system. For example, by going from one curve to the other as indicated by the arrow in Fig. 10, we find $`50\%`$ of contacts contribute approximately $`80\%`$ to the average contact force. This is a small distinction, and not nearly an order of magnitude difference between the two networks that one might expect if the strong forces dominated the weak phase.

A related question is the stability of the relative networks. The concept of fragility  suggests in the limit one can clearly distinguish the strong network from the weak network, the strong phase should be a minimally coordinated, particle network. For a 3D frictional packing this suggests a coordination number $`z=4`$ . To calculate the network-averaged coordination number of a subset of particles, the contacting neighbours of the chosen network need to be included. In Fig. 11 we draw a schematic for determining the coordination number given a sub-network of particles (denoted by the grey particles), knowing the list of network neighbours (white particles).

Computation of the coordination number $`z`$ for packings with different $`\mu `$, over a range of cut-off values $`f_{cut}`$ is shown in Fig. 12. The network-averaged coordination number of particle clusters, based on the forces that they carry, decreases monotonically from the bulk averaged coordination $`(f_{cut}=0)`$ to approximately $`z=1`$. It appears that $`f_{cut}2`$, represents some limit in the system in the sense that for $`f_{cut}>2`$, the average size of particle clusters contributing are particle pairs, i.e. the largest cluster that propagates large forces is only of size two.

Additionally, we provide examples of force network realizations. In Fig. 13 we show two force network configurations of a slice approximately four particle diameters thick taken from the centre of the large, frictional, periodic system ($`N=128000`$ and $`\mu =0.5`$). We only show bonds whose force is greater than twice the average. Figure 13(a) is the force network for the absolute normal contact forces F without depth normalisation for $`F_{cut}=2\overline{F}`$. This corresponds to a section through the middle of a wide sandpile. Figure 13(a) can be compared to the 2D experimental realization in Refs. . If we show all bonds, the force network is dense with many weak forces. This may be an indication of the relative sensitivity of the experimental visualisation technique which clearly cannot resolve the smallest forces. The number of large forces increases with depth giving a clear indication of the propagation of weight down the pile.

Figure 13(b) is the force network for the depth-normalized normal contact forces $`f`$ with $`f_{cut}=2\overline{f}`$. We find similar configurations for the cylindrically confined packing. Therefore this is the equivalent force network for a confined, frictional, unloaded packing. Because of the weight of the particles have been renormalized out of the force (mimicking walls that support forces), forces of all magnitudes are seen throughout the pack. In both cases, we find that extended force-bearing structures exist over a range of length scales, but do not necessarily transmit the largest forces only.

## IV Conclusions

We have demonstrated that large scale simulations of granular packings offer insight into the effects of particle friction on measurements of the distribution of particle-particle and particle-wall contact forces $`P(f)`$. Our detailed comparison between simulation and empirical fits, obtained from experiments , showed moderate agreement. However, we were only able to fit our data over the full range in $`f`$ after renormalising our data by neglecting the smallest forces in the system (using the partial data set). We reason that this is an appropriate way to account for the limited resolution in experiment. We also reiterate the fact that many simulation and numerical studies of force distributions do not show a clear exponential tail either at large $`f`$ and we believe this may partly be due to the resolution of very small forces in such computer experiments that affect the total normalisation parameters.

We were able to discern the influence that friction plays on $`P(f)`$ in the small force region. The fraction of particle-particle contacts that experience very small forces increases with friction even though the total number of contacts decreases with increasing $`\mu `$ . Excluded volume effects rather than the functional form of the force law appear to dominate the bulk behaviour of the system for dense packings. Our studies of very large systems show that the tails of $`P(f)`$ become marginally steeper with decreasing friction, i.e. $`\beta `$ in Eq. 3 increases as $`\mu `$ decreases. Our ongoing work on simulating confined packings will investigate some of these issues further .

The force-force spatial distribution function and contact point radial distribution function indicate that spatial correlations between the contact forces and the positions of the contact extends out only to approximately three particle diameters. This shows that force correlations dissipate quickly in the bulk and that the force transmission network propagates locally but becomes diffuse rapidly. On introducing a force cut-off scheme to analyse force-bearing structures, we found no clear evidence for distinct ‘weak’ and ‘strong’ force phases.

In general this discussion was only concerned with unloaded or unperturbed amorphous granular packings. In this sense we have provided information on the “reference state” of a granular material from a microstructural point of view. This state is rather insensitive to friction and is primarily determined by construction history . Although we have not investigated the perturbed state or response function of these systems , it is likely that particle properties play a much more significant role in the response of a granular system than in the static state. Some theoretical treatments  on force chain analysis may benefit from the information of this unperturbed system when calculating the resulting response of such a system under some force perturbation. In fact, the contact angle distribution in Fig. 8 seems to suggest that the ‘splitting angle’ $`\theta _s=90^{}\theta `$, in the language of Ref. , does seem lie predominantly in the range $`0<\theta _s<30^{}`$. Comparing the experimental visualisation in Ref.  and the theoretical model in Ref. , the force chain analysis can be thought of as a superposition of force chains on top of the background force network.

This work was supported by the Division of Materials Science and Engineering, Office of Science, U.S. Department of Energy. This collaboration was performed under the auspices of the DOE Center of Excellence for the Synthesis and Processing of Advanced Materials. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DE-AC04-94AL85000.