Eur. Phys. J. E (2014) 37: 81 DOI 10.1140/epje/i2014-14081-6
THE EUROPEAN PHYSICAL JOURNAL E
Regular Article
Ensemble theory for slightly deformable granular matter Ignacio G. Tejadaa Universit´e Paris-Est, Laboratoire Navier (UMR 8205), CNRS, ENPC, IFSTTAR, France Received 19 May 2014 and Received in final form 14 August 2014 c EDP Sciences / Societ` Published online: 25 September 2014 – a Italiana di Fisica / Springer-Verlag 2014 Abstract. Given a granular system of slightly deformable particles, it is possible to obtain different static and jammed packings subjected to the same macroscopic constraints. These microstates can be compared in a mathematical space defined by the components of the force-moment tensor (i.e. the product of the equivalent stress by the volume of the Voronoi cell). In order to explain the statistical distributions observed there, an athermal ensemble theory can be used. This work proposes a formalism (based on developments of the original theory of Edwards and collaborators) that considers both the internal and the external constraints of the problem. The former give the density of states of the points of this space, and the latter give their statistical weight. The internal constraints are those caused by the intrinsic features of the system (e.g. size distribution, friction, cohesion). They, together with the force-balance condition, determine which the possible local states of equilibrium of a particle are. Under the principle of equal a priori probabilities, and when no other constraints are imposed, it can be assumed that particles are equally likely to be found in any one of these local states of equilibrium. Then a flat sampling over all these local states turns into a non-uniform distribution in the force-moment space that can be represented with density of states functions. Although these functions can be measured, some of their features are explored in this paper. The external constraints are those macroscopic quantities that define the ensemble and are fixed by the protocol. The force-moment, the volume, the elastic potential energy and the stress are some examples of quantities that can be expressed as functions of the force-moment. The associated ensembles are included in the formalism presented here.
1 Introduction Granular materials are large configurations of discrete macroscopic particles of size larger than one micron [1] that display unusual behavior in comparison with that of solids, liquids or gases [1–3]. Since they are present in many different aspects of our daily lives, the establishment of a theoretical framework for their description has become an important issue. There are two important aspects that contribute to their unique properties: ordinary temperature plays no role, and interactions between grains are dissipative because of static friction and the inelasticity of collisions. One of the consequences of this is the possibility of obtaining different static packings of the same granular system under the same macroscopic constraints. If we focus on these packings, we can find different and inhomogeneous distributions of contacts, interparticle forces, voids, etc. [4–8]. Moreover, the mechanical behaviors of the packings are determined by their geometries, which are in turn adapted to the applied loads [9]. a
e-mail: [email protected]
Since granular systems consist of a lot of similar particles in which the interaction laws are rather known, a statistical mechanics approach might be feasible to describe them [10,11]. Edwards and Oakeshott were the first to do this [10]. Assuming that granular systems have entropy, they claimed that the volume plays the role of energy (V-ensemble [10,12–15]). There are other formalisms based on the energy of the whole system, such as [16–18]. It has been also suggested to consider the stress together with the volume (V–F-ensemble or full canonical ensemble) [19–21]. If the constraint of the volume is removed from the last, then the F-ensemble arises [22,23]. The force network ensemble theory [24–27] focuses on the statistical distribution of interparticle forces with a given contacts network. The aim of this paper is to propose a formalism for a general athermal ensemble theory that is valid to compare static packings of slightly deformable granular matter. Under some assumptions, a mathematical space is used to compare microstates and to obtain which the most probable distribution of some quantities is when some external constraints are imposed. Some of the features of this
Page 2 of 8
space are analyzed. Within this framework, some of the ensembles reported above, as well as the stress ensemble, can be analyzed.
2 Basic assumptions and definitions 2.1 Granular system, packings and protocol A granular system is a set of N discrete macroscopic particles that is determined by their features. For instance, for the systems analyzed here, these features are: the size and shape, the roughness and the mechanical properties of the material (e.g. density, constitutive relationships, strength criteria, etc.). For simplicity, without loss of generality, the formalism included in this paper is adapted for twodimensional systems, but the theoretical scheme is valid for three-dimensional systems too. The motion of each particle under the action of external and interaction forces is governed by classical mechanics. Only mechanical pairwise interactions are considered here. These include dissipative or viscous terms to represent inelastic collisions. Particles are deformable when they change their shape with boundary and body forces. Only small and elastic deformations are expected here. The contact theory (see [28]) can be used to find the interaction forces. In the case of two-dimensional elastic particles, k and l, the proposed relationship between the force (per unit of length) and the associated overlapping is linear (Hookean-like force): Fkl = −Kn δ kl xkl /xkl , with xkl = xl −xk and δ kl = D −xkl (for xkl ≤ D). The normal viscous force is proportional to γn (m/2)δ˙ kl , being γn the viscoelastic damping constant. Each static arrangement of a given granular system is referred to as packing. Static means that particles are force-balanced and at rest, so that state corresponds to a local energy minimum in the rugged potential energy landscape produced by the pairwise interactions. The packings may be in turn defined by their contacts and forces networks. Several forces networks may be associated with the same contacts network when the boundary forces are not fixed or when the contacts network is hyperstatic, which means that the number of contact points is higher than the strictly condition of static equilibrium (isotaticity). A protocol is a realization of different static packings of a system in a systematic manner (e.g. pouring, shaking, compressing, etc.). Although the obtained packings are static, it is actually a dynamic process in which kinetic energy is injected and dissipated. If microscopic degrees of freedom (e.g. phonons, electrons, etc.) were included in the model, then the energy dissipation would actually be a redistribution of grain-level energy to these microscopic degrees of freedom. This approach is thought for jammed states, although a precise definition of jamming is not attempted now. These states can be seen as those configurations in which each particle is in contact with its nearest neighbors in such a way that mechanical stability of a specific type is conferred to the packing [6]. Considering the jamming
Eur. Phys. J. E (2014) 37: 81
phase space [29], jammed packings of a granular system are close to the origin of the temperature-volume-stress space. However the limits of validity of this approach are not established yet. For practical purposes, all the packings analyzed here are highly hyperstatic. For each particle k in a packing Π, the corresponding region consisting of all points closer to its center than to the center of any other is itsso-called Voronoi cell. Its associated volume is V k , and k V k = V . In static packings, the equivalent stress field of each particle (the stress field averaged over its Voronoi cell) can be obtained from the mechanical forces that keep it k = (1/V k ) l Fikl rjkl where rikl is the in equilibrium: σij spatial position of the contact point between particles k k k and l (see [30]). The product Σij = σij · V k = l Fikl rjkl is referred to as force-moment or extensive pressure. The total force-moment and the average k stress of the whole system are given by Σij = k Σij and σij = Σij /V , respectively. Instead of the Cauchy components of the force-moment tensor Σij , it is interesting to use other variables that are functions of its invariants: I1 = Σxx + Σyy and 2 . The mean force-moment component I2 = Σxx Σyy − Σxy is defined as P = I1 /2 and expresses the level of the field. I 2 1 − I2 and the The deviatoric component is Q = 2 orientation Ω =
1 2
2Σ
xy arctan Σyy −Σ . P is the semi-sum of xx
principal force-moment values, i.e., the eigenvalues, and Q is equal to their semi-difference. The third component Ω ∈ [0, π) gives precisely the orientation of the field with respect to the external framework of reference, so it is not invariant. The relative deviatoric force-moment is a dimensionless quantity defined as Q = Q/P that gives an idea of how isotropic Q → 0.0 or anisotropic Q → 1.0 the field is. 2.2 Configurations On the level of the particles, many possible states of equilibrium exist (fig. 1). These states can be classified according to the number of forces, their relative angles, their value and direction, and the orientation of the set with respect to an external framework of reference. The term configuration is used here for a specific number of forces with a specific layout. Note that the orientation has been excluded in this definition (e.g. the states depicted in fig. 1C belong to the same configuration). It is worth noting that the same configuration α can often represent different states of equilibrium with different sets forces (fig. 2). Thus, it may represent different states of force-moment too. In the end, the possible configurations of a system are determined by the condition of force-balance l Fkl = 0 and by the internal constraints, i.e. those which are exclusively dependent on the features of the particles. If the system does not change, neither do its internal constraints nor do its possible configurations. Therefore the possible configurations of a system do not depend on the external
Eur. Phys. J. E (2014) 37: 81
Page 3 of 8
Fig. 2. Two possible local states of static equilibrium for a particle k that belong to the same configuration but have different sets of forces and represent different states of force-moment.
Fig. 3. An example of the same configuration with two different Voronoi cells. The second is only possible if the system is polydisperse. Fig. 1. Some examples of possible local states of static equilibrium for a particle k. A: With a different number of forces. B: With the same number of forces but a different layout. C: With the same layout but a different orientation.
constraints imposed to the packings that are being sampled. Some examples of internal constraints are: – In case of monodisperse systems (particles equally sized) the relative angle between two forces Δθkl = |θkl+1 − θkl | cannot be smaller than π/3. – In case of frictionless particles, all the forces must be central Fkl × xkl = 0. – In case of cohesionless materials all forces must be compressive. Fkl · xkl < 0.0. In configurations with a low number of forces (and especially when the system is polydisperse), it is possible to find local packings with the same layout but different Voronoi cell and thus different volumes (see fig. 3). However this effect is disregarded here. For simplicity, and assuming that packings are tight, a unique value of the volume Vα is supposed to be associated with each configuration1 . Nevertheless, in case of slightly deformable particles, this volume decreases as the forces are increased. 1 Otherwise some density of states functions should be introduced to consider all possible values between Vα,min and Vα,max .
It depends basically on stiffness (included in the contact law), and it is reasonable to assume that it only does with P (the variable that sets the force-moment level). Therefore Vα = Vα(P ) .
3 Elements of ensemble theory A general athermal ensemble theory is proposed for compressed, static and jammed packings of a given granular system of N slightly deformable particles that satisfy some macroscopic constraints. If N is an extremely large number, then it is possible to perform an analysis in the so-called thermodynamic limit. The specification of the number of particles and the actual values of the macroscopic quantities of a packing define its macrostate. Some macroscopic quantities, symbolized in general as M , are: the volume of the packing, V , the macroscopic stress, σij , the elastic potential energy that the packing stores, Ee , or the force-moment Σij = V σij . At grain level, however, a large number of possibilities still exist because at that level there will be in general a large number of different ways, W, in which the macrostate {M } can be realized (note that microscopic internal degrees of freedom are not included). Each of these different ways, i.e. each static packing (contacts and forces networks), specifies a microstate of the system. An ensemble E(N,M ) is precisely defined as a large number of copies
Page 4 of 8
of a system, each of which represents a possible microstate that the real system might be in or a large number of different packings of the same system satisfying the same macroscopic constraints. The exploration of microstates is not produced as time goes by (as it happens for thermal systems), because the energy scale kT is insignificant for granular media (at grain level) and because of the inelasticity of collisions and the non-conservative friction forces. Therefore a protocol is needed to explore the microstates. Identical systems are assumed to give the same average response to the same driving protocol (e.g. [31–34]). The postulate of equal a priori probabilities states that for a given macrostate of the system, and when there are no other constraints, after an appropriate protocol the system is equally likely to be in any one of these microstates. Appropriate protocol means precisely that both the forces and the contacts networks are constantly being changed and the system explores all its possible packings. It is assumed that microscopic internal degrees of freedom do not add new constraints (i.e. neither exclude microstates nor alter their statistical weights). In such case, the actual number of microstates W is just a function of the macroscopic variables: W = W(N,M ) . However, instead of trying to figure out what these whole packings are, we can focus on the local state of the particles. Then the following assumption is made: when all the possible packings of an ensemble, Π ∈ E(N,M ) , are considered, each particle k can be found in any equilibrium state that corresponds to a possible configuration α of the system, and there is not any a priori reason to favor any of these states more than any other. Furthermore, if the contact network were not changed during the driving process, it would be possible to establish some relationships between the values of the variables of the particles. Then the number of degrees of freedom of the packings would be equal to the number of free variables (those that can be modified to satisfy the macroscopic constraints). However, this is not the case. Here, all the possible packings of a system are being considered and, in consequence, there is no reason to establish any fixed relation between, for example, positions, neighbors or values of forces. In other words, for each contact network there are relationships that result on a known number of degrees of freedom (which depends on the hyperstaticity of the network and the variability of boundary forces), but when all of them are considered, particles are supposed to behave independently. There are no fixed relationships during the exploration of microstates. 3.1 Phase space The phase space is a mathematical space in which all the possible states of the system can be represented. Therefore the microstates can be compared there. Here the proposed coordinates are the configuration αk and the forcemoment state (P k , Qk , Ω k ) of each particle k. In this space some macroscopic constraints may be expressed as a function of these coordinates: M = M(αk ;P k ,Qk ,Ω k ) .
Eur. Phys. J. E (2014) 37: 81
Thus, this space is discrete in configurations and continuous in components of the force-moment tensor. The integration limits of the force-moment coordinates are those which guarantee that the stress field is compressive in any orientation. It can be proved that it happens when k k k k Σ k or, equivalently, Σxx , Σyy ∈ [0, ∞) and |Σxy | ≤ Σxx yy when P k ∈ [0, ∞), Qk ∈ [0, 1] and Ω k ∈ [0, π). In this space all the microstates of the system can be represented. For instance, it allows to locate particles that are locally packed at the same volume but are arranged according to different geometrical patterns, or particles that, although their arrangement follows the same geometrical pattern, have different values of the forces (and, consequently, of the force-moment). However, with this choice of coordinates, some information has been lost (e.g. the same configuration and force-moment state can be sometimes represented by different sets of values of the forces). It means that a flat sampling over force-moment coordinates is not taking into account all the possible states of equilibrium. This can be solved by using density of states functions gα(P,Q ,Ω) that are actually the result of a flat sampling over all possible local states of equilibrium. All the degrees of freedom are being considered within these functions. In particular, in two-dimensional problems the use of three components for the force-moment (P k , Qk , Ω k ) and one for the type of configuration αk does not absolutely mean that only 4N degrees of freedom are being considered. Although the establishment ab initio of gα(P,Q ,Ω) functions is complicated, this is not the objective, because they can be measured. Moreover, the definition of configurations {α1 , α2 , . . . , αr } proposed here is arbitrary. For example, all the configurations having the same number of forces could be gathered in an upper class (A1 = {α1 , . . . , αm }, A2 = {αm+1 , . . . , αn }, etc.) that does not distinguish between different layouts. Then there would be another classification {A1 , A2 , . . . , AZ } with different functions gA(P,Q ,Ω) . Even configurations may not be distinguished at all, but just using a global density of states function over the phase space g(P,Q ,Ω) . Using more configurations means that there are more functions to measure. However it allows to get more information about the most probable distributions (e.g. expected ratios between particles arranged as honeycombs with respect to those arranged with 5 forces). Although these functions still remain unknown, some symmetries can be found. For instance, consider the case of a monodisperse, frictionless and cohesionless system, and a particle k that is in a specific state of equilibrium under a set of forces {Fkl } and has an associated Voronoi cell of volume V k . It represent a force-moment state (P, Q , Ω). It can be easily proved that such state of equilibrium could also represent a force-moment state (P, Q , Ω + Θ), for any arbitrary Θ, just by rotating its layout an angle Θ. As configurations are defined independently of their spatial orientation, the density of states functions do not depend on the orientation of the field, so gα = gα(Ω) .
Eur. Phys. J. E (2014) 37: 81
On the other hand, if all the forces acting on the particle k are scaled with any arbitrary factor x, the particle remains in equilibrium (because the force-balance and the internal constraints are not affected by this scale transformation). If the variation of V k is disregarded, then the force-moment state changes from (P, Q , Ω) to (xP, Q , Ω). Therefore the possible number of states at level xP is at least the same as at level P . Indeed it usually grows with P , although the rate is uncertain. While a flat sampling over the force-moment variables Σxx , Σxy , Σyy leads to a growth that is proportional to P 2 (since dΣxx dΣxy dΣyy = 4P 2 Q dP dQ dΩ), a flat sampling over the values of the forces of each configuration does not. To illustrate this, consider a honeycomb configuration as the one shown in fig. 2. The equilibrium is produced by 6 forces {F k1 , . . . , F k6 } which follow a continuous uniform distribution U(0,Fmax ) . The minimum value is 0 because the system is cohesionless, while a maximum value is necessary because of the contact law (although usually F kl Fmax ). If the force-balance condition is imposed, then there are two additional constraints that make it possible to select randomly just four forces, say F k1 , F k2 , F k3 , F k4 ∈ U(0,Fmax ) , and compute k6 the other two, F k5 = F k1 + F k2 − F k4 and F = k3 k4 k1 k F + F − F . On the other hand, P = (R /2) l F kl or P = (Rk /2)(F k1 +2F k2 +2F k3 +F k4 ). Therefore, when a flat sampling over the space of forces is done, the statistical distribution of P for the honeycomb configuration results from the addition of independent random variables uniformly distributed. When F kl Fmax , as usual, i.e. for lower values of P , this distribution is proportional to P 3 . By contrast, in the case of a square lattice the number of states grows linearly with P (at low levels). And the same could be done with all possible configurations to compare them and to normalize the density of states functions. For those with five forces, it grows with P 2 , while the number of states is independent of P for configurations with three and two (opposite) forces. However, as the intention is to measure the functions, the important conclusion to extract is that the number of states grows monotonically with P and at rates that are different to P 2 . Furthermore, the distribution over Q values coming from a flat sampling over the space of forces is much more complex, but also feasible. This is, in the end, what makes a configuration more or less probable than others (when they are normalized by considering all possible configurations) and here is where the largest differences appear. For instance, in two-dimensional systems, the honeycomb represents by far more states than any other configuration. Under the approach presented here, this is the reason why these systems tend to crystallize. Considering all this, it seems convenient to express gα(Σij ) = hα(P ) · gα(Q ) . Under this scheme, gα(Q ) takes into account the internal constraints of equilibrium of each configuration while hα(P ) is the correction needed to represent properly the growth of the number of states with P .
Page 5 of 8
3.2 Canonical-like ensembles and equilibrium When a protocol is applied to a granular system, particles change their coordinates in the force-moment space with the cycles during the driving process, but it is not established how it happens (as, for example, Hamiltonian mechanics does in the phase space for classical ideal gases). Nevertheless the macroscopic constraints of the ensemble at the end of the driving cycles are fixed by the protocol. Then, the (thermodynamical) equilibrium is defined by reference to the macrostate (Boltzmann’s approach), which is uniquely given by a set of values of the macrovariables. The expected distribution in the equilibrium is obtained by maximizing Boltzmann’s entropy S ∝ ln W(N,M ) , i.e., by analyzing which the most probable distribution of an ensemble subjected to the same constraints is. Therefore, the problem is not mechanical but statistical, because the role of the dynamics during driving cycles is trivial for this definition. Moreover, since the exploration of mechanical stable states is not produced as time goes by (as thermal systems do), the ergodic hypothesis is not based in time averages but in ensemble averages (independent realizations of same process). Unfortunately, there is nothing like Liouville’s theorem to know how the evolution of the statistical distribution is when the packings are realized. This theory just gives the most probable distribution, but it does not necessarily remain constant with driving. It is worth mentioning that if the microscopic degrees of freedom were included, the total energy of all would be conserved and redistributed and that force equilibrium would happen when the total entropy (considering all the degrees of freedom) is maximal. Canonical ensembles are those in which some macroscopic quantities are externally controlled by a reservoir and particles can occupy individually different levels of such quantity. The most probable distribution in the phase space satisfying the constraints can be obtained (e.g. by applying Lagrange’s Multipliers method). When the constraints are based on additive and separable functions N (MN (αl ,P l ,Ql ,Ω l ) = k=1 M1(αk ,P k ,Qk ,Ω k ) ), the particles take independently the values of (αk , P k , Qk , Ω k ) and these values are identically distributed, the system is ideal. Here the second condition is accepted as true because there are no relationships between the coordinates of the particles that are kept during the driving process. The third condition is postulated. Under these assumptions, it is possible to set up canonical-like ensembles based on additive quantities that are conveniently expressed as a function of the coordinates of the phase space. The mathematical formalism of canonical ensembles for ideal systems can be applied. In particular, the whole partition function is equal to the partition function of a single particle to the power N , ZN (β) = [Z1(β) ]N , and the phase space of a single particle can be used to analyze the distribution function. Therefore the partition function of a canonical-like ensemble is expressed as
Page 6 of 8
Z1(β) =
Eur. Phys. J. E (2014) 37: 81
α
Pα,max
1
π
exp (−M1(α,P,Q ,Ω) β) 0
0
0
·hα(P ) · gα(Q ) · P · Q dP dQ dΩ. 2
(1)
Note that all the internal constraints of the system are implicitly included in hα(P ) and gα(Q ) . If these constraints change (for instance, by including friction or cohesive forces), so do these functions.
Finally, if a global density of states function is used (all configurations are gathered), then the expected partition function for P is expressed as Pmax exp (−P χP ) · h(P ) · P 2 dP, (6) Z1(χP ) = 0
while the statistical distribution of Q is directly ρ(Q ) = g(Q ) · Q .
4 Some ensembles for static and jammed granular matter All the elements necessary for an ensemble theory have been already presented. At this point it is possible to analyze the granular ensembles for static and jammed granular matter. The volume or V-ensemble was the first ensemble conceived for granular media [10]. In this canonical-like ensemble the and separable quantity is the volume additive k V(P l ) = k V(P k ) and the control parameter is the compactivity X. For this ensemble the relationship between the volume and the force-moment state of configurations Vα(P k ) is crucial. In the proposed phase space, the partition function of a single particle is expressed as Pα,max Gα exp (−Vα(P ) /X) · hα(P ) · P 2 dP, Z1(X) = 0
α
(2) where Gα has been defined as 1 gα(Q ) Q dQ . Gα =
(3)
0
In the force-moment ensemble the product of the external stress by the volume of the packing is externally controlled [22,23]. The control parameter is the angoricity (tensorial) Aij or its inverse χij . If the framework of reference is set according to the principal directions of the stress field, then χxy = 0.0. This is due to the fact that for cohesionless systems, Σxy is integrated between − Σxx Σyy and Σxx Σyy , and its average value will be only equal to 0.0 when χxy = 0.0 and exp (−χxy Σxy ) = 1.0. Then, defining χP = (χxx + χyy ) and χQ Ω = (χyy − χxx )/(χyy + χxx ) ∈ [0, 1), in the proposed phase space, the partition function of a single particle is expressed as Z1(χP ,χQ Ω ) = Pα,max 1 α
0
0
π
exp (−P χP [1 + Q cos 2ΩχQ Ω ])
0
·hα(P ) · gα(Q ) · P · Q dP dQ dΩ. 2
(4)
In case of isotropic compression χQ Ω = 0.0 and this expression is simplified as Pα,max Gα exp (−P χP ) · hα(P ) · P 2 dP. Z1(χP ) = α
0
(5)
(7)
Although the full canonical ensemble [19] was proposed before the force-moment ensemble, it is actually a combination of the V- and F-ensembles. It is useful to compare those packings of a granular system in which both the volume and the product of the stress by the volume are being externally controlled at the same time. Another ensemble is based on the elastic potential energy of the packings [18]. For this, relationships between the energy and the coordinates of the phase space have to be established. The quantity Mα(P,Q ,Ω) in eq. (1) is directly the potential energy stored by the static packings, i.e., Eeα(P,Q ) . In the stress ensemble the average stress is controlled but not the volume. As the stress is not an additive quantity, the ensemble is not canonical-like. As it is not a separable quantity, the space of a single particle cannot be used. The Boltzmann’s factor (for the case of isotropic compression) is given by l
k P k − V(P k) · lP l , (8) exp −ξp lV where ξp is the Lagrange’s multiplier that sets the average value of the stress.
5 Numerical simulation Some compression tests have been carried out with a classical molecular dynamics (MD) technique [35]2 to obtain different packings of the same system at the same macroscopic force-moment state. Then their statistical distributions of P and Q were obtained. The two-dimensional system consisted of N = 4900 equal-sized spheres of diameter D = 1.0 and mass m = 1.31. There was no gravity. The initial packings were completely ordered squared lattices of 70 × 70 particles, which were compressed by irregular and randomly generated walls. These were made of particles with the same features and behaved as rigid bodies. All the spheres interacted only through contact, with a viscous Hookean—or linear spring-dashpot—interaction normal to their lines of centers (see [36]). The Hookean constant used was Kn = 2.0 · 105 and the viscoelastic damping parameter γn = 5.5 · 102 . The time of viscoelastic collision of two particles is tcol = 8.5 · 10−3 . The timesteps employed were 1.0 · 10−5 (see [37]). 2
http://lammps.sandia.gov.
Eur. Phys. J. E (2014) 37: 81
Fig. 4. Histograms of P of 180 different packings at the same force-moment state.
An average force was applied to each particle of the irregular walls (perpendicularly to them). It was increased linearly up to F = 5·103 at T = 10.0. Then it was linearly decreased a ten percent at the same rate. The compression protocol was applied to the initial sample with 180 different irregular walls to obtain 180 different final packings at the same total force-moment state (P, Q 0). The histograms of P and Q values of these packings are shown in figs. 4 and 5, respectively3 . According to the theory proposed here, these are explained by eqs. (6) and (7), respectively. On the other hand, the values of the force-moment of three control particles, called A, B and C, were also sampled. Particles A and B were right in the middle of the initial packing and in contact with each other, while particle C was halfway between the particles A and C and the lower boundary. The mean values of P of the control samples were, respectively, P¯ A = 7.72·103 , P¯ B = 8.07·103 and P¯ C = 7.44 · 103 (the average value of P of each of the 180 packings was always around P¯ = 7.76 · 103 ). The stardard A deviations of the control samples were S180 = 2.70 · 103 , B 3 C 3 S180 = 3.15 · 10 and S180 = 1.67 · 10 . The correla2 tion coefficients of control samples were rAB = 0.135 and 2 = 0.003. In order to test whether they originate from rAC the same distribution, a Kruskal-Wallis one-way analysis of variance by ranks [38] was used, concluding that there is no evidence of stochastic dominance between the samples.
6 Discussion The 180 packings had different contacts and forces networks (particles were arranged with different configurations and different values of forces), but there does not
Page 7 of 8
Fig. 5. Histograms of Q = Q/P of 180 different packings at the same force-moment state.
seem to be much difference between their statistical distributions of force-moment. The small differences are due to the fact that there is nothing like Liouville’s theorem that guarantees the conservation of density in phase space. According to the theory proposed here, the most probable statistical distributions of P and Q are explained when appropriate density of states functions are included (eqs. (6) and (7)). However, it is still necessary to prove whether the functions obtained with this procedure are universal or they are linked to the employed protocol or the numerical method. As the obtained packings belong to the force-moment ensemble, the values of P of the control particles A, B and C should be independently and identically distributed under the hypothesis of an ideal system. Since low correlation coefficients between the control samples were found, this condition appears to be met, even though the sample of packings only covers a reduced part of the full space of solutions (the packings are at the same force-moment level but, for instance, particles A and B are always, presumibly, in contact since they were so in the initial sample) whereas the theory relies on the hypothesis that the protocol samples all the possible solutions with the same probability. The mean value of P of each of these particles over the 180 packings is reasonably close to the average value of P over the 4900 particles of each packing and they could originate from the same population. But results are not definitive and to go further on the verification of the hypothesis of ideal system, more simulations and nonparametric statistics analysis are needed.
7 Conclusion
3
The bottom and top of the box are the first and third quartiles, the band inside the box is the median and the ends of the whiskers represent the lowest/highest datum still within 1.5 times the interquartile range of the lower/upper quartile.
A general athermal ensemble theory is presented for slightly deformable granular matter. The approach is based on elaborations and developments of the original theory
Page 8 of 8
firstly introduced by Edwards and collaborators. It considers that the different static packings (described by their contact and forces networks) of a granular system that are compatible to some constraints, can be seen as distributions in a force-moment space. The coordinates of these space (P, Q , Ω) are derived from Cauchy’s tensorial components. When the constraint is such that an additive quantity MN (P k ,Qk ,Ω k ) is controlled, canonical-like ensembles arise. Moreover, if this quantity is of the form MN (P k ,Qk ,Ω k ) = k M1(P k ,Qk ,Ω k ) , and all the possible packings are explored by the protocol, particles are assumed to take independently and identically distributed values of force-moment and configuration. Then the system is ideal and it is possible to use the phase space of a single particle. In order to consider all the degrees of freedom of the problem, density of states functions have to be used. They are defined under the postulate of equal a priori probabilities and under the assumption that the protocol is able to explore all possible packings. If not, the measured functions would be linked to the protocol applied. The functions also take into account the fact that the points of the force-moment space are not uniformly filled when all the possible states on the level of the particles are considered. The functions can be defined globally or separately (according to the class of equilibrium or configuration) and can be measured. These functions are of two kinds: hα(P ) and gα(Q ) , being α the configuration. The hα(P ) functions express the rate of growth of possible states with P (which would be proportional to P 2 if it were a continuum medium). The gα(Q ) functions include all the internal constraints of the system: friction, cohesion, etc. These would change if the system does. On the other hand, the statistical weight of the points in the phase space is fixed by the ensemble. Examples of canonical and ideal ensembles are the volume, the elastic potential energy and the force-moment ensembles. A noncanonical and non-ideal ensemble is the stress ensemble. More research is needed to check whether or not the local states of equilibrium are independently and identically distributed, to measure the density of states functions, to check if they depend on the protocol and to disregard eventual computational effects. The author thanks Rafael Jimenez, Ricardo Brito, Diego Maza and Iker Zuriguel for helpful discussions and sincere encouragement. The author gratefully acknowledges support for his career from Labex Mod´elisation & Exp´erimentation pour la ´ Construction Durable and from Laboratoire Navier (Equipe ´ Multi-Echelle).
References 1. P.G. de Gennes, Physica A 261, 267 (1998). 2. H.J. Herrmann, Physica A 313, 188 (2002).
Eur. Phys. J. E (2014) 37: 81 3. H.M. Jaeger, S.R. Nagel, Rev. Mod. Phys. 68, 1259 (1996). 4. J.D. Bernal, J. Mason, Nature 188, 910 (1960). 5. S. Torquato, T.M. Truskett, P.G. Debenedetti, Phys. Rev. Lett. 84, 2064 (2000). 6. S. Torquato, F.H. Stillinger, Rev. Mod. Phys. 82, 2633 (2010). 7. F. Radjai, M. Jean, J.-J. Moreau, S. Roux, Phys. Rev. Lett. 77, 274 (1996). 8. T.S. Majmudar, R.P. Behringer, Nature 435, 1079 (2005). 9. J.-N. Roux, Phys. Rev. E 61, 6802 (2000). 10. S. Edwards, R. Oakeshott, Physica A 157, 1080 (1989). 11. M.P. Ciamarra, P. Richard, M. Schroter, B.P. Tighe, Soft Matter 8, 9731 (2012). 12. C. Song, P. Wang, Y. Jin, H. Makse, Nature 453, 629 (2008). 13. C. Song, P. Wang, Y. Jin, H. Makse, Physica A 389, 4497 (2010). 14. P. Wang, C. Song, Y. Jin, H. Makse, Physica A 390, 427 (2011). 15. C. Briscoe, C. Song, P. Wang, H. Makse, Physica A 389, 3978 (2010). 16. A. Coniglio, M. Nicodemi, Physica A 296, 451 (2001). 17. A. Fierro, M. Nicodemi, A. Coniglio, Phys. Rev. E 66, 061301 (2002). 18. I.G. Tejada, Physica A 390, 2664 (2011). 19. S. Edwards, Physica A 353, 114 (2005). 20. R. Blumenfeld, S. Edwards, J. Phys. Chem. B 113, 3981 (2009). 21. L.A. Pugnaloni, I. Sanchez, P.A. Gago, J. Damas, I. Zuriguel, D. Maza, Phys. Rev. E 82, 050301 (2010). 22. S. Henkes, C.S.O. Hern, B. Chakraborty, Phys. Rev. Lett. 99, 038002 (2007). 23. S. Henkes, B. Chakraborty, Phys. Rev. E 79, 061301 (2009). 24. J.H. Snoeijer, T.J.H. Vlugt, M. van Hecke, J. van Saarloos, Phys. Rev. Lett. 92, 054302 (2004). 25. J.H. Snoeijer, T.J.H. Vlugt, W.G. Ellenbroek, M. van Hecke, J.M.J. van Leeuwen, Phys. Rev. E 70, 061306 (2004). 26. J.H. Snoeijer, W.G. Ellenbroek, T.J.H. Vlugt, M. van Hecke, Phys. Rev. Lett. 96, 098001 (2006). 27. B.P. Tighe, A.R.T. van Eerd, T.J.H. Vlugt, Phys. Rev. Lett. 100, 238001 (2008). 28. K. Johnson, Contact Mechanics (Cambridge University Press, 1992). 29. A. Liu, S.R. Nagel, Nature 396, 21 (1998). 30. K. Bagi, Mech. Mater. 22, 165 (1996). 31. M.P. Ciamarra, A. Coniglio, M. Nicodemi, Phys. Rev. Lett. 97, 158001 (2006). 32. J.B. Knight, C.G. Fandrich, C.N. Lau, H.M. Jaeger, S.R. Nagel, Phys. Rev. E 51, 3957 (1995). 33. E.R. Nowak, J.B. Knight, E. Ben-Naim, H.M. Jaeger, S.R. Nagel, Phys. Rev. E 57, 1971 (1998). 34. M. Schroter, D.I. Goldman, H.L. Swinney, Phys. Rev. E 71, 030301 (2005). 35. S.J. Plimpton, J. Comput. Phys. 117, 1 (1995). 36. S. Luding, Collisions and Contacts between Two Particles, in Physics of Dry Granular Media, edited by H. Herrmann, J.-P. Hovi, S. Luding, Vol. 350 of NATO ASI Series (Springer, The Netherlands, 1998) pp. 285–304. 37. I. Tejada, R. Jimenez, Eur. Phys. J. E 37, 15 (2014). 38. W.H. Kruskal, W.A. Wallis, J. Am. Stat. Assoc. 47, 583 (1952).
THE EUROPEAN PHYSICAL JOURNAL E
Regular Article
Ensemble theory for slightly deformable granular matter Ignacio G. Tejadaa Universit´e Paris-Est, Laboratoire Navier (UMR 8205), CNRS, ENPC, IFSTTAR, France Received 19 May 2014 and Received in final form 14 August 2014 c EDP Sciences / Societ` Published online: 25 September 2014 – a Italiana di Fisica / Springer-Verlag 2014 Abstract. Given a granular system of slightly deformable particles, it is possible to obtain different static and jammed packings subjected to the same macroscopic constraints. These microstates can be compared in a mathematical space defined by the components of the force-moment tensor (i.e. the product of the equivalent stress by the volume of the Voronoi cell). In order to explain the statistical distributions observed there, an athermal ensemble theory can be used. This work proposes a formalism (based on developments of the original theory of Edwards and collaborators) that considers both the internal and the external constraints of the problem. The former give the density of states of the points of this space, and the latter give their statistical weight. The internal constraints are those caused by the intrinsic features of the system (e.g. size distribution, friction, cohesion). They, together with the force-balance condition, determine which the possible local states of equilibrium of a particle are. Under the principle of equal a priori probabilities, and when no other constraints are imposed, it can be assumed that particles are equally likely to be found in any one of these local states of equilibrium. Then a flat sampling over all these local states turns into a non-uniform distribution in the force-moment space that can be represented with density of states functions. Although these functions can be measured, some of their features are explored in this paper. The external constraints are those macroscopic quantities that define the ensemble and are fixed by the protocol. The force-moment, the volume, the elastic potential energy and the stress are some examples of quantities that can be expressed as functions of the force-moment. The associated ensembles are included in the formalism presented here.
1 Introduction Granular materials are large configurations of discrete macroscopic particles of size larger than one micron [1] that display unusual behavior in comparison with that of solids, liquids or gases [1–3]. Since they are present in many different aspects of our daily lives, the establishment of a theoretical framework for their description has become an important issue. There are two important aspects that contribute to their unique properties: ordinary temperature plays no role, and interactions between grains are dissipative because of static friction and the inelasticity of collisions. One of the consequences of this is the possibility of obtaining different static packings of the same granular system under the same macroscopic constraints. If we focus on these packings, we can find different and inhomogeneous distributions of contacts, interparticle forces, voids, etc. [4–8]. Moreover, the mechanical behaviors of the packings are determined by their geometries, which are in turn adapted to the applied loads [9]. a
e-mail: [email protected]
Since granular systems consist of a lot of similar particles in which the interaction laws are rather known, a statistical mechanics approach might be feasible to describe them [10,11]. Edwards and Oakeshott were the first to do this [10]. Assuming that granular systems have entropy, they claimed that the volume plays the role of energy (V-ensemble [10,12–15]). There are other formalisms based on the energy of the whole system, such as [16–18]. It has been also suggested to consider the stress together with the volume (V–F-ensemble or full canonical ensemble) [19–21]. If the constraint of the volume is removed from the last, then the F-ensemble arises [22,23]. The force network ensemble theory [24–27] focuses on the statistical distribution of interparticle forces with a given contacts network. The aim of this paper is to propose a formalism for a general athermal ensemble theory that is valid to compare static packings of slightly deformable granular matter. Under some assumptions, a mathematical space is used to compare microstates and to obtain which the most probable distribution of some quantities is when some external constraints are imposed. Some of the features of this
Page 2 of 8
space are analyzed. Within this framework, some of the ensembles reported above, as well as the stress ensemble, can be analyzed.
2 Basic assumptions and definitions 2.1 Granular system, packings and protocol A granular system is a set of N discrete macroscopic particles that is determined by their features. For instance, for the systems analyzed here, these features are: the size and shape, the roughness and the mechanical properties of the material (e.g. density, constitutive relationships, strength criteria, etc.). For simplicity, without loss of generality, the formalism included in this paper is adapted for twodimensional systems, but the theoretical scheme is valid for three-dimensional systems too. The motion of each particle under the action of external and interaction forces is governed by classical mechanics. Only mechanical pairwise interactions are considered here. These include dissipative or viscous terms to represent inelastic collisions. Particles are deformable when they change their shape with boundary and body forces. Only small and elastic deformations are expected here. The contact theory (see [28]) can be used to find the interaction forces. In the case of two-dimensional elastic particles, k and l, the proposed relationship between the force (per unit of length) and the associated overlapping is linear (Hookean-like force): Fkl = −Kn δ kl xkl /xkl , with xkl = xl −xk and δ kl = D −xkl (for xkl ≤ D). The normal viscous force is proportional to γn (m/2)δ˙ kl , being γn the viscoelastic damping constant. Each static arrangement of a given granular system is referred to as packing. Static means that particles are force-balanced and at rest, so that state corresponds to a local energy minimum in the rugged potential energy landscape produced by the pairwise interactions. The packings may be in turn defined by their contacts and forces networks. Several forces networks may be associated with the same contacts network when the boundary forces are not fixed or when the contacts network is hyperstatic, which means that the number of contact points is higher than the strictly condition of static equilibrium (isotaticity). A protocol is a realization of different static packings of a system in a systematic manner (e.g. pouring, shaking, compressing, etc.). Although the obtained packings are static, it is actually a dynamic process in which kinetic energy is injected and dissipated. If microscopic degrees of freedom (e.g. phonons, electrons, etc.) were included in the model, then the energy dissipation would actually be a redistribution of grain-level energy to these microscopic degrees of freedom. This approach is thought for jammed states, although a precise definition of jamming is not attempted now. These states can be seen as those configurations in which each particle is in contact with its nearest neighbors in such a way that mechanical stability of a specific type is conferred to the packing [6]. Considering the jamming
Eur. Phys. J. E (2014) 37: 81
phase space [29], jammed packings of a granular system are close to the origin of the temperature-volume-stress space. However the limits of validity of this approach are not established yet. For practical purposes, all the packings analyzed here are highly hyperstatic. For each particle k in a packing Π, the corresponding region consisting of all points closer to its center than to the center of any other is itsso-called Voronoi cell. Its associated volume is V k , and k V k = V . In static packings, the equivalent stress field of each particle (the stress field averaged over its Voronoi cell) can be obtained from the mechanical forces that keep it k = (1/V k ) l Fikl rjkl where rikl is the in equilibrium: σij spatial position of the contact point between particles k k k and l (see [30]). The product Σij = σij · V k = l Fikl rjkl is referred to as force-moment or extensive pressure. The total force-moment and the average k stress of the whole system are given by Σij = k Σij and σij = Σij /V , respectively. Instead of the Cauchy components of the force-moment tensor Σij , it is interesting to use other variables that are functions of its invariants: I1 = Σxx + Σyy and 2 . The mean force-moment component I2 = Σxx Σyy − Σxy is defined as P = I1 /2 and expresses the level of the field. I 2 1 − I2 and the The deviatoric component is Q = 2 orientation Ω =
1 2
2Σ
xy arctan Σyy −Σ . P is the semi-sum of xx
principal force-moment values, i.e., the eigenvalues, and Q is equal to their semi-difference. The third component Ω ∈ [0, π) gives precisely the orientation of the field with respect to the external framework of reference, so it is not invariant. The relative deviatoric force-moment is a dimensionless quantity defined as Q = Q/P that gives an idea of how isotropic Q → 0.0 or anisotropic Q → 1.0 the field is. 2.2 Configurations On the level of the particles, many possible states of equilibrium exist (fig. 1). These states can be classified according to the number of forces, their relative angles, their value and direction, and the orientation of the set with respect to an external framework of reference. The term configuration is used here for a specific number of forces with a specific layout. Note that the orientation has been excluded in this definition (e.g. the states depicted in fig. 1C belong to the same configuration). It is worth noting that the same configuration α can often represent different states of equilibrium with different sets forces (fig. 2). Thus, it may represent different states of force-moment too. In the end, the possible configurations of a system are determined by the condition of force-balance l Fkl = 0 and by the internal constraints, i.e. those which are exclusively dependent on the features of the particles. If the system does not change, neither do its internal constraints nor do its possible configurations. Therefore the possible configurations of a system do not depend on the external
Eur. Phys. J. E (2014) 37: 81
Page 3 of 8
Fig. 2. Two possible local states of static equilibrium for a particle k that belong to the same configuration but have different sets of forces and represent different states of force-moment.
Fig. 3. An example of the same configuration with two different Voronoi cells. The second is only possible if the system is polydisperse. Fig. 1. Some examples of possible local states of static equilibrium for a particle k. A: With a different number of forces. B: With the same number of forces but a different layout. C: With the same layout but a different orientation.
constraints imposed to the packings that are being sampled. Some examples of internal constraints are: – In case of monodisperse systems (particles equally sized) the relative angle between two forces Δθkl = |θkl+1 − θkl | cannot be smaller than π/3. – In case of frictionless particles, all the forces must be central Fkl × xkl = 0. – In case of cohesionless materials all forces must be compressive. Fkl · xkl < 0.0. In configurations with a low number of forces (and especially when the system is polydisperse), it is possible to find local packings with the same layout but different Voronoi cell and thus different volumes (see fig. 3). However this effect is disregarded here. For simplicity, and assuming that packings are tight, a unique value of the volume Vα is supposed to be associated with each configuration1 . Nevertheless, in case of slightly deformable particles, this volume decreases as the forces are increased. 1 Otherwise some density of states functions should be introduced to consider all possible values between Vα,min and Vα,max .
It depends basically on stiffness (included in the contact law), and it is reasonable to assume that it only does with P (the variable that sets the force-moment level). Therefore Vα = Vα(P ) .
3 Elements of ensemble theory A general athermal ensemble theory is proposed for compressed, static and jammed packings of a given granular system of N slightly deformable particles that satisfy some macroscopic constraints. If N is an extremely large number, then it is possible to perform an analysis in the so-called thermodynamic limit. The specification of the number of particles and the actual values of the macroscopic quantities of a packing define its macrostate. Some macroscopic quantities, symbolized in general as M , are: the volume of the packing, V , the macroscopic stress, σij , the elastic potential energy that the packing stores, Ee , or the force-moment Σij = V σij . At grain level, however, a large number of possibilities still exist because at that level there will be in general a large number of different ways, W, in which the macrostate {M } can be realized (note that microscopic internal degrees of freedom are not included). Each of these different ways, i.e. each static packing (contacts and forces networks), specifies a microstate of the system. An ensemble E(N,M ) is precisely defined as a large number of copies
Page 4 of 8
of a system, each of which represents a possible microstate that the real system might be in or a large number of different packings of the same system satisfying the same macroscopic constraints. The exploration of microstates is not produced as time goes by (as it happens for thermal systems), because the energy scale kT is insignificant for granular media (at grain level) and because of the inelasticity of collisions and the non-conservative friction forces. Therefore a protocol is needed to explore the microstates. Identical systems are assumed to give the same average response to the same driving protocol (e.g. [31–34]). The postulate of equal a priori probabilities states that for a given macrostate of the system, and when there are no other constraints, after an appropriate protocol the system is equally likely to be in any one of these microstates. Appropriate protocol means precisely that both the forces and the contacts networks are constantly being changed and the system explores all its possible packings. It is assumed that microscopic internal degrees of freedom do not add new constraints (i.e. neither exclude microstates nor alter their statistical weights). In such case, the actual number of microstates W is just a function of the macroscopic variables: W = W(N,M ) . However, instead of trying to figure out what these whole packings are, we can focus on the local state of the particles. Then the following assumption is made: when all the possible packings of an ensemble, Π ∈ E(N,M ) , are considered, each particle k can be found in any equilibrium state that corresponds to a possible configuration α of the system, and there is not any a priori reason to favor any of these states more than any other. Furthermore, if the contact network were not changed during the driving process, it would be possible to establish some relationships between the values of the variables of the particles. Then the number of degrees of freedom of the packings would be equal to the number of free variables (those that can be modified to satisfy the macroscopic constraints). However, this is not the case. Here, all the possible packings of a system are being considered and, in consequence, there is no reason to establish any fixed relation between, for example, positions, neighbors or values of forces. In other words, for each contact network there are relationships that result on a known number of degrees of freedom (which depends on the hyperstaticity of the network and the variability of boundary forces), but when all of them are considered, particles are supposed to behave independently. There are no fixed relationships during the exploration of microstates. 3.1 Phase space The phase space is a mathematical space in which all the possible states of the system can be represented. Therefore the microstates can be compared there. Here the proposed coordinates are the configuration αk and the forcemoment state (P k , Qk , Ω k ) of each particle k. In this space some macroscopic constraints may be expressed as a function of these coordinates: M = M(αk ;P k ,Qk ,Ω k ) .
Eur. Phys. J. E (2014) 37: 81
Thus, this space is discrete in configurations and continuous in components of the force-moment tensor. The integration limits of the force-moment coordinates are those which guarantee that the stress field is compressive in any orientation. It can be proved that it happens when k k k k Σ k or, equivalently, Σxx , Σyy ∈ [0, ∞) and |Σxy | ≤ Σxx yy when P k ∈ [0, ∞), Qk ∈ [0, 1] and Ω k ∈ [0, π). In this space all the microstates of the system can be represented. For instance, it allows to locate particles that are locally packed at the same volume but are arranged according to different geometrical patterns, or particles that, although their arrangement follows the same geometrical pattern, have different values of the forces (and, consequently, of the force-moment). However, with this choice of coordinates, some information has been lost (e.g. the same configuration and force-moment state can be sometimes represented by different sets of values of the forces). It means that a flat sampling over force-moment coordinates is not taking into account all the possible states of equilibrium. This can be solved by using density of states functions gα(P,Q ,Ω) that are actually the result of a flat sampling over all possible local states of equilibrium. All the degrees of freedom are being considered within these functions. In particular, in two-dimensional problems the use of three components for the force-moment (P k , Qk , Ω k ) and one for the type of configuration αk does not absolutely mean that only 4N degrees of freedom are being considered. Although the establishment ab initio of gα(P,Q ,Ω) functions is complicated, this is not the objective, because they can be measured. Moreover, the definition of configurations {α1 , α2 , . . . , αr } proposed here is arbitrary. For example, all the configurations having the same number of forces could be gathered in an upper class (A1 = {α1 , . . . , αm }, A2 = {αm+1 , . . . , αn }, etc.) that does not distinguish between different layouts. Then there would be another classification {A1 , A2 , . . . , AZ } with different functions gA(P,Q ,Ω) . Even configurations may not be distinguished at all, but just using a global density of states function over the phase space g(P,Q ,Ω) . Using more configurations means that there are more functions to measure. However it allows to get more information about the most probable distributions (e.g. expected ratios between particles arranged as honeycombs with respect to those arranged with 5 forces). Although these functions still remain unknown, some symmetries can be found. For instance, consider the case of a monodisperse, frictionless and cohesionless system, and a particle k that is in a specific state of equilibrium under a set of forces {Fkl } and has an associated Voronoi cell of volume V k . It represent a force-moment state (P, Q , Ω). It can be easily proved that such state of equilibrium could also represent a force-moment state (P, Q , Ω + Θ), for any arbitrary Θ, just by rotating its layout an angle Θ. As configurations are defined independently of their spatial orientation, the density of states functions do not depend on the orientation of the field, so gα = gα(Ω) .
Eur. Phys. J. E (2014) 37: 81
On the other hand, if all the forces acting on the particle k are scaled with any arbitrary factor x, the particle remains in equilibrium (because the force-balance and the internal constraints are not affected by this scale transformation). If the variation of V k is disregarded, then the force-moment state changes from (P, Q , Ω) to (xP, Q , Ω). Therefore the possible number of states at level xP is at least the same as at level P . Indeed it usually grows with P , although the rate is uncertain. While a flat sampling over the force-moment variables Σxx , Σxy , Σyy leads to a growth that is proportional to P 2 (since dΣxx dΣxy dΣyy = 4P 2 Q dP dQ dΩ), a flat sampling over the values of the forces of each configuration does not. To illustrate this, consider a honeycomb configuration as the one shown in fig. 2. The equilibrium is produced by 6 forces {F k1 , . . . , F k6 } which follow a continuous uniform distribution U(0,Fmax ) . The minimum value is 0 because the system is cohesionless, while a maximum value is necessary because of the contact law (although usually F kl Fmax ). If the force-balance condition is imposed, then there are two additional constraints that make it possible to select randomly just four forces, say F k1 , F k2 , F k3 , F k4 ∈ U(0,Fmax ) , and compute k6 the other two, F k5 = F k1 + F k2 − F k4 and F = k3 k4 k1 k F + F − F . On the other hand, P = (R /2) l F kl or P = (Rk /2)(F k1 +2F k2 +2F k3 +F k4 ). Therefore, when a flat sampling over the space of forces is done, the statistical distribution of P for the honeycomb configuration results from the addition of independent random variables uniformly distributed. When F kl Fmax , as usual, i.e. for lower values of P , this distribution is proportional to P 3 . By contrast, in the case of a square lattice the number of states grows linearly with P (at low levels). And the same could be done with all possible configurations to compare them and to normalize the density of states functions. For those with five forces, it grows with P 2 , while the number of states is independent of P for configurations with three and two (opposite) forces. However, as the intention is to measure the functions, the important conclusion to extract is that the number of states grows monotonically with P and at rates that are different to P 2 . Furthermore, the distribution over Q values coming from a flat sampling over the space of forces is much more complex, but also feasible. This is, in the end, what makes a configuration more or less probable than others (when they are normalized by considering all possible configurations) and here is where the largest differences appear. For instance, in two-dimensional systems, the honeycomb represents by far more states than any other configuration. Under the approach presented here, this is the reason why these systems tend to crystallize. Considering all this, it seems convenient to express gα(Σij ) = hα(P ) · gα(Q ) . Under this scheme, gα(Q ) takes into account the internal constraints of equilibrium of each configuration while hα(P ) is the correction needed to represent properly the growth of the number of states with P .
Page 5 of 8
3.2 Canonical-like ensembles and equilibrium When a protocol is applied to a granular system, particles change their coordinates in the force-moment space with the cycles during the driving process, but it is not established how it happens (as, for example, Hamiltonian mechanics does in the phase space for classical ideal gases). Nevertheless the macroscopic constraints of the ensemble at the end of the driving cycles are fixed by the protocol. Then, the (thermodynamical) equilibrium is defined by reference to the macrostate (Boltzmann’s approach), which is uniquely given by a set of values of the macrovariables. The expected distribution in the equilibrium is obtained by maximizing Boltzmann’s entropy S ∝ ln W(N,M ) , i.e., by analyzing which the most probable distribution of an ensemble subjected to the same constraints is. Therefore, the problem is not mechanical but statistical, because the role of the dynamics during driving cycles is trivial for this definition. Moreover, since the exploration of mechanical stable states is not produced as time goes by (as thermal systems do), the ergodic hypothesis is not based in time averages but in ensemble averages (independent realizations of same process). Unfortunately, there is nothing like Liouville’s theorem to know how the evolution of the statistical distribution is when the packings are realized. This theory just gives the most probable distribution, but it does not necessarily remain constant with driving. It is worth mentioning that if the microscopic degrees of freedom were included, the total energy of all would be conserved and redistributed and that force equilibrium would happen when the total entropy (considering all the degrees of freedom) is maximal. Canonical ensembles are those in which some macroscopic quantities are externally controlled by a reservoir and particles can occupy individually different levels of such quantity. The most probable distribution in the phase space satisfying the constraints can be obtained (e.g. by applying Lagrange’s Multipliers method). When the constraints are based on additive and separable functions N (MN (αl ,P l ,Ql ,Ω l ) = k=1 M1(αk ,P k ,Qk ,Ω k ) ), the particles take independently the values of (αk , P k , Qk , Ω k ) and these values are identically distributed, the system is ideal. Here the second condition is accepted as true because there are no relationships between the coordinates of the particles that are kept during the driving process. The third condition is postulated. Under these assumptions, it is possible to set up canonical-like ensembles based on additive quantities that are conveniently expressed as a function of the coordinates of the phase space. The mathematical formalism of canonical ensembles for ideal systems can be applied. In particular, the whole partition function is equal to the partition function of a single particle to the power N , ZN (β) = [Z1(β) ]N , and the phase space of a single particle can be used to analyze the distribution function. Therefore the partition function of a canonical-like ensemble is expressed as
Page 6 of 8
Z1(β) =
Eur. Phys. J. E (2014) 37: 81
α
Pα,max
1
π
exp (−M1(α,P,Q ,Ω) β) 0
0
0
·hα(P ) · gα(Q ) · P · Q dP dQ dΩ. 2
(1)
Note that all the internal constraints of the system are implicitly included in hα(P ) and gα(Q ) . If these constraints change (for instance, by including friction or cohesive forces), so do these functions.
Finally, if a global density of states function is used (all configurations are gathered), then the expected partition function for P is expressed as Pmax exp (−P χP ) · h(P ) · P 2 dP, (6) Z1(χP ) = 0
while the statistical distribution of Q is directly ρ(Q ) = g(Q ) · Q .
4 Some ensembles for static and jammed granular matter All the elements necessary for an ensemble theory have been already presented. At this point it is possible to analyze the granular ensembles for static and jammed granular matter. The volume or V-ensemble was the first ensemble conceived for granular media [10]. In this canonical-like ensemble the and separable quantity is the volume additive k V(P l ) = k V(P k ) and the control parameter is the compactivity X. For this ensemble the relationship between the volume and the force-moment state of configurations Vα(P k ) is crucial. In the proposed phase space, the partition function of a single particle is expressed as Pα,max Gα exp (−Vα(P ) /X) · hα(P ) · P 2 dP, Z1(X) = 0
α
(2) where Gα has been defined as 1 gα(Q ) Q dQ . Gα =
(3)
0
In the force-moment ensemble the product of the external stress by the volume of the packing is externally controlled [22,23]. The control parameter is the angoricity (tensorial) Aij or its inverse χij . If the framework of reference is set according to the principal directions of the stress field, then χxy = 0.0. This is due to the fact that for cohesionless systems, Σxy is integrated between − Σxx Σyy and Σxx Σyy , and its average value will be only equal to 0.0 when χxy = 0.0 and exp (−χxy Σxy ) = 1.0. Then, defining χP = (χxx + χyy ) and χQ Ω = (χyy − χxx )/(χyy + χxx ) ∈ [0, 1), in the proposed phase space, the partition function of a single particle is expressed as Z1(χP ,χQ Ω ) = Pα,max 1 α
0
0
π
exp (−P χP [1 + Q cos 2ΩχQ Ω ])
0
·hα(P ) · gα(Q ) · P · Q dP dQ dΩ. 2
(4)
In case of isotropic compression χQ Ω = 0.0 and this expression is simplified as Pα,max Gα exp (−P χP ) · hα(P ) · P 2 dP. Z1(χP ) = α
0
(5)
(7)
Although the full canonical ensemble [19] was proposed before the force-moment ensemble, it is actually a combination of the V- and F-ensembles. It is useful to compare those packings of a granular system in which both the volume and the product of the stress by the volume are being externally controlled at the same time. Another ensemble is based on the elastic potential energy of the packings [18]. For this, relationships between the energy and the coordinates of the phase space have to be established. The quantity Mα(P,Q ,Ω) in eq. (1) is directly the potential energy stored by the static packings, i.e., Eeα(P,Q ) . In the stress ensemble the average stress is controlled but not the volume. As the stress is not an additive quantity, the ensemble is not canonical-like. As it is not a separable quantity, the space of a single particle cannot be used. The Boltzmann’s factor (for the case of isotropic compression) is given by l
k P k − V(P k) · lP l , (8) exp −ξp lV where ξp is the Lagrange’s multiplier that sets the average value of the stress.
5 Numerical simulation Some compression tests have been carried out with a classical molecular dynamics (MD) technique [35]2 to obtain different packings of the same system at the same macroscopic force-moment state. Then their statistical distributions of P and Q were obtained. The two-dimensional system consisted of N = 4900 equal-sized spheres of diameter D = 1.0 and mass m = 1.31. There was no gravity. The initial packings were completely ordered squared lattices of 70 × 70 particles, which were compressed by irregular and randomly generated walls. These were made of particles with the same features and behaved as rigid bodies. All the spheres interacted only through contact, with a viscous Hookean—or linear spring-dashpot—interaction normal to their lines of centers (see [36]). The Hookean constant used was Kn = 2.0 · 105 and the viscoelastic damping parameter γn = 5.5 · 102 . The time of viscoelastic collision of two particles is tcol = 8.5 · 10−3 . The timesteps employed were 1.0 · 10−5 (see [37]). 2
http://lammps.sandia.gov.
Eur. Phys. J. E (2014) 37: 81
Fig. 4. Histograms of P of 180 different packings at the same force-moment state.
An average force was applied to each particle of the irregular walls (perpendicularly to them). It was increased linearly up to F = 5·103 at T = 10.0. Then it was linearly decreased a ten percent at the same rate. The compression protocol was applied to the initial sample with 180 different irregular walls to obtain 180 different final packings at the same total force-moment state (P, Q 0). The histograms of P and Q values of these packings are shown in figs. 4 and 5, respectively3 . According to the theory proposed here, these are explained by eqs. (6) and (7), respectively. On the other hand, the values of the force-moment of three control particles, called A, B and C, were also sampled. Particles A and B were right in the middle of the initial packing and in contact with each other, while particle C was halfway between the particles A and C and the lower boundary. The mean values of P of the control samples were, respectively, P¯ A = 7.72·103 , P¯ B = 8.07·103 and P¯ C = 7.44 · 103 (the average value of P of each of the 180 packings was always around P¯ = 7.76 · 103 ). The stardard A deviations of the control samples were S180 = 2.70 · 103 , B 3 C 3 S180 = 3.15 · 10 and S180 = 1.67 · 10 . The correla2 tion coefficients of control samples were rAB = 0.135 and 2 = 0.003. In order to test whether they originate from rAC the same distribution, a Kruskal-Wallis one-way analysis of variance by ranks [38] was used, concluding that there is no evidence of stochastic dominance between the samples.
6 Discussion The 180 packings had different contacts and forces networks (particles were arranged with different configurations and different values of forces), but there does not
Page 7 of 8
Fig. 5. Histograms of Q = Q/P of 180 different packings at the same force-moment state.
seem to be much difference between their statistical distributions of force-moment. The small differences are due to the fact that there is nothing like Liouville’s theorem that guarantees the conservation of density in phase space. According to the theory proposed here, the most probable statistical distributions of P and Q are explained when appropriate density of states functions are included (eqs. (6) and (7)). However, it is still necessary to prove whether the functions obtained with this procedure are universal or they are linked to the employed protocol or the numerical method. As the obtained packings belong to the force-moment ensemble, the values of P of the control particles A, B and C should be independently and identically distributed under the hypothesis of an ideal system. Since low correlation coefficients between the control samples were found, this condition appears to be met, even though the sample of packings only covers a reduced part of the full space of solutions (the packings are at the same force-moment level but, for instance, particles A and B are always, presumibly, in contact since they were so in the initial sample) whereas the theory relies on the hypothesis that the protocol samples all the possible solutions with the same probability. The mean value of P of each of these particles over the 180 packings is reasonably close to the average value of P over the 4900 particles of each packing and they could originate from the same population. But results are not definitive and to go further on the verification of the hypothesis of ideal system, more simulations and nonparametric statistics analysis are needed.
7 Conclusion
3
The bottom and top of the box are the first and third quartiles, the band inside the box is the median and the ends of the whiskers represent the lowest/highest datum still within 1.5 times the interquartile range of the lower/upper quartile.
A general athermal ensemble theory is presented for slightly deformable granular matter. The approach is based on elaborations and developments of the original theory
Page 8 of 8
firstly introduced by Edwards and collaborators. It considers that the different static packings (described by their contact and forces networks) of a granular system that are compatible to some constraints, can be seen as distributions in a force-moment space. The coordinates of these space (P, Q , Ω) are derived from Cauchy’s tensorial components. When the constraint is such that an additive quantity MN (P k ,Qk ,Ω k ) is controlled, canonical-like ensembles arise. Moreover, if this quantity is of the form MN (P k ,Qk ,Ω k ) = k M1(P k ,Qk ,Ω k ) , and all the possible packings are explored by the protocol, particles are assumed to take independently and identically distributed values of force-moment and configuration. Then the system is ideal and it is possible to use the phase space of a single particle. In order to consider all the degrees of freedom of the problem, density of states functions have to be used. They are defined under the postulate of equal a priori probabilities and under the assumption that the protocol is able to explore all possible packings. If not, the measured functions would be linked to the protocol applied. The functions also take into account the fact that the points of the force-moment space are not uniformly filled when all the possible states on the level of the particles are considered. The functions can be defined globally or separately (according to the class of equilibrium or configuration) and can be measured. These functions are of two kinds: hα(P ) and gα(Q ) , being α the configuration. The hα(P ) functions express the rate of growth of possible states with P (which would be proportional to P 2 if it were a continuum medium). The gα(Q ) functions include all the internal constraints of the system: friction, cohesion, etc. These would change if the system does. On the other hand, the statistical weight of the points in the phase space is fixed by the ensemble. Examples of canonical and ideal ensembles are the volume, the elastic potential energy and the force-moment ensembles. A noncanonical and non-ideal ensemble is the stress ensemble. More research is needed to check whether or not the local states of equilibrium are independently and identically distributed, to measure the density of states functions, to check if they depend on the protocol and to disregard eventual computational effects. The author thanks Rafael Jimenez, Ricardo Brito, Diego Maza and Iker Zuriguel for helpful discussions and sincere encouragement. The author gratefully acknowledges support for his career from Labex Mod´elisation & Exp´erimentation pour la ´ Construction Durable and from Laboratoire Navier (Equipe ´ Multi-Echelle).
References 1. P.G. de Gennes, Physica A 261, 267 (1998). 2. H.J. Herrmann, Physica A 313, 188 (2002).
Eur. Phys. J. E (2014) 37: 81 3. H.M. Jaeger, S.R. Nagel, Rev. Mod. Phys. 68, 1259 (1996). 4. J.D. Bernal, J. Mason, Nature 188, 910 (1960). 5. S. Torquato, T.M. Truskett, P.G. Debenedetti, Phys. Rev. Lett. 84, 2064 (2000). 6. S. Torquato, F.H. Stillinger, Rev. Mod. Phys. 82, 2633 (2010). 7. F. Radjai, M. Jean, J.-J. Moreau, S. Roux, Phys. Rev. Lett. 77, 274 (1996). 8. T.S. Majmudar, R.P. Behringer, Nature 435, 1079 (2005). 9. J.-N. Roux, Phys. Rev. E 61, 6802 (2000). 10. S. Edwards, R. Oakeshott, Physica A 157, 1080 (1989). 11. M.P. Ciamarra, P. Richard, M. Schroter, B.P. Tighe, Soft Matter 8, 9731 (2012). 12. C. Song, P. Wang, Y. Jin, H. Makse, Nature 453, 629 (2008). 13. C. Song, P. Wang, Y. Jin, H. Makse, Physica A 389, 4497 (2010). 14. P. Wang, C. Song, Y. Jin, H. Makse, Physica A 390, 427 (2011). 15. C. Briscoe, C. Song, P. Wang, H. Makse, Physica A 389, 3978 (2010). 16. A. Coniglio, M. Nicodemi, Physica A 296, 451 (2001). 17. A. Fierro, M. Nicodemi, A. Coniglio, Phys. Rev. E 66, 061301 (2002). 18. I.G. Tejada, Physica A 390, 2664 (2011). 19. S. Edwards, Physica A 353, 114 (2005). 20. R. Blumenfeld, S. Edwards, J. Phys. Chem. B 113, 3981 (2009). 21. L.A. Pugnaloni, I. Sanchez, P.A. Gago, J. Damas, I. Zuriguel, D. Maza, Phys. Rev. E 82, 050301 (2010). 22. S. Henkes, C.S.O. Hern, B. Chakraborty, Phys. Rev. Lett. 99, 038002 (2007). 23. S. Henkes, B. Chakraborty, Phys. Rev. E 79, 061301 (2009). 24. J.H. Snoeijer, T.J.H. Vlugt, M. van Hecke, J. van Saarloos, Phys. Rev. Lett. 92, 054302 (2004). 25. J.H. Snoeijer, T.J.H. Vlugt, W.G. Ellenbroek, M. van Hecke, J.M.J. van Leeuwen, Phys. Rev. E 70, 061306 (2004). 26. J.H. Snoeijer, W.G. Ellenbroek, T.J.H. Vlugt, M. van Hecke, Phys. Rev. Lett. 96, 098001 (2006). 27. B.P. Tighe, A.R.T. van Eerd, T.J.H. Vlugt, Phys. Rev. Lett. 100, 238001 (2008). 28. K. Johnson, Contact Mechanics (Cambridge University Press, 1992). 29. A. Liu, S.R. Nagel, Nature 396, 21 (1998). 30. K. Bagi, Mech. Mater. 22, 165 (1996). 31. M.P. Ciamarra, A. Coniglio, M. Nicodemi, Phys. Rev. Lett. 97, 158001 (2006). 32. J.B. Knight, C.G. Fandrich, C.N. Lau, H.M. Jaeger, S.R. Nagel, Phys. Rev. E 51, 3957 (1995). 33. E.R. Nowak, J.B. Knight, E. Ben-Naim, H.M. Jaeger, S.R. Nagel, Phys. Rev. E 57, 1971 (1998). 34. M. Schroter, D.I. Goldman, H.L. Swinney, Phys. Rev. E 71, 030301 (2005). 35. S.J. Plimpton, J. Comput. Phys. 117, 1 (1995). 36. S. Luding, Collisions and Contacts between Two Particles, in Physics of Dry Granular Media, edited by H. Herrmann, J.-P. Hovi, S. Luding, Vol. 350 of NATO ASI Series (Springer, The Netherlands, 1998) pp. 285–304. 37. I. Tejada, R. Jimenez, Eur. Phys. J. E 37, 15 (2014). 38. W.H. Kruskal, W.A. Wallis, J. Am. Stat. Assoc. 47, 583 (1952).
