Contact forces in a granular packing

We present the results of a systematic numerical investigation of force distributions in granular packings. We ﬁnd that all the main features of force transmission previously established for two-dimensional systems of hard particles hold in three-dimensional systems and for soft particles, too. In particular, the probability distribution of normal forces falls off exponentially for forces above the mean force. For forces below the mean, this distribution is either a decreasing power law when the system is far from static equilibrium, or nearly uniform at static equilibrium, in agreement with recent experiments. Moreover, we show that the forces below the mean do not contribute to the shear stress. The subnetwork of the contacts carrying a force below the mean thus plays a role similar to a ﬂuid surrounding the solid backbone composed of the contacts carrying a force above the mean. We address the issue of the computation of contact forces in a packing at static equilibrium. We introduce a model with no local simplifying force rules, that allows for an exact computation of contact forces for given granular texture and boundary conditions.


I. INTRODUCTION
Particle motion has been at the focus of experimental observation more often than contact forces. 1 Examples are hopper flow, granular avalanches, and size segregation. Moreover, purely kinematic models are naturally favored in most communities interested in granular materials. However, the crucial role played by the contact network and its inherent correlations does not appear in a straightforward way in the velocity field. In this respect, the inhomogeneities observed in the contact forces by means of the photoelastic method [2][3][4] or simulations similar to the one displayed in Fig.  1 are quite significant. They provide a faithful and direct image of the correlations inside the underlying geometry both at static equilibrium and in the course of its evolution. For this reason, a detailed knowledge of the mechanisms of force transmission should have a far-reaching impact on our understanding of dynamical phenomena, flow regimes, and plastic thresholds in granular materials. 5 In this paper, we study the probability distribution of contact forces both in two-dimensional ͑2D͒ and threedimensional ͑3D͒ systems. The robustness of the results with respect to material parameters, boundary conditions and numerical procedures will be emphasized. We show that in both 2D and 3D systems the distribution of interparticle normal forces is a power law with a negative exponent or nearly uniform for the forces below the mean force, and a decreasing exponential for larger forces, the values of the exponents depending on the average coordination number or the extent of departure from static equilibrium. Then, we introduce a model of contact forces in a packing at static equilibrium. In this model, the contact network is explicitly accounted for, so that no local force rule is needed like in other models of force transmission in granular packings. This model allows for an efficient computation of contact forces when a static solution exists. The force distributions obtained by this method are similar to those observed in contact dynamics and molecular dynamics simulations. Finally, we will discuss the relevance of the statistical properties of the force network to the overall description of force transmission in terms of the stress tensor.

II. SIMULATION METHODS
We used two different simulation methods: contact dynamics ͑CD͒ and molecular dynamics ͑MD͒. The CD method is based on the hard-particle approximation with three material parameters: the coefficient of friction and the coefficients of normal and tangential restitution. [6][7][8] The Coulomb friction law and the condition of particle impenetrability are prescribed without any regularization. In the MD method, a regularized form of the friction law ͑commonly a viscoplastic formulation of it͒ is used, and a repulsive interaction potential is introduced in order to take care of the hard core mutual exclusion of the particles. 9,10 Depending on the stiffness, the particles can be made more or less soft. Several schemes exist for the implementation of the method. For the MD simulations reported in this paper, we used a fifth-order predictor-corrector scheme. 11 The CD method is naturally more efficient than the MD method, and it provides a more precise prescription of the Coulomb friction law. Assuming that the particles are perfectly rigid, the CD method handles the nonsmooth motion ͑involving velocity jumps͒ of particles due to multiple collisions ͑transmission of momenta through the whole contact network͒ by means of the coefficients of restitution and a relaxation procedure similar to the Gauss-Seidel method. In the limit of dilute systems, the CD method becomes equivalent to the so-called event-driven method widely used for the simulation of granular gases. The phenomena in which elasticity is involved, such as wave propagation, cannot be simulated by the CD method. Such problems should rather be tackled by the MD method provided a fine time resolution is employed.
As far as the contact forces are concerned, we expect similar results from the two methods for stiff particles. The reason is that the contact forces in a noncohesive granular packing are basically determined by the contact network and the external forces acting on the particles. Differences may arise as a result of the organization of the contact network itself following the method used. For example, the coordination number, i.e., the average number of particles touching a particle, seems to decrease with the particle stiffness in MD simulations.
Another source of discrepancy is the ''nonsmooth'' feature of the dry friction law: The sliding velocity and the friction force at a contact belong to a set of admissible values which cannot be represented as a mathematical function. 6,7,12 This property makes a granular packing at static equilibrium hyperstatic, i.e., there are more forces to be determined than balance equations. 12 In the CD method, these indeterminacies are removed by the dynamics and force history. For example, in a pile at static equilibrium, the forces are partially determined by the preparation process. In the MD method, there is no indeterminacy because of the regularization of contact laws ͑contact elasticity and viscous friction around zero sliding velocities͒. Our simulations suggest that only very weak forces are influenced either by the history in CD simulations or by the force schemes in MD simulations. However, most of time, those weak forces are close to the limit of numerical precision. We may assume that they provide a physical limit of precision for our numerical results.

A. Two-dimensional systems
We simulated the same system of 4000 circular particles, once by the CD method and then by the MD method. Particle radii were uniformly distributed between 3.7 and 7.5 mm. The gravity and the particle-wall coefficient of friction were set to zero. The particle-particle coefficient of friction and the coefficient of normal restitution were 0.4 and 0.2, respectively. The particles were contained in a rectangular box. One of the walls was loaded by a constant force. The system was biaxially compressed by the inward motion of one of the walls. Then the compression was stopped and the system was allowed to relax to the static equilibrium. Figure 2 shows the probability density P of normal forces at static equilibrium for the CD and MD simulations. We see that the two distributions collapse on the same curve except when the force decreases toward zero, where we observe a small decrease of P in MD simulations compared to CD simulations. For forces N below the mean ͗N͘, referred to as ''weak forces,'' P is nearly uniform. For forces larger than the mean force, the ''strong forces,'' P decreases exponentially, The exponent ␤ equals roughly 1.4. The nearly uniform distribution of weak forces is an interesting feature of the distribution. The weak forces, comprising nearly 60% of contacts, are more frequent than the mean force itself.

B. Three-dimensional systems
We simulated two systems of 4000 particles contained in cylindrical box with two different coefficients of friction 0.1 and 0.4. The samples were prepared by a uniaxial compression under the action of a constant load applied on the upper lid. Figure 3 shows the distribution P of normal forces at static equilibrium. The two distributions collapse nearly on the same curve with the same value 1.4 of ␤ as in the above 2D systems. Figure 4 displays the force network in a thin layer parallel to cylinder axis. Strong ''force chains'' can easily be distinguished as in a 2D packing. The strongest chains have a linear aspect and they are mostly parallel to the axis of compression ͑vertical͒.
Our simulated system is similar to the experimental samples of Mueth et al., 13 who used a refined carbon paper technique for measuring the normal forces on the interior surfaces of a cylindrical vessel filled with spherical particles. They found a similar result for the distribution of normal forces with almost the same value 1.4 of ␤ as in our simulations, although they measured only wall-particle forces and the smallest forces they measured were only one decade below the mean force. We checked that the distribution P in our simulations is the same at the contacts along the confining walls as inside the packing.
Furthermore, Mueth and co-workers proposed the functional form which fits the distribution P over the whole range of the measured normal forces. Figure 5 shows that the form ͑2͒ with bϭ0.6 and ␤ϭ1.35 fits our data excellently, as well, except for N→0. Actually, a slight increase in P was observed in the experiments as N decreased toward zero. We have observed a similar increase in our simulations, but we checked that it was related to the numerical precision in the range of very weak forces, i.e., below 0.1% of the mean force. As argued by Mueth et al., 13 the above function for the range of weak forces provides a fit essentially indistinguishable from a power law N Ϫ␣ which we proposed following our early simulations, 14 as long as ␣ is positive and close to zero. However, as we shall see below, as soon as the coor- dination number decreases due, for example, to fast shearing, the distribution of weak forces is well described by a decreasing power law N Ϫ␣ and the functional form ͑2͒ is not suited any more.
Although the fitting form ͑2͒ suggests a continuous transition from weak forces to strong forces, it does capture the very different statistical behaviors of weak and strong forces. This difference is further materialized in the way the weak forces take part in the overall shear stress, geometrical anisotropy and frictional dissipation in the contact network. 15 Some aspects of such a ''bimodal'' character of the force network will be discussed in Sec. VI.

IV. DYNAMIC FORCES
The contact network plays the same role for the randomization of forces in a dense packing of particles as the binary collisions in a granular gas for the particle momenta. 13,16 In static equilibrium, the mean coordination number z is controlled by the confining pressure p, the coefficient of friction , and the dynamical processes from which the equilibrium state is reached. Nevertheless, whatever the values of p and , the requirement of static equilibrium for each particle together with geometrical disorder imply a mean coordination number larger than 3 in a 2D packing. In the CD simulations with circular particles, we find a value of z ranging from 3.2 to 3.8.
In contrast, when a granular packing is forced to deform, z may decrease below 3. This means that a subset of particles has either two contacts or no contacts at all. The forces are then less efficiently redistributed than at static equilibrium. This implies long-range force correlations within the contact network. Figure 6 shows an example of normal forces in a system of 3000 particles biaxially sheared by the fast motion of one of the walls. We observe a strong concentration of forces on contact lines ranging over the whole system and a ramified structure around them. These lines are naturally unstable and dynamical in nature. In fact, the whole force net-work is constantly modified in the course of deformation. What is the probability distribution of these ''dynamic'' forces?
We simulated three different systems, which we will refer to as A, B, and C. All of them are biaxially sheared with the same shear speed, but the confining pressure p decreases from system A to system C. Hence, the rate of energy injection and thus the dynamic motion of particles increases from A to C, and at the same time z decreases. The coordination numbers, considering only force-transmitting contacts, are 2.4, 2.1, and 1.6 in systems A, B, and C, respectively. We found that in all cases the forces larger than the average ͗N͘ are exponentially distributed, but the exponent ␤ is reduced from 1.4 in A to 1 in C. On the other hand, our data, represented on the log-log scale in Fig. 7, show that the weak forces have a power law distribution with an exponent ␣ increasing in absolute value from A to C. This distribution can be roughly written in the following form: where k is the normalization factor given by Moreover, considering the mean force ͗N͘ as the point of crossover between the two parts of the distribution, we get the following relation between the exponents: Note that the nearly uniform distribution of static forces is recovered by setting ␣ϭ0 in Eq. ͑3͒. Then, from Eq. ͑5͒ we get ␤ϭͱ2Ӎ1.4, which is the value we found for the distribution of static forces. Of course, it is desirable to find a single form fitting the whole range of forces for the distribution of dynamic forces. However, the writing in Eq. ͑3͒ reasonably describes the basic trends of the distributions observed in Fig. 7.
The power law divergence of the number of forces as the force decreases toward zero remains to be understood, as well as the nearly uniform distribution of weak static forces. Figure 6 suggests a self-similar branching process for weak forces. The absence of a force scale in the weak network means that the corresponding contacts are partially ''screened'' by archlike structures inducing local force correlations.

V. COMPUTING STATIC FORCES
A simple way to model the contact forces is to consider the components of the forces acting on each particle along a fixed direction in space and require that these components are exactly balanced on each particle. The granular medium can then be thought of as composed of parallel layers perpendicular to that direction, the particles of each layer touching those of adjacent layers. The idea behind the so-called ''q model'' 17 was to model the inherent disorder of the contact network by assuming that the incoming forces from a layer iϪ1 into the layer i are randomly redistributed into the layer iϩ1. This particular local force rule is the simplest one, but it can obviously be replaced by other rules in view of a better matching of the results to the data. The merit of this model is that it is amenable to an analytical treatment. Asymptotically, i.e., for i→ϱ, the model predicts a purely exponential force distribution P(F)ϰe ϪF/F 0 . This is an interesting result in itself compared to numerical and experimental observations although the model is based on crude approximations.
Of course, a layer-by-layer ''propagation'' of static forces through a local force rule is not realistic since it replaces a boundary value problem by a diffusive initial value problem in space. Moreover, the invariance by rotation requires a vectorial formulation of the problem. Even then, these conditions may be quite insufficient to provide the right distribution of forces. The reason is clear from the numerical results: the strong contacts form chains and loops, i.e., local correlations and nonlocal interactions, which may not easily be captured in a nondiffusive model. In other words, in order to tackle the problem of force distribution it is essential to account ab initio for the granular texture, i.e., the contact network, instead of introducing an arbitrary local force rule.
Indeed, we show here that the forces can be computed for a granular packing without any other physical input than the requirement of static equilibrium for each particle, given the contact network and external forces. In order to simplify the presentation, we consider only the normal forces, though the method can be generalized with no difficulty to include friction forces with the Coulomb friction law.
Consider two particles in the bulk of a granular medium touching at the contact c. From the equations of static equilibrium for each particle, the following relation is easily derived: where N is the force at the contact c, N in is the sum of the components of all forces pushing the two particles against one another along the contact normal ͑incoming forces͒, and N out is the sum of the components of all forces tending to separate the two particles from one another along the contact normal ͑outgoing force͒. This is a ''transfer equation'' for the contact c in the sense that it relates the force N acting between the two particles to in and out forces exerted by other particles. We thus transform the system of equations of static equilibrium for the particles into a set of transfer equations for the contacts. Physically, Eq. ͑6͒ means that a contact force in a granular medium is given by a difference. As a result, however large the values of in and out forces, the force at a contact can be arbitrarily weak. This is quite plausible in view of the nearly uniform probability distribution of weak forces. For a given contact, the in force is a compressive force, whereas the out force plays the same role as a tensile force. The limit of N in ϭN out implies Nϭ0, which corresponds to a full screening of the contact by the surrounding contacts.
The transfer equation suggests a very simple relaxation method for the computation of forces. Start with a set of arbitrary values N of all contact forces. From these values and the known external forces, compute the values of N in and N out and hence NЈϭ0.5(N in ϪN out ) for each contact. Then, do the following: and start again with NϭNЉ at each contact. For a given precision ␦Ͼ͉NЉϪNЈ͉, the iteration of this process converges to a solution which fulfills the transfer equation for each contact and hence the equations of static equilibrium for each particle. We found that the relaxation is smooth and fast. If the iteration does not converge, there is no static solution. Different aspects of this method, such as the influence of the underlying contact network and the initial values of forces, will be presented elsewhere. We note that other methods have been proposed as alternatives to the q model for computing static forces. [18][19][20] For the computation of dynamic forces, the transfer equation should be extended to include particle accelerations. Then, a similar method may be employed to calculate both forces and accelerations. The CD method itself may be viewed as an implementation of the above idea.
In Fig. 8, we have shown an example of the computed forces in a system of 15 000 disks of the same size forming a triangular contact network. In order to randomize the network, we removed a subset of contacts by requiring that the force is zero at those contacts. This ''substitutional disorder'' is strong enough to allow for an inhomogeneous distribution of forces. The only boundary force is applied on one of the walls for which a transfer equation is written in the same spirit as for the particles. Furthermore, in order to ensure a full equivalence between the set of equations of static equilibrium and the set of transfer equations, we consider the zero forces at nontransmitting contacts as external forces. We get three transfer equations per particle, since on a triangular network each particle has six contacts and each contact is shared by two particles. In this way, the two equations of static equilibrium for each particle transform into three transfer equations on three contacts with different directions. When a contact is considered as nontransmitting, this does not reduce the number of equations since the zero force at the contact is treated as an external force, but the number of unknown forces to be determined is thus reduced.
The computation time depends naturally on the precision required. For a precision of the order 0.001 times the mean force, the system shown in Fig. 8 can be computed in not more than half an hour on a workstation of moderate power. We remark that the friction forces can easily be included in this method and that the computation of the forces does not necessarily require a regular underlying contact network. In contrast to the q model, we found a uniform distribution for the weak forces and a decreasing exponential distribution for the strong forces.

VI. FORCES AND THE AVERAGE PRESSURE
What do we learn about the average pressure p in a granular medium from the force distributions? The pressure in a granular medium is a statistical quantity. From the expression of the stress tensor, the following relation between the average pressure and the mean force is easily derived: 16 where n c is the number of contacts per unit volume, and ͗l͘ is the average distance between the centers of touching particles. To establish the above expression, the weak correlation between the contact force N and the intercenter distance l has been neglected.
Equation ͑8͒ is similar to the expression of kinetic pressure in a molecular gas, namely pϭn p k B T, where n p is the number of particles per unit volume, k B is the Boltzmann constant, and T is the temperature. The distribution P of contact forces plays the same role with respect to the static pressure in a granular system as the Maxwell-Boltzmann distribution P B of particle velocities in a gas with respect to the kinetic pressure. The important difference between the two distributions is that, on one hand P is much wider than P B , and on the other hand, in contrast with P B , P shows no peak. The absence of a central tendency in the statistics of contact forces means that the pressure p does not by itself provide a sufficient average information about the full range of contact forces.
In this respect, given that the mean force is approximately the turning point from weak forces to strong forces in the distribution P, we may at least distinguish the contribu-tions p w and p s of weak and strong forces, respectively, to the overall pressure pϭp w ϩ p s . From the expression of P, we get In static equilibrium, for ␣ϭ0 and ␤ϭ1.4, we have p w Ӎ0.29p. This distinction between ''weak'' and ''strong'' pressures would remain purely formal if the weak forces contributed exactly the same way as strong forces to the mechanical properties of a granular medium. Our simulations by the CD and MD methods show that this is not the case. Figure 9 shows the contribution of the forces below to the stress ratio q/p as a function of in a simple shear deformation simulated by the MD method. A similar behavior was observed by the CD method. 15 We see that the weak forces give almost no contribution to the shear stress. The whole shear stress is sustained by the strong network. In other words, the weak pressure is mechanically similar to a fluid pressure. This means that, quite apart from its importance as a physical phenomenon in its own right, the local inhomogeneity of contact forces is not indifferent to the macroscopic description of a granular material. From a statistical point of view, this behavior implies that a granular medium is composed of two distinct and complementary phases. The ''weak phase'' behaves like a fluid, whereas the ''strong phase'' constitutes the solid backbone of the medium.

VII. CONCLUSION
We presented a set of numerical observations about the distribution of contact forces in a granular packing. Both contact dynamics and molecular dynamics methods were used for the simulations. We showed that the forces below the average force ͑weak forces͒ have a generic power law distribution, whereas the forces larger than the average force ͑strong forces͒ have a decreasing exponential distribution. The power law distribution of weak forces becomes a nearly uniform distribution at static equilibrium, in agreement with FIG. 9. Contribution of the forces below to the stress ratio q/ p as is increased from zero to the maximal force in a simple shear simulated by the MD method. the experiments of Mueth et al. 13 These robust features of force distribution are the same in both two-and threedimensional packings.
We introduced a simple model that allows for an efficient computation of static forces. The equations of static equilibrium are first transformed into a condition of force transfer across the contact between two touching particles. Then, a relaxation scheme is used over the transfer equations in order to compute the contact forces. We applied this method to a regular packing with a prescribed contact disorder, and found the same force distribution as observed in simulations. This model implies no local force rule assumed in other models of force transmission in granular packings.
Finally, we discussed the relevance of these findings to the mean pressure and shear stress in a granular packing. We argued that the contributions of weak and strong forces should be distinguished. A careful analysis of the data shows, indeed, that the weak forces do not contribute to the shear stress. The weak pressure is analogous to a fluid pressure, whereas the strong forces bear the whole shear stress of the medium.
Coming to the perspectives of this work, it is obvious that several aspects of the problem, such as the distribution and activation of friction forces and the behavior of dynamic forces in a gradual crossover from static equilibrium to a fully fluidized regime, are to be pursued. We think that the two-phase behavior of the contact forces is not only a property of the forces, but more fundamentally a property of the contact network itself. The analysis of the anisotropy of the contact network in connection with the forces shows indeed interesting correlations. 15 An investigation of the local cor-relations of the contact network may provide a key for understanding the origins of the force distributions as observed in numerical simulations and experiments.