<figure><img src="image_1.png" /></figure>

HAL
open science

# Multicontact dynamics of granular systems

Farhang Radjai

► To cite this version:

Farhang Radjai. Multicontact dynamics of granular systems. Computer Physics Communications, 1999, 121-122, pp.294-298. 10.1016/S0010-4655(99)00337-9. hal-01425239

HAL Id: hal-01425239

https://hal.science/hal-01425239

Submitted on 3 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multicontact dynamics of granular systems

Farhang Radjai

LMGC, UMR CNRS, Université Montpellier II, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France

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.

## 1. 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.

## 2. 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(\theta)$ with direction $\theta$ is usually well fitted by a truncated Fourier expansion [4]:

$$E(\theta) = \frac{c}{\pi} \{1 + A \cos 2(\theta - \theta_c)\}, \quad (1)$$

where $c$ is the total number of contacts, the parameter $A$ represents the amplitude of anisotropy, and $\theta_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 $\theta_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

<figure><img src="image_3.png" /><figcaption>Fig. 1. Polar diagram of the number of contacts $E(\theta)$ as a function of their direction $\theta$ in the system of Fig. 2. The solid line is a fit according to Eq. (1).</figcaption></figure>

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.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Velocity field during simple shearing of a system of 4070 particles.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Normal forces in the system of Fig. 2. Line thickness is proportional to the normal force.</figcaption></figure>
### 3. 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) *Unilaterality 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 $\mu$ 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 $[-\mu N, \mu 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–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

<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

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.

### 4. 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 $\langle N \rangle$, and a decreasing exponential law for larger forces [6,11]:

$$P_N(\xi) = \begin{cases} k\xi^{-\alpha} & \xi < 1, \\ ke^{\beta(1-\xi)} & \xi > 1, \end{cases} \quad (2)$$

where $\xi = N/\langle N \rangle$ 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 $\alpha'$ and $\beta'$. One might write the distribution law of forces as an exponentially decreasing function with a power-law prefactor, as

<figure><img src="image_7.png" /><figcaption>Fig. 5. Probability density $P_N$ of normal forces $\xi$ normalized with respect to the mean force.</figcaption></figure>
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 $\langle N \rangle = \int_0^\infty NP_N(N) dN$ is the point of
separation between the two parts of the distribution,
the exponents $\alpha$ and $\beta$ are related by $\beta^2 = (1 -
\alpha)(2 - \alpha)$. These exponents change only with the level
$\eta$ of kinetic agitation and the average coordination
number. $\alpha$ increases from zero with $\eta$. In a quasi-
static flow, the distribution of weak forces is nearly
uniform, i.e. $\alpha \simeq 0$ and $\beta \simeq 1.4$ (for tangential forces
$\alpha' \simeq 0.3$ and $\beta' \simeq 1$). Interestingly, our numerical
simulations show the same distribution (2) also in 3D
packings with exponents $\alpha \simeq 0$ and $\beta \simeq 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.

**5. 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

<figure><img src="image_8.png" /><figcaption>Fig. 6. Amplitude of anisotropy *A* of the ξ-network as a function of the force cutoff ξ in the system of Fig. 2; see text for the definition of the ξ-network.</figcaption></figure>

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* trans-
mission 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.

## 6. 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 force-texture 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.

## References

[1] Excellent reviews of the state of art can be found, for example, in: Physics of Granular Media, H.J. Herrmann, J.-P. Hovi,

S. Luding (Eds.) (Kluwer Academic Publishers, Dordrecht, 1998).

[2] I. Goldhirsch, G. Zanetti, Phys. Rev. Lett. 70 (1993) 1619.

[3] S. McNamara, W.R. Young, Phys. Rev. E 53 (1996) 5089.

[4] M. Oda, Soils Foundations 12 (2) (1972) 2.

[5] P. Dantu, in: Proc. 4th Int. Conf. Soil Mechanics and Foundations Engineering (Butterworths Scientific Publications, London, 1957).

[6] F. Radjai, M. Jean, J.J. Moreau, S. Roux, Phys. Rev. Lett. 77 (1996) 274.

[7] J.J. Moreau, Eur. J. Mech. A/Solids 13 (4) (Suppl.) (1994) 93.

[8] M. Jean, in: Mechanics of Geometrical Interfaces, A.P.S. Selvadurai, M.J. Boulon (Eds.) (Elsevier Science B.V., Amsterdam, 1995).

[9] F. Radjai, in: Physics of Granular Media, H.J. Herrmann, J.-P. Hovi, S. Luding (Eds.) (Kluwer Academic Publishers, Dordrecht, 1998).

[10] J. Schäfer, S. Dippel, D. Wolf, J. Phys. I (France) 6 (1996) 5.

[11] F. Radjai, D. Wolf, S. Roux, M. Jean, J.J. Moreau, in: Powders and Grains 97, R.P. Behringer, J.T. Jenkins (Eds.) (Balkema, Rotterdam, 1997) p. 211.

[12] D.M. Mueth, H.M. Jaeger, S.R. Nagel, Force distribution in a granular medium, preprint.

[13] F. Radjai, D. Wolf, M. Jean, J.J. Moreau, Phys. Rev. Lett. 80 (1998) 61.