A procedure to join the force and volume ensemble statistical descriptions of granular media

Granular media consist of a large number of discrete particles interacting mostly through contact forces that, being dissipative, jeopardizes a classical statistical equilibrium approach based on energy. Instead, two independent equilibrium statistical descriptions have been proposed: the Volume Ensemble and the Force Network Ensemble. Hereby, we propose a procedure to join them into a single description, using Discrete Element simulations of a granular medium of monodisperse spheres in the limit state of isotropic compression as testing ground. By classifying grains according to the number of faces of the Voronoï cells around them, our analysis establishes an empirical relationship between that number of faces and the number of contacts on the grain. In addition, a linear relationship between the number of faces of each Voronoï cell and the number of elementary cells proposed by T. Aste and T. Di Matteo in 2007 is found. From those two relations, an expression for the total entropy (volumes plus forces) is written in terms of the contact number, an entropy that, when maximized, gives an equation of state connecting angoricity (the temperature-like variable for the force network ensemble) and compactivity (the temperature-like variable for the volume ensemble). So, the procedure establishes a microscopic connection between geometry and mechanics and, constitutes a further step towards building a complete statistical theory for granular media in equilibrium.


Introduction
Granular media, like sand, coffee grains or mineral rocks, consist of a large number of discrete particles interacting mostly through contact forces [1]. Since their macroscopic properties arise from microscopic interactions, they should be perfect systems for applying statistical mechanics; however, the interactions are dissipative, which makes them hard to describe by using energy and temperature as main variables [2]. Instead, two main statistical equilibrium approaches have been proposed: the Edwards' Volume Ensemble [3,4], which accounts for all grains' configurations in mechanical equilibrium inside a given volume, and the Force Network Ensemble, developed by Snoeijer, Tighe and coworkers [5][6][7], which takes into account all sets of contact forces in mechanical equilibrium with a given external stress. Nevertheless, a relationship between the compactivity χ and the angoricity α −1 (i.e. the temperature-like variables for those ensembles) has not been established.
Hereby, we propose how those two ensembles could be combined for the limit state of isotropic compression of a set of monodisperse spheres. This is done by looking at the statistical distribution of volumes for Voronoï cells with the same number of faces s and the statistical distribution of pressures for grains with the same number * e-mail: jsreyl@unal.edu.co * * e-mail: jdmunozc@unal.edu.co * * * e-mail: woquendop@unal.edu.co of contacts z. Next, an empirical relationship connecting s and z is established that allows to find a connection between χ and α −1 by minimizing the total entropy against the coordination number z.

Aste and Di Mateo's volume ensemble
Edwards' Volume Ensemble takes the volume as the main thermodynamic variable, with an entropy accounting for all possible geometric configurations inside a total volume V T obeying volume exclusion and mechanical stability with external constraints, i.e. all possible jammed states [4,8]. In a illuminating proposal by Aste and Di Mateo [9] the volume V T is divided into C statistically independent elementary cells of volume v i > v min , where v min is a minimum volume due to steric exclusion. Because in isotropic compression the shape of such elementary cells should not be relevant, their volumes will be the degrees of freedom (DOF) for the system. Thus, we can define a phase space of C dimensions where each coordinate is the excess volume v i = v i − v min for each elementary cell. Then, it is possible to derive [9,10] that the probability for a singe elementary volume to take the specific value v is where v = V T /C is the average volume per cell. So, the entropy for a single cell can be calculated as ; and since the entropy is an extensive quantity, the entropy for all C independent cells is where Λ is a reference length 1 . The temperature-like variable 1 χ = ∂S ∂V T ⇒ χ = V T C − v min coincides with Edward's compactivity and indicates how much free volume is there per elementary cell. Although we don't know for sure what those elementary cells are, by performing a Voronoï tessellation and assuming that each tessellation cell contains approximately the same number of k v elementary cells, it follows that Voronoï cells volumes V voro should exhibit a k − gamma distribution [11], with V min = k v v min , a distribution that has been found in many experimental and computational experiments [12,13], (with k v ≈ 12 for monodisperse cells in the limit state of isotropic compression [10]). This way, the average V voro and the variance σ 2 V of that distribution can be used to calculate the compactivity of the system χ as well as the number of elementary cells per tessellation cell k v , Furthermore, Oquendo et. al [10] derived a state equation between the compactivity and the packing fraction φ for the system which is satisfied with A ≈ 0.04481d 3 for the limit state of isotropic compression and V min /v grain ≈ 1.3250 for monodisperse spherical packings [14,15].

Force network ensemble
A second major approach to the statistical mechanics granular media is the Force Network Ensemble (FNE) [5][6][7], which is based on the hyperstaticity of dense states, i.e., that the number of contacts per grain z is larger than the one z iso necessary to solve the equilibrium and external stress constraints uniquely and, therefore, there are multiple valid force networks f for the system. The number of degrees of freedom for the contact force variables is In the canonical ensemble under isotropic compression at pressure p (that is, when the stress tensor is σ α β = pδ α β ), each force network appears with a probability where α −1 is the angoricity [3] and Ω(p) is the number of force networks fulfilling those constraints. Since Ω(p) indicates the content of the solution space spanning k f dimensions and p can be considered as the typical diameter of one dimension [5], then Ω(p) ∝ p k f , and the following equation of state fulfills: The FNE has been widely studied by theoretical, experimental and computational works [16][17][18]. In a recent Monte Carlo study of the FNE on a single grain at constant angoricity, Cardenas et. al. [19] found that the pressure per grain, calculated as the sum of the normal forces per particle p i = z i j=1 | F n i j |, exhibits a k-gamma probability distribution, where k → k f for low angoricities (α −1 < 10 −2 ), satisfying the equation of state (7). This result suggests that the pressure can be considered as the sum of k f independent force variables f i with exponential distribution P( f i ) = αe −α f i . Thus, the entropy for one force variable will be S f 1 = 1 + ln α −1 p 0 , with p 0 a reference pressure 2 , and the total FNE entropy would be

Our work
To study a possible connection between the pressure per grain and the volume of the Voronoï cells in the limit state of isotropic compression we performed DEM simulations with the software LIGGGHTS [20][21][22] on a 3D cubic packing of 80000 monodisperse spheres. The contact forces are Hertzian plus a dissipative term as in Ref. [20], with time step ∆t = 10 −5 s. The limit configurations are obtained through four sequential procedures : A random insertion of particles inside a cubic box with initial length L i = 0.3m to create a dilute granular media with packing fraction φ ≈ 0.2, a particle grow for t grow = 0.5s until φ ≈ 0.3 in a process analogue to the Lubachevsky-Stillinger's algorithm [23], a relaxation for t relax = 0.2s until the particles take fixed positions and their kinetic energy is dissipated and an isotropic compression of the walls for t compress = 6s with inertial number I =˙ d ρ m P wall = 5 × 10 −4 (where˙ is the deformation ratio, d the particle diameter, ρ m the mass density of the system and P wall the external pressure applied on all walls). The simulation was run five times: one for each external pressure P wall = 13.9, 27.8, 55.5, 111.1 and 222.2Pa, and the wall velocity v wall was set to satisfy the low inertial number. 2 We choose p 0 = eα −1 iso from Eq. 10 with S f k f = 0, i.e. an isostatic state. and q = (with σ 1,2,3 are the eigenvalues of the stress tensor 3 ) to assure an anisotropy ratio q p < 10 −2 for the limit state. The analysis was then performed on all particles at least three diameters away from the walls, and excluding rattlers. All other material properties are listed in Table 1.
Let us take P wall = 55.5Pa as an example. The set of all Voronoï cells follows a k-gamma distribution with k v ≈ 14 and χ/d 3 = 0.0146, as in Ref. [10]. Nevertheless, not all Voronoï cells are the same. When divided into subsets by their number of faces, they still exhibit kgamma distributions, all with almost the same compactivity (χ/d 3 ∼ 0.0085(4) for subsets with more than 1000 cells, which is different from the global one), but with different ks (Fig. 1, left). Indeed, the results show a linear relation between the number of elementary cells per Voronoï cell k v and the number of faces s of such cell, k v = 2.9(1)s − 19(2) 4 ( Fig. 1, right). In a similar manner, the pressures per grain follows a k-gamma distribution with k = 1.78 and α −1 = 4.78 × 10 −3 (Fig. 2, left). When classified by the number of contacts z, they still show k-gammas with almost the same angoricity (α −1 ∼ 0.0036(2) for subsets with more than 1000 grains, which is different from the global one), but with different ks, which in the limit of low angoricity should obey Eq. 3 The granular stress tensor was computed from the microscopic interactions as σ α,β = 1 2V i j f i j,α r i j,β where r i j is the vector joining the centers of grains i and j and f i j is the interaction force between them. 4 Values between parentheses are 1 sigma error bars (7). We observe the linear relationship k f = 0.63(9)z − 0.5(3), independent of the pressure applied on the walls on the range 13.9Pa < P wall < 222Pa (Fig. 2, right). The connection between both ensembles is established by linking the number of faces s for the Voronoï cells and the number of contacts z for the particles. The joint distribution of s and z is almost the same for all P wall on the established range (Fig. 3,left). Although a particle with a given number of contacts can have a Voronoï cell with different number of faces (and vice-versa), for every z there is a preference for a given s. Indeed their averages follow linear relations (Fig. 3, right). We propose s = −0.30(1)z + 15.9 (9), where the negative sign implies that adding contacts reduces the number of faces. By assuming that pressure and volume are uncorrelated variables (a reasonable assumption for very rigid spheres), the total entropy S Total will be the sum of those for the volume and the force network ensembles (Eq. (2,10)), where C/N = k v is the number of elementary cells per Voronoï cell and k f /N = k f is the number of force degrees of freedom per particle. Using the linear expressions before, the parameters k v and k f can be written as function of the number of contacts: with M v = 0.87 (8), B = 27 (7), M f = 0.63 (9)) and D = 0.5 (3). Then, the total entropy can be written as a function of the coordination number z, The force entropy grows with the number of contacts, while the volume entropy decreases; therefore, in a statistical equilibrium state we could expect entropy to remain almost constant (i.e. maximized) with the number of contacts. Figure 4 (left) shows that this is indeed the case, with a total variation around 3% for the five values of P wall . Finally, maximizing the total entropy against the number of contacts z, gives an expression relating angoricity and compactivity, that is, an equation of state for a granular medium of rigid monodisperse spheres under low pressure isotropic compressions. Fig. 4 (right), shows that the logarithms of the angoricities for the subsets classified by z and the logarithms of the compactivities for the subsets classified by s fulfill a linear relation, as expected from Eq. (14). The positive slope of 2.94±1.21 is close enough to the expected value M v /M f = 1.38 to encourage us for developing more research in the future.

Conclusions
Summarizing, this work proposes a procedure to join the force network ensemble on a single grain and the volume ensemble of the Voronoï cell around it for a monodisperse spherical granular medium in the limit of isotropic compression, just by establishing an empirical relation between the number on contacts on the grain with the number of faces of the Voronoï cell. Using this relation we express the total entropy of the system as a function of the number of contacts, which is expected to be approximately constant in equilibrium. From there we propose a state equation between the average pressure per particle and the packing fraction, to be verified in future works. The procedure constitutes a further step to reach the holy grail of granular media, i.e. to develop a comprehensive statistical theory of granular media joining the volume ensemble and the force networks, clearing a way towards a theoretical description of dense granular media in static equilibrium.