Contact forces distribution for a granular material from a Monte Carlo study on a single grain
Manuel A. C\'ardenas-Barrantes, Jose Daniel Mu\~noz, William F., Oquendo

TL;DR
This study introduces a Monte Carlo method to analyze force distributions on a single grain in granular materials, revealing gamma-distributed pressures and an equipartition relation, advancing understanding of force networks at the grain level.
Contribution
Develops a Monte Carlo sampling technique for the force ensemble on a single grain, demonstrating gamma-distributed pressures and an equipartition theorem in granular systems.
Findings
Pressure per grain follows a gamma distribution.
Average pressure satisfies an equipartition relation.
Force variables are independent and exponentially distributed.
Abstract
The force network ensemble is one of the most promising statistical descriptions of granular media, with an entropy accounting for all force configurations at mechanical equilibrium consistent with some external stress. It is possible to define a temperature-like parameter, the angoricity {\alpha}^{-1}, which under isotropic compression is a scalar variable. This ensemble is frequently studied on whole packings of grains; however, previous works have shown that spatial correlations can be neglected in many cases, opening the door to studies on a single grain. Our work develops a Monte Carlo method to sample the force ensemble on a single grain at constant angoricity on two and three-dimensional mono-disperse granular systems, both with or without static friction. The results show that, despite the steric exclusions and the constrictions of Coulomb's limit and repulsive normal forces,…
| a) | |||||
| 5.0 | 10.0 | 20.0 | 50.0 | ||
| 3 | 1.04 [1.02] | 1.01 [1.00] | 0.98 [0.99] | 0.97 [0.98] | |
| 4 | 1.04 [1.04] | 1.04 [0.99] | 0.99 [0.98] | 0.96 [0.97] | |
| 5 | 1.02 [1.04] | 1.05 [1.00] | 1.00 [0.99] | 0.96 [0.96] | |
| 6 | 1.04 [1.05] | 1.05 [1.00] | 1.00 [0.99] | 0.95 [0.96] | |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
∎
11institutetext: Manuel A. Cárdenas Barrantes22institutetext: 22email: [email protected]
1 Simulation of Physical Systems Group, Department of Physics, Universidad Nacional de Colombia, Carrera 30 No. 45-03, Ed. 404, Of. 348, Bogota D.C., Colombia.
33institutetext: 2 Department of Mathematics, Physics, and Statistics, Faculty of Engineering, Universidad de la Sabana, Km 7 Autopista Norte de Bogota, Chía, Colombia.
Contact forces distribution for a granular material from a Monte Carlo study on a single grain
Manuel A. Cárdenas-Barrantes 1
Jose Daniel Muñoz 1
William F. Oquendo 2
(Received: date / Accepted: date)
Abstract
The force network ensemble is one of the most promising statistical descriptions of granular media, with an entropy accounting for all force configurations at mechanical equilibrium consistent with some external stress. It is possible to define a temperature-like parameter, the angoricity , which under isotropic compression is a scalar variable. This ensemble is frequently studied on whole packings of grains; however, previous works have shown that spatial correlations can be neglected in many cases, opening the door to studies on a single grain. Our work develops a Monte Carlo method to sample the force ensemble on a single grain at constant angoricity on two and three-dimensional mono-disperse granular systems, both with or without static friction. The results show that, despite the steric exclusions and the constrictions of Coulomb’s limit and repulsive normal forces, the pressure per grain always show a gamma distribution with scale parameter and shape parameter close to , the number of degrees of freedom in the system. Moreover, the average pressure per grain fulfills an equipartition theorem in all cases (in close parallelism with the one for an ideal gas). These results suggest the existence of independent random variables (i.e. elementary forces) with identical exponential distributions as the basic elements for describing the force network ensemble at low angoricities under isotropic compression, in analogy with the volume ensemble of granular materials.
Keywords:
Force network ensemble Contact force distribution Single grain Monte Carlo method Pressure Angoricity
1 Introduction
Granular media, like sand, rice, coffee grains, soils and powders, are relevant in everyday’s life and of paramount importance on industrial applications. However, there is still no general theoretical description for them. Most of our knowledge is ciphered in empirical constitutive equations relating the stresses, the strains and the energies inside the material Jacques2000 through many parameters whose microscopic origins are not always clear. But still there is a beacon of hope in statistical mechanics. Granular media are systems of many particles with a relevant interest in macroscopic quantities like volume fractions, strains and external stresses; so, they would be perfect candidates for a statistical mechanics analysis Henkes2009 ; cowin1979 . The task, nevertheless, is not easy, because granular media are dissipative systems. One major statistical mechanics approach is the force network ensemble Hecke2004 ; Vlugt2011 ; Ellenbroek2004 , which considers all possible configurations on a fixed contact network among grains that fulfill the requirements of mechanical equilibrium in agreement with some imposed external stress. Many configurations are possible, because the number of variables (contact forces) usually is greater than the constraints (equations of mechanical equilibrium per grain and values of external stress), in what is called hyperstaticity. A force entropy is defined by counting all possible force networks compatible with the external stress, and a temperature-like quantity, the angoricity is defined to relate that entropy with that external stress Edwards1998 . For the case of isotropic compression, this external stress is characterized by a single parameter, the pressure , and the angoricity is a scalar variable, .
There are many theoretical, computational and experimental works on the force network ensemble. For instance, F. Radjai and coworkers Radjai1990 found by contact dynamics simulations that forces below the mean distribute like a power law, and forces above, like an exponential decay. Also, Metzger and Donahue Metzger2005 suggested that forces follow a Bose-Einstein distribution. Later on, Thige, Snoejder, Vlugh and coworkers derive from dimensional considerations an equipartition relation for frictionless grains, which describes the expected value of the pressure as a function of the excess correlation number and the angoricity Vlugt2011 . By using a clever Monte Carlo algorithm (the wheel moves) to sample a 2D triangular array of monodisperse disks under isotropic compression Schaeffer2005 , they found that the correlation length in such a system is of the order of a single diameter Ellenbroek2004 , and that correlations among grains can be neglected in many cases Vlugt2011 . This exciting result suggests that the main behavior of the entire system could be understood from the behavior of the forces on a single grain. Indeed, the granocentric model proposed by Maxime Clusel and co-workers Clusel2009 ; Corwin2010 ; Puckett2011 manages to capture the essential properties of the dense granular material, such as the global density, from the statistics of the available space and the ratio of contacts to neighbors around a single grain, obtained from experimental measurements.
The present work develops a Monte Carlo procedure to sample the force network ensemble on a single grain, which is assumed to be part of a mono disperse granular system, both in 2D and 3D, with either frictionless or frictional interactions, and uses this method to investigate the ensemble. The method does not fix a pressure on the grain, but an angoricity (canonical ensemble), and only accepts mechanically stable configurations obeying steric exclusions and using non-cohesive normal forces and frictional forces (if any) below the Coulomb threshold. In all cases studied, the pressure follows a gamma distribution and the angoricity and pressure per grain fulfill an equipartition relation that, when extrapolated, reproduces the functional relation found by Thige and Vlugh for the whole system Vlugt2011 . This would validate the hypothesis that, under certain conditions, the statistics of forces on each grain can be considered as independent and shows that steric exclusions and force constraints do not have a measurable effect on these results.
2 The force network ensemble
Following Thige and Vlugh Vlugt2011 , let us consider a granular sample of grains in a two-dimensional box, in quasi-mechanical equilibrium with some external stress. The statements of mechanical equilibrium per grain and consistency with an external stress are given by
[TABLE]
where is the contact force acting on grain i by grain j, is the vector from the center of grain to the center of grain and is the volume (two-dimensional volume) of the packing. Eqs. (2) can be expressed in matrix form as
[TABLE]
Both arrays and are defined by the packing’s geometry. The vector gathers the independent components of the stress tensor, , and vector includes all components of the interacting forces among grains (that is, the number of variables). Besides, one says that a contact network is isostatic if there is only one force network satisfying Eq. (2), and hyper-static otherwise. Under isotropic pressure, , the density of states for a fixed value of pressure can be obtained as
[TABLE]
where the integral runs over all contact forces. The three factors , and assure for the mechanical constraints. The Heaviside step function assures that the contact forces are repulsive ones.
The hyperstaticity can be characterized by the mean excess coordination number, . For instance, let us consider frictionless systems in 2D. In this case, the total number of equations is . Thus, the system is hyper-static when , where ( for large systems) is the coordination number of the isostatic system (assuming no periodic boundary conditions). Since each contact has a single variable shared by two grains, the number of excess variables is just .
When , any force network satisfying Eq. (2) can be written as Schaeffer2005
[TABLE]
where is a particular solution to Eq. (2), are base vectors spanning the null space of solutions and are the components of in such a base. For 2D frictionless systems, for instance, the dimension of the space of all possible force networks is .
The maximum entropy principle can also be used to build up a canonical ensemble, with entropy
[TABLE]
where takes values of 1 or 0 if the force network fulfills or not the mechanical constraints. Under isotropic compression, the external stress is characterized by a single parameter, the external pressure . A force with magnitude appears in this canonical ensemble with probability , with the partition function and , a scalar Lagrange multiplier, whose inverse is called angoricity Hecke2004 . The average pressure can be computed as usual,
[TABLE]
Because the space of force networks is a convex polytope in dimensions, with the linear dimension of the polytope, one can assume Vlugt2011 ; therefore, by Eq. (6),
[TABLE]
in the thermodynamic limit.
From a microscopic point of view, one can define a local pressure on a single grain as the sum of all normal forces acting on it. The density of states for this pressure is
[TABLE]
where , and assure for mechanical equilibrium, external stress and normal repulsive forces, as before. For a single grain with frictionless contacts in dimensions, for instance, is a vector of components (that is, the number of variables is ) and the number of mechanical constraints is (one per dimension, to assure equilibrium, plus one to fix ). Thus, the configuration space at fixed pressure has dimensions Vlugt2011 , and one could assume
[TABLE]
which has been verified by molecular dynamics simulations Hecke2004 .Then, from the equipartition function for a single grain,
[TABLE]
one obtains Eq. (6),
[TABLE]
Similar results follow for frictional contacts or 3D ensembles, but with other values of . These are the relations we would like to verify by using Monte Carlo.
3 2D-Monte Carlo sampling on a single grain
3.1 2D Frictionless grains
Consider a single 2d-grain with z contacts, which is assumed to be part of a two-dimensional mono-disperse set of frictionless disks. The Monte Carlo sampling starts by setting the angular positions () for the z contacts. The main restriction here is steric exclusion, which also limits the coordination number z to a maximum of six. In addition, because we are just interested in hyperstatic systems, the minimum coordination number z is three, otherwise the system would be completely defined by the constraint of mechanical equilibrium. We consider two different possibilities for the contacts: either at fixed positions, equally spaced across the circle (called fixed) or randomly chosen from all possible configurations obeying steric exclusion (that is, with angles larger than between consecutive contacts) and using contact configurations per run. Next, we assign an initial set of non-cohesive forces for the contacts satisfying the constraints Eq. (2) of mechanical equilibrium. This is done by choosing contacts at random and choosing for them. The normal forces for the other two contacts are obtained from the mechanical equilibrium equations Eq. (2),
[TABLE]
The pressure on the grain is computed as .
Once the initial configuration is set, we start a Metropolis sampling scheme at constant angoricity, as follows: A new force configuration is generated by choosing contacts at random and changing the value of its normal force in a small amount , randomly generated on the interval ( for our simulations). The other two normal forces are modified to assure that mechanical equilibrium Eq. (3.1) is fulfilled. Then, the restrictions of non-cohesive forces are checked: If any new force is negative, the force configuration is rejected and a new one is randomly generated; otherwise, the new value of the pressure is computed. At this point, we introduce the Metropolis acceptance rate: If the pressure change , the move to the new configuration is accepted; otherwise, it is accepted only if a random number is less or equal than . This process is repeated as many times as necessary to reach equilibrium, i.e. when the pressure starts to fluctuate around a mean value. Typical equilibrium times are around time steps for , and lower for higher angoricities. Sampling starts after two equilibrium times. Once equilibrium is reached, the time correlation function for the pressure was computed and the correlation time, , measured. Typical correlation times were around for and smaller for higher angoricities. Samples are taken every .
Fig. 1 shows the pressure distribution for different values of angoricity and contact numbers (with both fixed and random contact positions). All curves are well fitted by Gamma distributions,
[TABLE]
where , known as the shape parameter, and , known as the scale parameter, are estimated from the average pressure and the variance in the histogram as
[TABLE]
By plotting against (Fig. 1, inset) we found that the parameter reaches a maximum for low and starts to decay. For low angoricities (large values) , which is the number of degrees of freedom in this one-grain system. From this result and Eq. (13), we conclude that the the pressure on a single frictionless grain distributes like
[TABLE]
This is the same relation Eq. (9) found by Tighe B P, Vlugt T. in Vlugt2011 for a bidimensional force packing.
Plotting against (Fig. 2) shows a linear relationship between these two quantities for all values of angoricity,
[TABLE]
and the relation is even better for low angoricities. Equation (16) is a clear representation of an equipartition theorem, in perfect parallel with the equipartition relation for an ideal gas, where the average pressure plays the role of the energy, and the angoricity, the one of temperature. This result can be extrapolated to the whole packing if one takes the number of degrees of freedom as the dimensionality of the space of hyperstatic solutions, , recovering the equipartition relation Eq. (11) proposed by Thighe and Vlugh.
3.2 2D Grains with static friction
To extend the previous analysis to a more realistic situation, grains with static friction, only few changes must be introduced. First, there are now two variables per contact: the normal and tangential components of the contact force. Second, the conditions for mechanical equilibrium are now three: two for the forces plus one for the torques,
[TABLE]
[TABLE]
[TABLE]
Third, there is an extra restriction: The absolute value of the tangential force at each contact cannot surpass the Coulomb limit, . The torque condition in the initial configuration is assured by starting with all and solving for two normal random forces, as before. The equilibrium conditions Eq. (17) are assured at each Monte Carlo step by solving for three random variables, instead of two. The small variations in the tangential forces are randomly generated from the interval to avoid an extreme number of rejections, that is ten times smaller than the one for the normal forces. The Coulomb restriction is included when checking for non-cohesive normal forces. Anything else is unchanged.
We see again a gamma distribution for the probability in all cases (Fig. 3). The parameter k shows the same behavior as before, but reaching a different value for low angoricities, that is the new number of degrees of freedom in the system. Thus, the probability function for a grain is now proportional to
[TABLE]
which extends the relationship found by Tighe and Vlugt Vlugt2011 Eq. (9) to the frictional case.
Plotting the average pressure times the inverse angoricity against the coordination number (Fig. 4) evidences again an equipartition equation relating these two quantities,
[TABLE]
where just the number of degrees of freedom has changed. Furthermore, Eq. (19) shows to be the same for all simulated friction coefficients.
4 3D-Monte Carlo sampling on a single grain
The Monte Carlo procedure on a single three-dimensional grain runs very similar to the two-dimensional case. The first difference is that the contacts are chosen on the sphere, either in regular configurations or at random positions. In this case, once the first contact is set, the second one is accepted only if steric exclusions are respected, and so on. The maximal coordination number imposed by steric exclusions is ; nevertheless, obtaining random configurations for is very unlikely; thus larger coordination numbers were studied on regular configurations only.
4.1 3D Frictionless grains
In frictionless systems there is one variable per contact, and the constraints for mechanical equilibrium are now three,
[TABLE]
where is the -th component of the normal force . Thus, the isostatic limit would be . Nevertheless, such a configuration is only possible for a single unstable configuration: a planar symmetrical array of three spheres. Stable configurations are reached only for coordination numbers . The Monte Carlo procedure starts by fixing contact positions, either regular or at random. For the initial configuration, normal forces are chosen randomly and set to , whereas the other three are computed from Eq. (20). If all computed forces are repulsive, the configuration is accepted, the new pressure is computed and the Metropolis step is completed as before.
The pressure shows again gamma distributions (Fig. 5). At low angoricities, the parameter reaches a new value, , that is the number of degrees of freedom in the system. Thus Eq. (13), the probability function for a grain is now proportional to
[TABLE]
and, therefore Eq. (6), the system fulfils an equipartition relationship
[TABLE]
This is exactly what is observed in our numerical simulations (Fig. 6).
4.2 3D Grains with static friction
Introducing friction requires few changes. There are now three variables per contact : the magnitude of the normal force and the two components of the tangential force . The conditions for mechanical equilibrium are now six: three for the forces plus three for the torques:
[TABLE]
[TABLE]
where indexes the three components of each vector and goes from the centre of the grain to the -th contact point. The initial condition is chosen the same as in the 3D frictionless case (i.e. initial tangential contact forces are set to zero). Also, new configurations are only accepted if both all normal forces are repulsive ones and the absolute value of the tangential force at each contact does not surpass the Coulomb limit, . Except by these small changes, each Monte Carlo step is completed as before.
The pressure histogram shows again a gamma distribution (Fig. 7). At low angoricities (high values of ) the parameter goes to , i.e. the new number of degrees of freedom in the system. Thus, the probability function for a grain is proportional to
[TABLE]
and the system fulfills an equipartition relationship
[TABLE]
almost the same for different friction coefficients (Fig. 8).
5 Discussion and Conclusions
We developed a Monte Carlo method to sample the ensemble of force configurations on a single grain, both with or without friction. The grain is assumed to be part of a monodisperse granular medium, either in two or three dimensions, under isotropic compression. The contacts on the grains are chosen either at random or in a regular configuration, always respecting steric exclusions. The set of mechanical stable force configurations is sampled by a Metropolis Monte Carlo algorithm at constant angoricity. New configurations are accepted only if all normal forces are repulsive and all tangent forces are below the Coulomb’s static frictional limit.
Our results show, first, that the pressure on the single grain (that is, the sum of all normal forces) follows a gamma distribution in all cases. This coincides with the proposal of T. Aste an T. Di Matteo Aste_Emer2008 for the distribution of Voronoï cells in a granular media in the frame of Edward’s statistical mechanics of volumes Aste_Emer2008 , which has also been confirm for grains under isotropic compression by Oquendo and co-workers Oquendo2016 . Actually, the procedure we used to find and are the same employed there. The insets in Fig. 1, 3, 5 and 7 show that, in all cases, the scale parameter in the gamma distribution is very similar to the number of degrees of freedom for the forces on the grain and, that in the limit for low angoricities (high values of ). This result could be observed as a natural consequence of the sampling Metropolis Monte Carlo procedure, which moves variables and accepts or rejects new configurations to fulfill a detailed-balance among individual configurations with exponential probability distributions. Nevertheless, the point here is that this result is not altered neither by steric exclusions nor by the restrictions on the forces (non-cohesive forces and Coulomb’s static limit for the tangential forces). Moreover, the gamma distribution we postulate here for the force network ensemble is compatible with all previous descriptions Vlugt2011 ; Tighe2010 . Second, we also found in all cases (Fig. 2, 4, 6 and 8) that the average pressure fulfills an equipartition-like relation . These two results are valid for all systems studied, either with or without friction and both in two or three dimensions. These findings confirm for a single grain the proposal by Tighe and Vlugt Vlugt2011 for the whole system, postulating an equipartition-like relation and that , with the number of degrees of freedom in the system.
Actually, this two results combine to give interesting consequences. First, because in a gamma distribution the shape and scale parameters are related by Eq. (14), it implies that scale parameter and angoricity are equal, . Table 1 shows that this relation is fulfilled in all studied cases, both in 2D and 3D. This result suggest a novel method to measure the angoricity in monodisperse granular media under isotropic compression, just by looking at the scale parameter of the distribution of pressures per grain in the ensemble Eq. (14), resembling the propose by T. Aste and T. Di Matteo for computing the compressibility of Edward’s volume ensemble for granular media from the distribution of volumes for Voronoï or Delaunay cells Aste_Emer2008 mentioned before.
It is well known DevroyeNon1986 that the sum of independent random variables with the same exponential distribution
[TABLE]
shows itself a Gamma distribution Eq. (13). That suggest that the pressure on a single grain could be considered as the sum of independent random variables, one for each degree of freedom, with distribution Eq. (26). Let us consider the 2D frictionless case, where all contact forces are normal. With contacts we would expect independent random variables with the same exponential distribution. If two contacts and are fixed at random, the variables
[TABLE]
are the degrees of freedom of the system. We found that, when positive, those variables follow exactly the exponential distribution Eq. (26) (Fig. 9). In contrast, the contact forces theirselves show, for low angoricities, a Gamma distribution with shape parameter (Fig. 9, inset), as the pressure per grain does (Fig. 1) and, Gamma distributions grow as power laws and decay as exponentials, in agreement with previous results Radjai1990 ; Metzger2005 . So, the variables would act as fundamental forces, similar to the fundamental volumes proposed by Aste and Di Matteo for the volume ensemble in granular media Aste_Emer2008 . If it were the general case, the force network ensemble in the limit of low angoricities would be considered as a set of independent random variables with distribution Eq. (26) and, all results by Tighe and Vlugt Vlugt2011 mentioned before could be derived. Moreover, by assuming that the magnitude of such elementary forces would be the only relevant variables in the system (as a first approximation for a monodisperse granular system under isotropic compression), the entropy for the force network ensemble would be
[TABLE]
with a reference pressure. If this were the general case would be an interesting topic for future research.
This work shows that the probability distribution and the equipartition relation proposed by Thige and Voigt for the pressure per grain Vlugt2011 can be obtained by extending the same assumptions to the force network ensemble on a single grain, and gives a further ground for that assumptions. Furthermore, our results suggest that, at low angoricities, not only grains can be considered independent form each other but, even more, that more elementary independent variables could exist, all with the same basic exponential distribution. These results constitute a further step in the understanding of the force network ensemble for granular media.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1(1) Aste, T., Di Matteo, T.: Emergence of gamma distributions in granular materials and packing models. Phys. Rev. E 77 (2008). DOI 10.1103/Phys Rev E.77.021309
- 2(2) Clusel, M., Corwin, E., Siemens, A., Brujic, J.: A granocentric model for random packing of jammed emulsions. Nature 460 , 611 (2009)
- 3(3) Corwin, E., Clusel, M., Siemens, A., Brujic, J.: Model for random packing of polydisperse frictionless spheres. Soft Matter 6 , 2949 (2010)
- 4(4) Cowin, S., Satake, M.: Continuum mechanical and statistical approaches in the mechanics of granular materials. J. Rheol. 23 (23), 243 (1979). DOI http://dx.doi.org/10.1122/1.549526
- 5(5) Devroye, L.: Non-Uniform Random Variable Generation, 1 edn. Springer-Verlag (1986)
- 6(6) Duran, J.: An Introduction to the Physics of Granular Materials, 2 edn. Springer (2000)
- 7(7) Edwards, S., Grinev, D.: Statistical mechanics of vibration-induced compaction of powders. Ph. Rev. E 58 (4), 4758 (1998). DOI http://dx.doi.org/10.1103/Phys Rev E.77.021309
- 8(8) Ellenbroek, G., Hecke, M., Snoeijer, H., Vlugt, H., Leeuwen, J.: Ensemble theory for force networks in hyperstatic granular matter. Ph. Rev. E 70 (6), 061,306 (2004). DOI http://dx.doi.org/10.1103/Phys Rev Lett.92.054302
