Multicontact dynamics of granular systems

We present numerical results concerning force distributions and force-texture correlations in granular media. We show that two modes of stress transmission coexist in a granular assembly, and discuss the robustness of this behavior with respect to material parameters, boundary conditions, and simulation methods


Introduction
Granular materials are composed of macroscopic particles that interact only through contacts. Although non cohesive granular materials exhibit fluid-like and solid-like states, they behave in a very different way from molecular solids and fluids. Recent active research has uncovered interesting and often counterintuitive effects in the dynamics of granular media, and has raised deep questions as to the relevance of familiar concepts to this field [1].
The aim of this short contribution is to discuss the origins and some important features of multicontact states of granular systems, i.e. those states that involve a network of contacts through which the momenta propagate more efficiently than by the motion of individual particles.

Multicontact states
Fundamentally, three microscopic features play a key role in the dynamics of granular media: mutual exclusion of particles, dissipative nature of interactions due to plastic deformation during collisions or sliding friction, and dynamically-induced randomness as the only source of thermodynamic disorder that sets by itself the energy scale. As a consequence of these microscopic features, arguments based on entropy may not work any more. An open dissipative system can in principle self-organize and give rise to structures. In fact, dense clusters appear even in strongly excited granular systems [2]. This process may lead to the collapse of particles into multicontact sub-structures, where the interparticle contacts form a connected network through which the momenta can propagate over distances far larger than the particle size [3]. It is obvious that the evolution of such structures can not be described in terms of a sequence of binary collisions and thus thermodynamic description and any event-driven algorithm for their simulation breaks down. Each collision in the multicontact state is a multiple collision due to the propagation of momenta in the contact network.
The simplest example of a fully multicontact system is a granular packing in quasi-static flow or simply in static equilibrium. In traditional fields, such as soil mechanics, a longstanding effort has been spent in order to rationalize the mechanical response of granular materials to shear for specified boundary con-ditions. Elasto-plastic models describe in many respects the strain-stress behavior of granular materials in quasi-static flow. This macroscopic approach deals with granular medium as a continuum whose mechanical behavior is described in terms of a stress tensor, a strain-rate tensor and some internal variables. In the 60s, however, two important discoveries revived interest in the microscopic aspects of granular materials: (1) The contact network is generically anisotropic.
This means that the contact normals are not randomly directed in space. In a properly sheared packing, the density of contacts E(θ) with direction θ is usually well fitted by a truncated Fourier expansion [4]: where c is the total number of contacts, the parameter A represents the amplitude of anisotropy, and θ c is the direction for which the maximum of E is reached. Starting a system with an initially isotropic texture (A = 0), A increases with the shear strain and θ c coincides with the major principal direction of the strain-rate tensor. Fig. 1 shows a polar diagram of E for a system of 4070 particles subject to simple shear. The velocity field is displayed in Fig. 2. (2) In contrast to the highly uniform density of a close-packed assembly, the distribution of contact forces is highly heterogeneous [5]. The forces can be as large as six times the average force in static equilibrium and nearly 60% of contacts carry a force below the average force [6]; see Fig. 3. These observations hint at a physics very different from that of dilute systems and, they also suggest that the mechanical state of the system may be governed by several internal variables associated with the texture and the modes of force transmission inside the system.

Numerical methods
For the investigation of multicontact states, we performed extensive numerical simulations of systems containing as large as 4000 circular particles. For these simulations, we used two different methods: contact dynamics (CD) and molecular dynamics (MD). The CD method prescribes two contact laws [7,8]: (1) Unilaterallity condition: If the separation velocity v n of two particles is positive, then the normal force N is zero. If, on the other hand, v n is zero, i.e. the two particles stay in contact, then N can have a positive indefinitely large value so as to prevent interpenetration. This is shown as a graph, Signorini's graph, in Fig. 4(a). (2) Coulomb's friction law: If the sliding velocity v t is nonzero, then the friction force T resists sliding and its value is given by the coefficient of friction μ times the normal force N . If, on the other hand, v t is zero (non-sliding contact), then T can take any value in the interval [−μN, μN]. This is shown in Fig. 3(b) as a graph. Both Signorini's condition and Coulomb's law are nonsmooth in the sense that the two conjugate variables {v n , N} or {v t , T } belong to a continuous set of acceptable values which can not be represented as a mathematical function. It can be shown that dynamics removes locally this degeneracy and allows to calculate all forces and velocities with no resort to any regularization trick [7][8][9]. In the MD method, "regularized" forms of the above laws are implemented [10]. In the simplest MD algorithms, the vertical branches of the above graphs (v n = 0 and v t = 0) are replaced by steep lines.
In this way, the main physical input in the two methods is essentially the same. Both methods have Fig. 4. (a) Signorini's graph relating normal force N and separation velocity v n at a contact between two particles; (b) Coulomb's graph relating friction force T and sliding velocity v t at a contact. been applied with much success for the investigation of granular phenomena. We note, however, that the CD method does not imply a fine time resolution and is thus faster than the MD method. Moreover, it allows for a more precise handling of friction than in the MD method. On the other hand, the CD method does not prescribe the elasticity of contacts (particles are infinitely stiff) and hence cannot be used for the investigation of the phenomena, such as sound propagation, that involve the deformation of particles. The most significant results of our simulations can be summarized around two main topics: statistical distribution of forces and force-texture correlations.

Force distributions
Our simulations by MD and CD methods show that the statistical distribution P N of normal forces is generically a power law, with a negative exponent, for the forces below the mean normal force N , and a decreasing exponential law for larger forces [6,11]: where ξ = N/ N is the force normalized by the mean force and k is the normalization factor. One example is shown in Fig. 5. The distribution P T (T ) of tangential forces has the same form as Eq. (2) with different exponents α and β . One might write the distribution law of forces as an exponentially decreasing function with a power-law prefactor, as well. This alternative writing, however, does not fit the transition part between the two behaviors better than Eq. (2). The writing of the force distribution in the form of Eq. (2), which completely separates the two parts of the distribution, is motivated by the fact that, as we shall see below, the contacts that carry a force below the mean force, to which we will refer as "weak contacts" in the following, do not play the same role with respect to the overall resistance of the system to shear as the contacts carrying a force larger than the mean force ("strong contacts"). As the average force N = ∞ 0 NP N (N) dN is the point of separation between the two parts of the distribution, the exponents α and β are related by β 2 = (1 − α)(2 − α). These exponents change only with the level η of kinetic agitation and the average coordination number. α increases from zero with η. In a quasistatic flow, the distribution of weak forces is nearly uniform, i.e. α 0 and β 1.4 (for tangential forces α 0.3 and β 1). Interestingly, our numerical simulations show the same distribution (2) also in 3D packings with exponents α 0 and β 1.4 at static equilibrium. These exponents change very slightly with the coefficient of friction. The size-dispersion of particles does not seem to influence the force distributions either. This distribution with almost the same values of exponents has been recently confirmed by careful experiments by means of the carbon paper technique for forces of a granular packing at the contacts with the walls of its container [12].
The distribution P N of normal forces in Eq. (2) plays the same role with respect to the static pressure of a granular system as the Maxwell-Boltzmann distribution of particle velocities in a gas with respect to the kinetic pressure. The two distributions have the same exponential tail, but the important difference between the two distributions is that in P N the weak forces are at least as frequent as the average force.

Force-texture correlations
In the above statistical analysis of contact forces, we considered all forces independently of directions of the contacts through which they are transmitted. There is, however, an interesting correlation between forces and contact directions, which determine the anisotropy of the texture. A simple way to study this correlation is to consider the subset of contacts which carry a force below a given cutoff ξ . We shall refer to this subset as the "ξ -network". The variation of a quantity such as the anisotropy evaluated for the "ξ -network" as ξ is varied from 0 to the maximal force in the system, allows then to estimate its correlation with the contact force. Fig. 6 shows the amplitude of anisotropy A as a function of ξ in a simply sheared system at a given moment of deformation. Surprisingly, the direction of anisotropy is orthogonal to the major principal direction of the strain-rate tensor (negative values) as long as ξ is below the average force [13]. In other words, the direction of the anisotropy of the weak network is orthogonal to the axis of compression, whereas that of the strong network is parallel.
It is also possible to evaluate the contribution of the ξ -network to the overall stress inside the system. We found that the shear stress Q for ξ < 1 is negligibly small compared to the total deviation load sustained by the system. Those forces only contribute 28% of the average pressure in the medium. This means that the weak network behaves essentially like an interstitial liquid. The whole shear stress is sustained by the strong network. We also found that the whole dissipation due to sliding friction takes place inside the weak network, while all contacts in the strong network are non-sliding.
These properties, which are well reproduced both by CD and MD simulations, indicate a bimodal transmission of stress in multicontact systems. The texture is composed of two phases, a weak phase and a strong phase, that are intimately correlated. This means that for the description of stress transmission in a granular medium, two stress tensors are needed.

Conclusion
We discussed features of the multicontact states of granular systems at variance with the dilute states. After a brief account of the major physical inputs in numerical simulations of granular systems, we presented two sets of numerical observations concerning the statistical distribution of contact forces and the forcetexture correlations. These observations suggest that the texture of a granular material is composed of two complementary phases and raise the problem of a theoretical understanding of the origins of this bimodal behavior.