The contact dynamics (CD) method, also called nonsmooth contact dynamics (NSCD), is a discrete element method (DEM) for the simulation of granular materials. It emerged from a mathematical formulation of nonsmooth dynamics and the subsequent algorithmic developments by J. J. Moreau and M. Jean (Moreau [1977, 1983, 1988a,b, 1993, 1994], Jean and Pratt [1985], Jean [1988], Jean and Moreau [1992], Jean et al. [1994], Jean [1995, 1999]). The fundamental difference between this method and the common DEM or molecular dynamics (MD) approach lies in the treatment of small length and time scales involved in the dynamics of granular media. In MD-type DEM, pioneered by P. Cundall, the particles are treated as rigid bodies but the contacts between particles are assumed to be compliant and obeying a viscoelastic behavior in which the local strain variables are the relative particle positions and displacements (Cundall [1971], Cundall and Strack [1979], Thornton and Yin [1991], Herrmann [1993], Thornton [1993], Pöschel and Buchholtz [1995], Thornton [1997], Luding [1998], Matuttis et al. [1999], McNamara and Herrmann [2004], Garca and Medina [2007], Gilabert et al. [2007], Richefeu et al. [2007]). The time-stepping schemes used for the numerical integration of the equations of motion imply thus a fine resolution of the small time and length scales involved in contact interactions. In the CD method, these small scales are neglected and their effects absorbed into contact laws together with a nonsmooth formulation of particle dynamics described at a larger scale than small elastic response times and displacements.

The CD method has been applied to investigate granular materials (Moreau [1997], Radjai et al. [1996b, 1998], Bratberg et al. [2002], Radjai and Roux [2002], Staron et al. [2002], Nouguier-Lehon et al. [2003], Renouf et al. [2004], McNamara and Herrmann [2004], Taboada et al. [2005], Saussine et al. [2006], Azéma et al. [2007], Ries et al. [2007]), as well as other mechanical systems composed of rigid bodies with frictional contact interactions such as masonry and tensegrity structures (Acary and Jean [1998], Nineb et al. [2006]). The results prove often to be in good agreement with experimental observation, and for static and plastic shear properties with MD simulations (Radjai et al. [1995], Moreau [1997], Radjai et al. [1997a, 1999], Lanier and Jean [2000], Radjai and Roux [2004]). The differences between the two methods arise mainly from the scales of description and the presence of elastic parameters in the MD method (Radjai et al. [1997a], McNamara and Herrmann [2004, 2006]).

We present here the CD method as a consistent model of nonsmooth and multicontact granular dynamics expressed in contact coordinates. Nonsmoothness refers to various degrees of discontinuity in local or global characteristics of a dynamical system. The mathematical concepts and tools for the treatment of nonsmooth dynamics were developed in relation with mechanical problems involving unilateral constraints and in the context of convex analysis; see Brogliato [1999] for a detailed history. The multicontact feature is present in static states and in dense flows of granular materials where spatial correlations occur at large length scales and impulse dynamics is mixed with smooth particle motions at different time scales (Bershadskii [1994], Lvoll et al. [1999], Puglisi et al. [2002], Radjai and Roux [2002], Silbert et al. [2002], Staron et al. [2002], Pouliquen [2004], Majmudar and Behringer [2005], Agnolin and Roux [2007], Behringer et al. [2007], Olsson and Teitel [2007]).

While explicit integration schemes used in the molecular dynamics approach are rather straightforward to implement in a numerical code, the underlying model of the CD method seems to be more complex and thus less accessible to computer implementation. Here, we would like to highlight the intuitive aspects of the CD method and to illustrate its simple but powerful logics in treating multiple contacts.

The dynamics of a homogeneous granular flow involves at least three different
scales: 1) small time and length scales characterizing contact interactions, 2)
intermediate scales associated with particle rearrangements and shear rate, and 3)
large length scales related to geometrical correlations at even larger scales. The
elastic response time τ_{e} of the particles and their contacts is far below the
rearrangement time characterized by the mean time interval τ_{c} between successive
events (collisions, …). In the same way, the relative stiffness E∕p, where E is the
Young modulus of the particles and p the mean confining stress, is normally high,
and hence the elastic displacements λ_{e} are tiny compared to the mean particle size
d.

Various DEM algorithms differ in the treatment of small scales. The common
denominator is that the particles are considered as rigid bodies and thus the
material behavior is attributed to the contact zones, with local strain variables
derived from rigid-body degrees of freedom (translations and rotations) of the
particles. When the small scales are numerically resolved, as in MD-like or
smooth DEM, the particle motions are smooth (twice differentiable) and the
equations of dynamics are integrated with the help of force laws governing
particle interactions. The latter expresses the contact force (,) between a
pair of touching particles as a (mono-valued) function of their contact
displacement (a displacement vector with respect to a reference state or
defined geometrically) and its time derivative . With an explicit integration
scheme, the time step δt must be small compared to τ_{e} for the sake of
numerical stability. Numerical efficiency often imposes the use of simple
force laws (linear elastic with viscous damping and a Coulomb friction
law).

An interesting issue is whether one can or should neglect the small sub-particle scales and turn to a model of granular dynamics at the scale of particle rearrangements. In other words, when only plastic granular flow is of interest, the contact elastic strains may become irrelevant. Numerically, this also means that large time steps can then be used and materials with very stiff particles simulated. Another motivation is that, when deformable particles are considered and modeled, e.g. within a finite-element approach, it is more plausible to consider their mutual contacts as rigid, i.e. as purely unilateral constraints. This is conform to the requirement of non-interpenetration between solid bodies.

Such a particle-scale model of granular dynamics implies, however, velocity jumps.
Let δt be the time resolution (a time step in numerical simulations). Since δt ≫ τ_{e},
the net effect of fast collisional dynamics in a multicontact system over the time
interval δt is a finite change of particle velocities and impulsive forces through the
contact network as well as a net dissipation due to contact plasticity or
viscosity. This is similar to shock dynamics in a dilute gas governed by
binary collisions except that in the latter the impulse is confined to a single
contact at a time whereas in a multicontact system the whole system is
affected. This nonsmooth dynamics can not be described by second order
equations of motion since the accelerations are no more defined. The model
should thus necessarily be formulated at the velocity level with built-in
discontinuities.

The CD method is based on a particle-scale model of granular dynamics with two major ingredients. First, the interactions between particles are described by contact laws instead of force laws. In this framework, a contact is a dynamic and unilateral element with a coarse-grained behavior over the considered time interval δt. Secondly, the equations of dynamics are adapted to account consistently for both smooth and nonsmooth motion in a time-stepping scheme. Below, we discuss in detail these features of the CD method.

Let us consider two particles i and j and a potential contact point α between them. We assume that a unique contact plane (line in 2D) tangent to the two particles at α can be geometrically defined so that the contact can be endowed with a local reference frame defined by a normal unit vector and a tangential unit vector lying on the tangent plane. The orientation of the axes is a matter of convenience. In the following, the subscripts n and t refer to normal and tangential components, respectively, with regard to the contact frame.

A potential contact point between two particles has the following dynamic
content. As long as the distance δ_{n} between the two particles remains positive
(corresponding to a gap), no force is activated and the normal force f_{n} is
identically zero. But when δ_{n} = 0, a nonnegative (repulsive) normal force f_{n} is
mobilized at the contact point and it can take indefinitely large values depending
on the forces acting on the two particles; see Fig. 1. These conditions define
a complementary relation, called Signorini’s conditions, between δ_{n} and
f_{n}:

(1) |

Only one of the two alternatives is true. This mutual disjunction between the two statements is fundamental as it implies that none of the two variables can be reduced to a mono-valued function of the other. In other words, Signorini’s conditions define a degenerate, i.e. multi-valued, interaction law. It can be represented as a graph displayed in Fig. 2. Symbolically, we represent this complementarity relation as follows:

δ_{n}f_{n}.
| (2) |

A contact is persistent if both δ_{n} = 0 and u_{n} = _{n} = 0. The normal force
vanishes at a breaking contact, i.e. when δ_{n} = 0 and u_{n} > 0. Hence, Signorini’s
complementarity relation can be developed as follows:

(3) |

In this form, Signorini’s conditions contain a kinematic constraint: For a contact,
i.e. for δ_{n} = 0, we have

u_{n}f_{n}.
| (4) |

The full dynamic content of Signorini’s conditions appears for the second derivative
_{n} = _{n} with respect to time where the relative acceleration _{n} and f_{n} will obey
Signorini’s complementarity relation. However, in a nonsmooth formulation where
velocity jumps are expected, we have to stay at the velocity level and the relation
(3) should further be reinterpreted according to the coarse-graining time, as we
shall see below.

The three alternatives in (3) lead to the condition of a complete contact law formulated by Moreau as follows (Moreau [1994]):

(5) |

In this writing, δ_{n} ≤ 0 corresponds to the condition of a geometrical contact
which was expressed in (3) and (1) by an equality. The first and third
conditions are the same as in (3). The second condition is necessary to
ensure that the particle motions will respect the unilateral nature of a
contact. Up to numerical precision, this condition in CD simulations prevents
the violation of unilateral constraints in the absence of a repulsive force
law.

In the context of nonsmooth motion, the time derivative _{n}(t) is not unique.
Assuming bounded variation, we thus distinguish between the left-limit velocity
u_{n}^{-} and the right-limit velocity u_{
n}^{+} which in a time-stepping formulation should
be considered as the contact velocities at times t and t + δt, respectively.
With a finite time resolution δt, the actual velocity is immaterial since
only left and right velocities (and the velocity jumps) enter the dynamics;
see section 4. In analogy with a binary collision, the main unknown of
the problem is u_{n}^{+} that we would like to calculate given the approach
velocity u_{n}^{-} attributed to the beginning of the time step in a time-stepping
scheme. The contact model depends on the choice of the velocity u_{n} and
the nature of the force f_{n} involved in the complementarity relation (3) or
the relations (5) of a complete contact law in the context of nonsmooth
evolution.

In a granular flow, the normal velocity u_{n} at a contact evolves during a time
interval δt either as a result of smooth motion or due to impulses generated
by collisions. These impulses propagate through the contact network so
that a contact may experience several successive impulses during δt. Such
events can be resolved for a sufficiently small time increment δt or they
may be tracked according to an event-driven scheme. The event-tracking
strategy is, however, numerically inefficient, of limited applicability and in
contradiction with the scope of the CD method based on coarse-grained
dynamics.

Moreau recognized that the approximation of the contact force f_{n} during δt is a
measure problem in the mathematical sense (Moreau [1994, 2004]). A static or
regular force f^{s} is the density of the measure f^{s} dt with respect to dt. In contrast,
an impulse p generated by a collision has no density with respect to dt. In other
words, the forces at the origin of the impulse are not resolved at the scale δt. In
practice, however, we can not differentiate these contributions in a coarse-grained
dynamics, and the two contributions should be summed up to a single measure
df′_{n}. Then, the contact force is defined as an average of this measure over
δt:

f_{n} = ∫
_{t}^{t+δt}df′_{
n}
| (6) |

With this definition, the equations of dynamics can be formulated as measure
differential equations, as we shall see in section 4. Except in static equilibrium
where the forces are of static origin, f_{n} generally depends on time resolution. For
large time steps, the fast dynamics, e.g. successive collisions or rearrangement
events, is partially filtered out and the remaining effect appears as an average force
p∕δt, in addition to the static part f_{n}^{s}. For small enough time steps, the fast
dynamics prevails and the corresponding contribution p∕δt increases.

As to the velocity u_{n} involved in Signorini’s condition (3), a simple and physically
motivated choice is to assume that u_{n} is a weighted mean between u_{n}^{-} and u_{
n}^{+}:
u_{n} = η u_{n}^{-} + (1 - η)u_{
n}^{+}, where η is a material parameter characterizing the
contact. The value η = 0 seems plausible as it corresponds to a complementarity
Signorini relation between the mean force f_{n} and the right-limit velocity u_{n}^{+} which
is the main unknown of the problem. However, for u_{n} = 0 this choice leads to
u_{n}^{+} = 0. But, according to (3) nonzero forces occur only for u_{
n} = 0. Hence,
this choice implies that, independently of the dynamics of the system,
nonzero forces will occur only at persistent (u_{n}^{+} = 0) contacts. This choice
eludes thus the treatment of collisions with nonzero normal restitution
coefficient e_{n} for which the contact is not persistent (u_{n}^{+} > 0) although the
corresponding force or impulsion is nonzero. For η≠0, a binary shock implies
-u_{n}^{+}∕u_{
n}^{-} = η∕(1 - η). Identifying this ratio with e_{
n}, one gets η = e_{n}∕(1 + e_{n}).
Hence, we set

u_{n} = .
| (7) |

Notice that, in contrast to the common definition of restitution coefficient for
binary collisions, Eq. (7) is not a kinematic relation between u_{n}^{-} and u_{
n}^{+}; It is
simply the expression of a weighted mean velocity that is assumed to obey
Signorini’s conditions u_{n}f_{n}. Its exact role in the CD model will become clear
in combination with the equations of dynamics.

The Coulomb law of dry friction is, by definition, a complementarity relation
between the friction force f_{t} and the sliding velocity u_{t} at a contact point between
two particles:

(8) |

where μ is the coefficient of friction and it is assumed that the unit tangential
vector t points in the direction of the sliding velocity so that _{t} ⋅ = u_{t}. The graph
of this complementarity relation is displayed in Fig. 3. Like Signorini’s conditions,
this is a degenerate law that can not be reduced to a mono-valued function between
u_{t} and f_{t}. For a concise expression of Coulomb’s law, we will use the following
notation:

u_{t}f_{t}.
| (9) |

As in the case of Signorini’s conditions, it is useless to consider higher order
developments (for u_{t} = 0 and _{t} = 0) in the context of nonsmooth motion. This
can be useful only in a system or for system parameters for which the
contact network does not evolve and the contact status (sliding/nonsliding,
force-transmitting/nontransmitting) fully characterizes the mechanical state. In the
MD method, the second condition (the vertical branch) of Coulomb’s friction law
(8) is replaced either by a tangential elastic relation or by a viscous law. This
regularization transforms the complementarity relation into a mono-valued function
at the price of introducing small local strain or strain-rate parameters into the
problem.

In the framework of the CD model, Coulomb’s law needs to be reformulated to
incorporate velocity discontinuities and impulsions for a large integration time δt.
The points discussed with regard to Signorini’s conditions in section 3.1 apply also
to Coulomb’s conditions. The force f_{t} at a contact represents the average effect of
static and impulsive forces experienced by the contact during the time lag δt. The
sliding velocity u_{t} entering Coulomb’s law (8) is thus a mean velocity. Along the
same lines as for normal contact reactions, a simple model consistent with
tangential restitution is to set

u_{t} = .
| (10) |

where e_{t} ∈ [-1, 1] is the tangential restitution coefficient. This expression of u_{t} is
assumed to obey Coulomb’s law u_{t}f_{t}.

The complementarity relations (3) and (8) together with the expressions (7) and
(10) of the mean velocity define a coarse-grained model of frictional contact with
three material parameters μ, e_{n} and e_{t}. This model can be used to treat frictional
contact problems for deformable or undeformable bodies. It has the unique feature
of combining static and impulsive forces in the same framework. The rich
content of this model can be appreciated only in interplay with the equations
of dynamics, as we shall see below for granular media composed of rigid
particles.

The rigid-body motion of the particles in a granular assembly is governed by
Newton’s equations under the action of imposed external bulk or boundary forces
_{ext}, and the contact reaction forces ^{α} exerted by neighboring particles at the
contact points α. The contact reactions may involve contact torques but torque
transmission will not be considered here. We also consider mainly a 2D system
since our focus here is on the method rather than its applications to various
geometries. An absolute reference frame with unit vectors (, ŷ) is assumed, and we
set ẑ = ×ŷ. Each particle is characterized by its mass m, moment of inertia I,
mass center coordinates , mass center velocity , angular coordinates θ, and
angular velocity ωẑ. For a smooth motion (twice differentiable), the equations of
motion of a particle are

(11) |

where = ∑
_{α} ^{α} and = ẑ ⋅∑_{
α} ^{α} × ^{α} where ^{α} is the contact vector joining
the center of mass to the contact α and _{ext} represents the moment of external
forces.

For a nonsmooth motion with time resolution δt involving impulses and velocity discontinuities, an integrated form of the equations of dynamics should be used. As discussed in section 3, every force sums up static and impulsive actions so that it should naturally be represented as a measure d′ rather than a density with respect to time. It has a density with respect to dt only for a smooth motion. Hence, the equations of dynamics should be written as an equality of measures:

(12) |

where d′ = ∑
_{α}d ′^{α} and d′ = ẑ ⋅∑_{
α} ^{α} × d ′^{α}. These measure differential
equations can be integrated over δt with the definition of as an approximation of
the integral over d′:

∫
_{t}^{t+δt}d′ = δt
| (13) |

In the same way, we set ∫
_{t}^{t+δt}d′ = δt. With these definitions, the integration
of Eq. (12) over δt yields

(14) |

where (^{-},ω^{-}) and (^{+},ω^{+}) are the left-limit and right-limit velocities of the
particle, respectively. In this form, the accelerations are replaced by velocity
jumps ^{+} -^{-} and ω^{+} - ω^{-}, and Newton’s equations take the form of an
equality between the change of momenta and the average impulse during
δt.

Although Eqs. (14) may be considered as resulting from a direct integration of smooth equations (11) over a short time lag, the formulation of dynamics in terms of measures is a necessary step in a rigorous mathematical formulation undertaken by Moreau and other precursors of nonsmooth mechanics. For example, from this formulation, it becomes clear that, according to (13), is a coarse-grained force. This means that, its dynamic content depends on the time resolution δt. Clearly, in a quasi-static granular flow, the contact forces reflect basically the confining stresses which are regular static forces applied on the system, so that the value of the time step in CD simulations will have negligibly small effect on the calculated forces. In contrast, for a faster flow and in dynamic transients during quasi-static flow, the value of time step matters. This should not, however, be considered as a numerical artifact but rather as a physical effect reflecting the nonsmooth character of granular flow. Such effects of the integration time are also observed in MD simulations. For example, it was shown that the velocity distributions and their spatial correlations obtained from the integration of particle displacements crucially depend on the integration time (Radjai and Roux [2002]).

The equations of dynamics can be written in a compact form for a set of N_{p}
particles by using matrix representation. The particles are labelled with
integers i ∈ [1,N_{p}]. The forces and force moments F_{x}^{i},F_{
y}^{i},^{i} acting on
the particles i are arranged in a single high-dimensional column vector
represented by a boldface letter F belonging to ℝ^{3Np}. In the same way,
the external bulk forces F_{ext,x},F_{ext,y},_{ext} applied on the particles and
the particle velocity components U_{x}^{i},U_{
y}^{i},ω^{i} are represented by column
vectors F_{ext} and U, respectively. The particle masses and moments of
inertia define a diagonal 3N_{p} × 3N_{p} matrix denoted by M. With these
notations, the equations of dynamics (14) are cast into a single matrix
equation:

M(U^{+} -U^{-}) = δt(F + F_{
ext})
| (15) |

In the MD method, the equations of motion of the particles are integrated with the help of force laws which express the contact forces as a function of the actual mechanical state (particle positions and velocities) of the system. In the CD method, the contact forces are not explicit functions of the state. Hence, the forces and velocities should be determined at the same time. To this end, as the contact laws are expressed in contact variables, we need to express the equations of dynamics in the same variables.

The contacts are labelled with integers α ∈ [1,N_{c}], where N_{c} is the total
number of contacts. Like particle velocities, the contact velocities u_{n}^{α} and
u_{t}^{α} can be collected in a column vector u ∈ ℝ^{2Nc}. In the same way, the
contact forces f_{n}^{α} and f_{
t}^{α} are represented by a vector f ∈ ℝ^{2Nc}. We
would like to transform the equations of dynamics from F and U to f and
u.

The formal transformation of matrix equations (15) is straightforward. Since the contact velocities u are linear in particle velocities U, the transformation of the velocities is an affine application:

u = G U | (16) |

where G is a 2N_{c} × 3N_{p} matrix containing basically information about the
geometry of the contact network. A similar linear application relates f to
F:

F = H f | (17) |

where H is a 3N_{p} × 2N_{c} matrix. We will refer to H as contact matrix. It contains
the same information as G in a dual or symmetric manner. It can easily be shown
that

H = G^{T }
| (18) |

where G^{T } is the transpose of G. This property can be inferred from the equivalence
between the power F ⋅U developed by ’generalized’ forces F and the power f ⋅u
developed by the bond forces f. In general, the matrix H is singular and, by
definition, its null space has a dimension at least equal to 2N_{c} - 3N_{p}. A schema of
this transformation from particle dynamics to transfer equations is displayed in Fig.
4.

The matrix H^{iα} can be decomposed into two matrices H_{
n}^{iα} and H_{
t}^{iα} such
that

(19) |

and

F^{i} = ∑_{
α}(H_{n}^{iα} f_{
n}^{α} + H_{
t}^{iα} f_{
t}^{α})
| (20) |

Using these relations, Eqs. (15) can be transformed into two equations for each contact α:

(21) |

We now can make appear explicitly linear relations between the contact variables from Eqs. (21) and definitions (7) and (10). We set

_{k1k2}^{αβ} = ∑_{
i,j}H_{k1}^{T,αi}M^{-1,ij}H_{
k2}^{jβ},
| (22) |

where k_{1} and k_{2} stand for n or t. With this notation, Eqs. (21) can be rewritten as

The coefficients _{k1k2}^{αα} for each contact α can be calculated as a function of the
contact geometry and inertia parameters of the two partners 1_{α} and 2_{α} of
the contact α. Let _{i}^{α} be the contact vector joining the center of mass of
particle i to the contact α. The following expressions are easily established:

An alternative representation of Eqs. (23) and (24) is

The two offsets a

In order to solve the system of 2N_{c} transfer equations (in 2D) with the
corresponding complementarity relations, we proceed by an iterative method which
converges to the solution simultaneously for all contact forces and velocities. We
first consider a single-contact situation to which we will refer as the local SC
problem (SC standing for Signorini-Coulomb). It consists of the determination of
contact variables f_{n}^{α}, f_{
t}^{α}, u_{
n}^{α} and u_{
t}^{α} at a single contact given the values of the
offsets a_{n}^{α} and a_{
t}^{α} at the same contact. Formally, by combining the transfer
equations (26) and (27) with the complementarity relations (4) and (9), it is easily
shown that the local SC problem is equivalent to the following relations:

The solution of this problem is given by intersecting the lines representing transfer
equations with Signorini’s and Coulomb’s graphs; see Fig. 6. The intersection
occurs at a unique point due to the positivity of the coefficients _{k1k2}^{αα} (positive
slope). In other words, the dynamics removes the degeneracy of the contact
laws. Notice, however, that the two intersections can not be established
separately when _{nt}^{αα}≠0. To find the local solution, one may consider
the intersection of transfer equations with the force axis, i.e. by setting
u_{n} = u_{t} = 0. This yields two values g_{n}^{α} and g_{
t}^{α} of f_{
n}^{α} and f_{
t}^{α}, respectively:

In order to illustrate how this works, let us consider a simple contact problem
between two disks of masses m_{1} and m_{2} and radii R_{1} and R_{2} lying on the x-axis
and subjected to constant external forces F^{1} and -F^{2} , respectively, where
F^{1} ≥ F^{2} ≥ 0; see Fig. 7. We also assume that the two particles come in touch with
a relative velocity u_{n}^{-} = U_{
x}^{2} -U_{
x}^{1} ≤ 0. The normal coefficient of restitution is e_{
n}
and the inverse reduced mass is _{nn} = (m_{1} + m_{2})∕(m_{1}m_{2}) since for a
disk c_{t} = ⋅ = 0; see equation (25). From equations (28) and (29), we
have

a_{n} = -(1 + e_{n}) u_{n}^{-} +
| (36) |

The transfer equation is

_{nn}f_{n} = (1 + e_{n})u_{n} + a_{n}
| (37) |

The intersection of this line with Signorini’s graph occurs on the vertical branch if
a_{n} > 0, and this is the case. Hence, u_{n} = 0, and we have u_{n}^{+} = -e_{
n}u_{n}^{-}
and

f_{n} = a_{n} = -(1 + e_{n}) u_{n}^{-} +
| (38) |

This is an interesting expression as it shows the two origins of the total normal
force experienced by the contact. The first term is an impulsive force induced by
the impact velocity u_{n}^{-} and averaged over the time lag δt whereas the second term
is a static force induced by the applied external forces F^{1} and F^{2}. The first
term vanishes only if the impact velocity is zero, i.e. if the two particles
initially touch. In this case, we have u_{n}^{+} = 0 so that the two particles stay in
contact. The particle velocities can be calculated from the equations of
motion of each particle and the value of f_{n}: U^{1+} = U^{1-} + δt(F^{1} - f_{
n})∕m_{1}
and U^{2+} = U^{2-} + δt(-F^{2} + f_{
n})∕m_{2}. In the case U^{1-} = U^{2-} = 0, we get
U^{1+} = U^{2+} = (F^{1} -F^{2})∕(m_{
1} + m_{2}) which corresponds to the velocity of the center
of mass of the two particles.

In a multicontact system, the terms b_{n}^{α} and b_{
t}^{α} in the offsets a_{
n}^{α} and a_{
t}^{α} depend
on the forces and velocities at contacts β≠α; see equations (28), (28), (30) and (31).
Hence, the solution for each contact depends on all other contacts of the system
and it must be determined simultaneously for all contacts. We will refer to this
problem as the global SC problem.

An intuitive and robust method to solve the global SC problem is to search the
solution as the limit of a sequence {f_{n}^{α}(k),f_{
t}^{α}(k),u_{
n}^{α}(k),u_{
t}^{α}(k)} with α ∈ [1,N_{
c}].
Let us assume that the transient set of contact forces {f_{n}^{α}(k),f_{
t}^{α}(k)} at the
iteration step k is given. From this set, the offsets {a_{n}^{α}(k),a_{
t}^{α}(k)} for all contacts
can be calculated through the relations (28) and (29). The local SC problem can
then be solved for each contact α for these values of the offsets, yielding an
updated set of contact forces {f_{n}^{α}(k + 1),f_{
t}^{α}(k + 1)}. This correction
procedure is equivalent to the solution of the following local SC problem:

Remark that this force update procedure does not require the contact
velocities u_{n}^{α}(k + 1),u_{
t}^{α}(k + 1)} to be calculated as the offsets depend
only on the contact forces. The set {f_{n}^{α}(k),f_{
t}^{α}(k)} evolves with k by
successive corrections and it converges to a solution satisfying the transfer
equations and complementarity relations at all potential contacts of the system.
The iteration can be stopped when the set {f_{n}^{α}(k),f_{
t}^{α}(k)} is stable with
regard to the force update procedure within a prescribed precision criterion
ε_{f}:

< ε_{f} ∀α.
| (41) |

Finally, from the converged contact forces, the particle velocities {^{i}} can be
computed by means of the equations of dynamics (14). The efficiency of this
iterative scheme depends on how the information (imposed boundary and bulk
forces or velocities) propagates through the system. Clearly, the information
propagates faster if the updated contact force (f_{n}^{α}(k + 1),f_{
t}^{α}(k + 1)) at a contact
α is used to update also the values of the offsets (a_{n}^{β}(k),a_{
t}^{β}(k)) at the
neighboring contacts β that will be treated next. This sweeping of the contact
network might be optimized to some extent depending on the expected
solution.

The iterative procedure depicted above provides a robust method which proves
efficient in the context of granular dynamics. Note that the information is treated
locally and no large matrices are manipulated during iterations. Fig. 8 shows
several snapshots of the normal force network in a packing of disks in the course
of iterations. The left and bottom walls are immobile while the top and
right walls are free to move and subjected to the same confining stress
p. The sample is dense and its time evolution is not of interest here; we
rather calculate the forces and velocities by solving iteratively the global SC
problem for this system. In the time-stepping scheme this is equivalent to the
resolution of the problem over a single time step. The initial values of
the contact forces are set to zero, i.e. f_{n}^{α}(0) = 0 and f_{
t}^{α}(0) = 0 for all
α. We see that the information (boundary forces) propagates from the
walls smoothly throughout the system with increasing number of iterations
(130, 365, 1000 and 2000). If the iterations were stopped before the full
propagation of the information, the solution would involve strong force
gradients.

The evolution of the probability distribution function of normal forces during
iterations is shown in Fig. 9. When normalized by the average normal force ⟨f_{n}⟩ at
the converged state, the distribution appears to increase in width and it becomes
stable beyond 1000 iterations but carries the signature of the well-known
exponential shape of strong contact forces just after a few iterations. The number
N_{i} of necessary iterations to converge is strongly dependent on the precision ε_{f}
but not on δt. Fig. 10 shows N_{i} as a function of ε_{f} for a given global SC
problem in a packing of 2500 particles with initialization of the forces to zero.
We see that N_{i} diverges quite fast as ε_{f} is decreased. It should, however,
be noted that the number of iterations is substantially reduced when the
iteration is initialized with a globally correct guess of the forces. This is
the case in an evolution problem where the forces at each time step can
be initialized with the forces computed at the preceding step; see section
6.

The number of iterations is partially controlled by the special features of force
distribution in granular media. As a result of arching, the probability density
function P(f_{n}) of normal forces shows no or weak central tendency. A small peak
is observed when the system is isotropic. But in all cases the number of
forces tends either to a finite value P(0) as f_{n} → 0 or diverges (Radjai
et al. [1997b], Mueth et al. [1998], Radjai et al. [1999], Antony [2001], Silbert
et al. [2002], Kruyt [2003], Metzger [2004], Majmudar and Behringer [2005], Richefeu
et al. [2006], Azéma et al. [2007], Kruyt and Antony [2007]). Hence, a huge
number of contacts (nearly 60% in a weakly polydisperse packing) carry small
forces whereas the strong forces are less in number but exponentially distributed. It
is important to remark that the weak forces cannot simply be neglected as they
play an organic role in the overall stability of the assembly Radjai et al. [1998].
Due to the absence or weakness of central tendency, an increasing number of weak
contacts should be resolved as ε_{f} is decreased, requiring thus an increasingly larger
number of iterations. In sum, for a prescribed precision ε_{f}, the forces below ε_{f} ⟨f_{n}⟩
are not correctly calculated. This provides a physical criterion which may
be employed to calibrate the required precision in dealing with a specific
problem.

A high precision over forces might be required when the static equilibrium of a packing is studied. During a quasistatic flow, a high precision is necessary in the investigation of transition between successive states separated by small time intervals δt. The trajectory is all the more determinist as the precision on forces, and hence on velocities, is high. For longer evolutions, when some degree of ergodicity may be assumed, it might be physically justified to consider lower precision in each step at the level of the resolution of the global SC problem. The choice of the precision can be combined with a more or less large time-step δt.

It should also be briefly mentioned that the uniqueness of solution of the global SC
problem in a multicontact system is not proved. We have 3N_{p} equations of
dynamics and 2N_{c} complementarity relations. The unknowns of the problem are
3N_{p} particle velocities and 2N_{c} contact forces. Since there are equal numbers of
unknowns and dynamic or contact relations, and since the solution of the local SC
problem is unique, the global solution could have been unique, too. The point is
that the 2N_{c} contact complementarity relations are not equations but
inequations. Thus, the extent of indeterminacy of the solution reflects all
possible combinations of contact forces accommodating the complementarity
relations.

In fact, the CD method provides a suitable framework to study the indeterminacy
of granular assemblies by searching self-equilibrated solutions with external forces
and particle velocities set to zero in the expressions of a_{n}^{α} and a_{
t}^{α}; Eqs. (28) and
(29). For a given configuration of particles, the standard SC problem can thus be
solved by setting external forces and particle velocities to zero and for an
ensemble {f_{n}^{α}(0),f_{
t}^{α}(0)} varied according to a statistical distribution. The
ordering of the sequence of updated contacts during the iterative procedure
may also affect the solution. It happens that for a generically disordered
granular system, self-equilibrated solutions can not be found unless in small
partially or totally ordered systems. This implies that the CD treatment of a
granular assembly involves practically no indeterminacy. Thus, if slightly
different solutions are found according to the iteration method employed,
they should be attributed to the lack of precision below the convergence
criterion ε_{f}, affecting mainly the very weak contacts. The indeterminacy can
be studied also analytically by considering the null space of the contact
matrix H (Radjai et al. [1996a], McNamara et al. [2005]). The degree of
indeterminacy may be high, but it does not imply significant force variability
since the solutions are restrained by the complementarity relations. This is
also consistent with a common observation in CD simulations that the
force fluctuations, when they occur, are basically of impulsive nature and
not spurious effects that could be attributed to the presence of scale-free
self-stresses.

Finally, a few words about the role of the restitution coefficient e_{n} is useful. In a
binary collision, the contact opening is fully controlled by e_{n} as in the simple
example discussed in section 5.1. The situation is different in a multicontact system
where a contact can open even if u_{n}^{-} = 0. In other words, contact opening is
governed by the multicontact dynamics and not only by the local restitution. In
contrast, contact closing, i.e. the creation of a persistent contact, depends only on
the local restitution coefficient. This is because, according to equation
(7), with e_{n} > 0 the conditions u_{n}^{+} = 0 and u_{
n}^{-} < 0 imply u_{
n} < 0 which
contradicts Signorini’s conditions. Hence, in the CD method a collision
between two particles leads to the creation of a persistent contact only if
e_{n} = 0.

The problem of multiple shocks is a complex issue that should be treated at very
small time scales of elastic waves, and this is out of the scope of the contact
dynamics model which was designed to find mechanically admissible solutions of
multicontact dynamics at much larger time scales. Experimental observation,
however, suggests that in a multicontact system multiple shocks may occur and
dissipate the whole kinetic energy at very short times so that the effective
restitution coefficient should indeed be set to zero (Luding et al. [1994a,b], Dippel
et al. [1996], da Cruz et al. [2005]). This effect is the reminiscence of inelastic
collapse (McNamara and Young [1992]). When the time increment δt is large
enough compared to the characteristic time of successive collisions, the choice
e_{n} = 0 is relevant. In dense systems, this choice allows for the opening and closing
of contacts due to multicontact dynamics. In a loose assembly, the time scales are
different and the effective restitution coefficient might be nonzero for coarse-grained
dynamics.

It should also be mentioned that the restitution coefficients e_{n} and e_{t} in
binary collisions are not pure material parameters but depend on the surface
geometry and impact velocity (Lun and Savage [1986], Kuwabara and
Kono [1987], Smith and Liu [1992], Foerster et al. [1994], Labous
et al. [1997], Luding [1998], Brilliantov et al. [2004, 2007]). There is no
fundamental difficulty in taking such effects into account in a CD algorithm. This is
also true for the coefficient of friction μ which may depend on the sliding velocity.
It should, however, be remarked again that the manifestation of physical effects of
this kind is different under multicontact conditions with respect to binary
interactions.

The global SC problem may well occur as an event at particular instances of a granular flow. The iterative resolution method presented above can then be applied to calculate the contact forces and particle velocities at those instances. Hence, the global SC problem may be embedded in an event-driven algorithm which can be a useful approach in certain circumstances. The global SC problem may also be considered for the estimation of the stability of a static assembly of particles. In the CD method, the global SC problem is associated with a time-stepping scheme.

In order to set up this scheme, we need to come back to the contact laws and
remark that the first condition of the Signorini relation in (3) is the only condition
referring to space coordinates. Actually, the SC problem was formulated
at the velocity level for both dynamics and contact laws, and the first
Signorini condition was accounted for by assuming that only the potential
contacts, where δ_{n} ≤ 0, were involved in the SC problem. In other words,
the contact network is defined explicitly from particle positions in a SC
problem and it will no more evolve during the coarse-graining time interval
δt.

Another important feature of the global SC problem is the prospective
character of the treatment expressed by the second condition of a complete
contact law (5). This means that the right-limit velocities are calculated such
that the complementarity relations will not be violated by the subsequent
motion of the particles. This condition is achieved through an appropriate
definition of the mean velocities u_{n} and u_{t} in Eqs. (7) and (10) entering the
complementarity relations. In other words, the treatment is fully implicit and
the right-limit velocities {^{i+},ω^{i+}} may be used to increment particle
positions.

The above two remarks devise the following time-stepping scheme. Let t and t + δt
be the considered time interval. The configuration {^{i}(t)} and particle velocities
{^{i}(t),ω^{i}(t)} are given at time t. The latter play the role of left-limit velocities
{^{i-},ω^{i-}} in the global SC problem. The contact network {α,^{α}, ^{α}} is set up
from the configuration at time t or from an intermediate configuration {_{m}^{i}}
defined by

_{m}^{i} ≡^{i}(t) + ^{i}(t).
| (42) |

When this configuration is used for contact detection, other space-dependent
quantities such as the inverse mass parameters _{k1k2}^{αα} and external forces _{
ext}^{i}
should consistently be defined for the same configuration and at the same time
t + δt∕2.

Then, the global SC problem is solved iteratively over the contact network and the
right-limit particle velocities {^{i+},ω^{i+}} are calculated. The latter correspond to
the velocities at the end of the time step t + δt:

This scheme is unconditionally stable due to the implicit discretization.
Hence, no damping parameters at any level are needed. For this reason, the
time step δt can be large. The real limit imposed on the time step is the
occurrence of cumulative numerical errors leading to undesired excess overlaps
between the particles. Of course, such overlaps have no dynamic significance
in the CD method, but they falsify the geometry and thus the evolution
of the system. By construction, if there are initially no overlaps between
the particles, the contact dynamics ensures that no overlaps will occur in
the course of evolution. But the particle positions are secondary in the
CD method and they are updated from the integration of the velocities.
Even for a high precision ε_{f} , the calculated velocities and contact forces
by solving the SC problem involve some numerical imprecision that may
lead to excessive overlaps at a number of contacts in the long run. Hence,
to avoid such effects, one requires both a sufficiently high precision and
rather small time steps. In high-quality simulations of shear flow with
≃ 10^{4} particles, the typical value of the time step is ≃ 10^{-4} s. This value is
by several orders of magnitude larger than the usual values of the time
step in MD simulations. It should also be borne in mind that the time
step δt is not a precision criterion for the evolution problem in the CD
method. The precision is mainly controlled by ε_{f}. The time step should rather
be considered as a coarse-graining parameter for nonsmooth dynamics.
It should be reduced if the impulse dynamics at small time scales is of
interest.

An important feature of granular dynamics is the occurrence of highly nonlinear and subtle transitions at small time scales. These short-time phenomena include sharp impulsive transitions, collective rearrangements and frictional interlocking. As a result, the dynamics over successive time steps needs variable treatment. For a rather smooth evolution less effort should be consumed than when a major rearrangement event takes place. In the CD simulations, the number of iterations represents the required effort at each step. For the same level of precision, the number of iterations varies considerably in the course of time stepping. One example is shown in Fig. 11. It can be checked that, if the initializing forces at each step are retrieved from the last step, the number of iterations is correlated with dynamic events. Hence, in the CD method subtle rearrangement events are always calculated with the required high precision. In some cases, a maximum number of iterations may be imposed in order to increase efficiency, but this is equivalent to reducing precision when a larger number of iterations are necessary. Moreover, incomplete relaxation leads to wave-like perturbations in time evolution (Unger et al. [2002]).

The basics of the CD (contact dynamics) method for discrete element simulation of granular materials were presented. This method can be viewed as the algorithmic formulation of nonsmooth granular dynamics at the scale of particle rearrangements where small elastic response times and displacements are neglected. Two major issues arising in this context are treated in the framework of the CD method: 1) The contact laws expressed as complementarity relations between the contact forces and velocities and 2) The nonsmooth motion involving velocity jumps with impulsive unresolved forces as well as smooth motion with well-resolved static forces. A consistent description of the dynamics at the velocity level leads naturally to an implicit time-stepping scheme together with an explicit treatment of the evolution of the contact network. Most concepts developed in this formulation are rather intuitive although they rely on a sound mathematical background of nonsmooth dynamics and convex analysis.

The focus of this text was on the intuitive features of the CD method with respect
to subtle collective phenomena involved in the multicontact dynamics of granular
media. In particular, we discussed the role of the coarse-graining time δt, the
precision issue in the iterative procedure for the resolution of the global SC
(Signorini-Coulomb) problem with respect to force distributions, and the
relevance of the coefficients of restitution e_{n} and e_{t} and their interpretation in
a multicontact system. The CD method is characterized by its unique
feature of bringing together in the same formalism two limit regimes of
granular dynamics: 1) collisional regime governed by binary shocks and
incomplete energy restitution and 2) static regime governed by multiple
contacts, geometrical disorder, force balance and dynamic rearrangements.
Hence, this method provides a suitable framework for the investigation of
dense granular flows where smooth evolutions are intermingled with sharp
transitions.

The CD method can also be considered as an adequate framework for the numerical treatment of frictional contact problems. Indeed, the Coulomb friction and perfectly rigid contact condition are implemented in an exact form, i.e. without introducing artificial penalization parameters or damping. Given a contact network, all kinematic constraints implied by contact laws are simultaneously taken into account together with the equations of dynamics in order to determine the velocities and contact forces in the system. This global SC problem is solved by an iterative process pertaining to the Gauss-Seidel iterative method that consists of solving a single SC problem at each contact, and successively and iteratively updating the forces until a convergence criterion is fulfilled (Jeffreys and Jeffreys [1988]). The method is thus capable of dealing properly with the nonlocal character of the momentum transfers resulting from the impenetrability of the particles. It can be employed to study stiff systems for which smooth MD-like methods require small time steps for numerical stability and the stiffness matrix may become ill-conditionned as the contact network evolves.

The CD method is unconditionally stable due to its inherent implicit time integration scheme. The uniqueness of the solution at each time step is not guaranteed for perfectly rigid particles. However, the variability of admissible solutions is generally below numerical precision. The variability resulting from numerical precision can be reduced and the calculations significantly accelerated by initializing the iterative procedure at each step with the forces computed in the preceding step.

The basic algorithm presented in this paper can be (and has been) extended to deal with richer contact laws, various particle shapes and more efficient resolution of the global SC problem in 2D and 3D. Several variants of the time-stepping scheme can be found in (Jean [2001]) and a parallelization strategy in (Renouf et al. [2004]). The contact laws can be supplemented with a complementarity relation between a torque and a contact spin variable (Bratberg et al. [2002]). Using such a complementarity relation, the rolling friction is easily implemented in this framework. Adhesion forces can be introduced by a simple shift of the complementarity relations:

A particular attention should be paid to the origin of the contact forces in the CD method. An example is the uniaxial compression of a dense granular material by imposing a constant velocity on a wall. In the MD simulation of this problem, the displacement of the wall causes mainly the elastic deformation of the particles and the contact forces increase accordingly. In the CD simulation of the same problem, since the contact laws involve no force scale and no static boundary or bulk force are applied, the force scale is fixed by the imposed velocity through the impulsive terms in Eqs. (28) and (29). Since no rearrangements can occur due to a too high density, the corresponding kinetic energy is not dissipated. This energy increases adiabatically, the contact forces increase proportionally to the displacement and the particles interpenetrate. In contrast, if the uniaxial compression is controlled by an increasing boundary force, the contact forces increase in proportion to the applied force as in MD simulations, and, in the absence of particle rearrangements, the contact reaction forces balance exactly the driving force so that the packing stays indefinitely in static equilibrium. In the MD approach, the static forces are fully encoded in the particle positions with a scale given by the stiffness. In the CD approach, there is no such force scale and thus the static force scale should be defined externally.

Finally, a future work on the coarse-grained granular dynamics might help to get more insight into the dissipative and impulsive phenomena that are modeled by their average effects in the CD method. The suggestion is to apply the same methodology as in (Radjai and Roux [2002]) in order to evaluate the influence of time resolution. One can first perform a standard MD simulation of a shear flow, then take averages over the dynamics for a given integration time δt and finally compare this coarse-grained dynamics with CD simulations of the same flow with a time step δt. Repeating this procedure for several values of δt, will allow us to appreciate the CD method and compare it in a meaningful way with the MD method.

V Acary and M Jean. Numerical simulation of monuments by the contact dynamics method. In DGEMN-LNEC-JRC, editor, Monument-98, Workshop on Seismic Performance of Monuments, pages 12–14. LNEC, 1998.

Ivana Agnolin and Jean-Noël Roux. Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks. Phys Rev E Stat Nonlin Soft Matter Phys, 76(6-1):061302, Dec 2007.

S. J. Antony. Evolution of force distribution in three-dimensional granular media. Phys Rev E, 63(1 Pt 1):011302, Jan 2001.

E. Azéma, F. Radjai, and G. Saussine. Quasistatic rheology and force transmission in a packing of polyhedral particles. Submitted., 2008.

Emilien Azéma, Farhang Radjai, Robert Peyroux, and Gilles Saussine. Force transmission in a packing of pentagonal particles. Phys. Rev. E, 76(1 Pt 1):011301, Jul 2007.

R. P. Behringer, Karen E Daniels, Trushant S Majmudar, and Matthias Sperl. Fluctuations, correlations and transitions in granular materials: statistical mechanics for a non-conventional system. Philos Transact A Math Phys Eng Sci, Aug 2007.

A. Bershadskii. Stochastic density waves of granular flows: strong-intermittent dissipation fields with self-organization. Physica A, 210: 386–390, 1994.

I. Bratberg, F. Radjai, and A. Hansen. Dynamic rearrangements and packing regimes in randomly deposited two-dimensional granular beds. Phys. Rev. E, 66:031303–1, 2002.

Nikolai Brilliantov, Clara Saluea, Thomas Schwager, and Thorsten Pschel. Transient structures in a granular gas. Phys Rev Lett, 93(13):134301, Sep 2004.

Nikolai V Brilliantov, Nicole Albers, Frank Spahn, and Thorsten Pschel. Collision dynamics of granular particles with adhesion. Phys Rev E Stat Nonlin Soft Matter Phys, 76(5 Pt 1):051302, Nov 2007.

B Brogliato. Nonsmooth mechanics. Springer, London, 1999.

P. Cundall. A computer model for simulating progressive large scale movements of blocky rock systems. In Proceedings of the Symposium of the International Society of Rock Mechanics, Vol. 1, 132-150, 1971.

P. A. Cundall and O. D. L. Strack. A discrete numerical model for granular assemblies. Géotechnique, 29(1):47–65, 1979.

Frdric da Cruz, Sacha Emam, Michal Prochnow, Jean-Nol Roux, and Franois Chevoir. Rheophysics of dense granular materials: discrete simulation of plane shear flows. Phys Rev E Stat Nonlin Soft Matter Phys, 72(2 Pt 1):021309, Aug 2005.

S. Dippel, G. G. Batrouni, and D. E. Wolf. Collision-induced friction in the motion of a single particle on a bumpy inclined line. Phys.Rev. E, 54: 6845, 1996.

S. F. Foerster, M. Y. Louge, H. Chang, and K. Allia. Measurements of the collision properties of small spheres. Phys. Fluids, 6(3):1108–1115, 1994.

Xavier Garca and Ernesto Medina. Acoustic response of cemented granular sedimentary rocks: molecular dynamics modeling. Phys Rev E Stat Nonlin Soft Matter Phys, 75(6 Pt 1):061308, Jun 2007.

F. A. Gilabert, J-N. Roux, and A. Castellanos. Computer simulation of model cohesive powders: influence of assembling procedure and contact laws on low consolidation states. Phys Rev E Stat Nonlin Soft Matter Phys, 75 (1 Pt 1):011303, Jan 2007.

H. J. Herrmann. Molecular dynamics simulations of dry granular media. In The first Nisshin Engineering Particle Technology International Seminar: Discrete Particle Simulations in Powder Technology, page 8, Osaka, Japan, 1993.

M. Jean. Unilateral contact and dry friction: time and space variables discretization. Arch. of Mech., Warszawa, 40(1):677–691, 1988.

M. Jean. Frictional contact in rigid or deformable bodies: numerical simulation of geomaterials, pages 463–486. Elsevier Science Publisher, Amsterdam, 1995.

M. Jean. The non smooth contact dynamics method. In Special issue on modeling contact and friction. 1999.

M Jean. Micromécanique des matériaux granulaires, chapter Simulation numérique discre, pages 199–324. Hermes, Paris, 2001.

M. Jean and J. J. Moreau. Unilaterality and dry friction in the dynamics of rigid body collections. In Proceedings of Contact Mechanics International Symposium, pages 31–48, Lausanne, Switzerland, 1992. Presses Polytechniques et Universitaires Romandes.

M. Jean and E. Pratt. A system of rigid bodies with dry friction. International Journal Eng. Sci., pages 497–513, 1985.

M. Jean, F. Jourdan, and B. Tathi. Numerical dynamics for the simulation of deep drawing. In Proc. of IDDRG 1994, Handbooks on Theory and Engineering Applications of Computational Methods, May 16-17 1994.

H. Jeffreys and B. S. Jeffreys. Methods of Mathematical Physics. Cambridge, England, 1988.

N. P. Kruyt. Contact forces in anisotropic frictional granular materials. International Journal of Solids and Structures, 40(13-14):3537–3556, 2003.

N. P. Kruyt and S. J. Antony. Force, relative-displacement, and work networks in granular materials subjected to quasistatic deformation. Phys. Rev. E, 75(5):051308–8, May 2007.

G. Kuwabara and K. Kono. Restitution coefficient in a collision between two spheres. Japanese Journal of Applied Physics, 26(8):1230–1233, 1987.

L. Labous, A. D. Rosato, and R. Dave. Measurements of collision properties of spheres using high-speed video analysis. Phys. Rev. E, 56: 5715, 1997.

J. Lanier and M. Jean. Experiments and numerical simulations with 2d disks assembly. Powder Technology, 109:206–221, 2000.

S. Luding. Collisions and contacts between two particles. In H. J. Herrmann, J.-P. Hovi, and S. Luding, editors, Physics of dry granular media - NATO ASI Series E350, page 285, Dordrecht, 1998. Kluwer Academic Publishers.

S. Luding, E. Clément, A. Blumen, J. Rajchenbach, and J. Duran. Anomalous energy dissipation in molecular dynamics simulations of grains: The “detachment effect”. Phys. Rev. E, 50:4113, 1994a.

S. Luding, E. Clément, A. Blumen, J. Rajchenbach, and J. Duran. Studies of columns of beads under external vibrations. Phys. Rev. E, 49(2): 1634, 1994b.

C. K. K. Lun and S. B. Savage. The effects of an impact velocity dependent coefficient of restitution on stresses developed by sheared granular materials. Acta Mechanica, 63:15–44, 1986.

G. Lvoll, K. J. Mly, and E. G. Flekky. Force measurements on static granular materials. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics, 60(5 Pt B):5872–5878, Nov 1999.

T. S. Majmudar and R. P. Behringer. Contact force measurements and stress-induced anisotropy in granular materials. Nature, 435(7045): 1079–1082, Jun 2005.

H.-G. Matuttis, S. Luding, and H. J. Herrmann. Discrete element methods for the simulation of dense packings and heaps made of spherical and non-spherical particles. Powder Technology, 1999. in press.

S. McNamara and W. R. Young. Inelastic collapse and clumping in a one-dimensional granular medium. Phys. Fluids A, 4(3):496, 1992.

Sean McNamara and Hans Herrmann. Measurement of indeterminacy in packings of perfectly rigid disks. Phys. Rev. E, 70(6 Pt 1):061303, Dec 2004.

Sean McNamara, Ramn Garca-Rojo, and Hans Herrmann. Indeterminacy and the onset of motion in a simple granular packing. Phys. Rev. E, 72(2 Pt 1):021304, Aug 2005.

Sean C McNamara and Hans J Herrmann. Quasirigidity: some uniqueness issues. Phys. Rev. E, 74(6 Pt 1):061303, Dec 2006.

Philip T Metzger. Granular contact force density of states and entropy in a modified edwards ensemble. Phys Rev E Stat Nonlin Soft Matter Phys, 70(5 Pt 1):051303, Nov 2004.

J. J. Moreau. Evolution problem associated with a moving convex set in a hilbert space. Journal of differential equations, 26:347–374, 1977.

J. J. Moreau. Liasons unilatérales sans frottement et chocs inélastiques. Comptes Rendus de l’Académie des Sciences, 296:1473–1476, 1983.

J. J. Moreau. Bounded variation in time. In P.D. Panagiotopoulos and G Strang, editors, Topics in Nonsmooth Mechanics, pages 1–74. Bikhäuser, Basel, 1988a.

J. J. Moreau. New computation methods in granular dynamics. In Powders & Grains 93, page 227, Rotterdam, 1993. A. A. Balkema.

J. J. Moreau. Numerical investigation of shear zones in granular materials. In D. E. Wolf and P. Grassberger, editors, Friction, Arching, Contact Dynamics, pages 233–247, Singapore, 1997. World Scientific.

J.J. Moreau. Unilateral contact and dry friction in finite freedom dynamics, volume 302 of International Centre for Mechanical Sciences, Courses and Lectures. Springer, Vienna, 1988b.

J.J. Moreau. Some numerical methods in multibody dynamics: Application to granular materials. European Journal of Mechanics A/Solids, supp.(4):93–114, 1994. Formulation mathematiques tire du livre Contacts mechanics.

J.J. Moreau. An introduction to unilateral dynamics. In M. Frémond and F. Maceri, editors, Novel approaches in civil engineering, number 14 in Lecture Notes in Applied and Computational Mechanics, pages 1–46. Springer-Verlag, 2004.

D. M. Mueth, H. M. Jaeger, and S. R. Nagel. Force distribution in a granular medium. Phys. Rev. E, 57:3164, 1998.

S. Nineb, P. Alart, and D. Dureisseix. Approche multi-échelle des systèmes de tenségrité. Revue Européenne de Mécanique Numérique, 15: 319–328, 2006.

C. Nouguier-Lehon, B. Cambou, and E. Vincens. Influence of particle shape and angularity on the behavior of granular materials: a numerical analysis. Int. J. Numer. Anal. Meth. Geomech, 27:1207–1226, 2003.

Peter Olsson and S. Teitel. Critical scaling of shear viscosity at the jamming transition. Phys Rev Lett, 99(17):178001, Oct 2007.

T. Pöschel and V. Buchholtz. Molecular dynamics of arbitrarily shaped granular particles. J. Phys. I France, 5(11):1431–1455, 1995.

O. Pouliquen. Velocity correlations in dense granular flows. Phys. Rev. Lett., 93(24):248001–4, December 2004.

A. Puglisi, A. Baldassarri, and V. Loreto. Fluctuation-dissipation relations in driven granular gases. Phys Rev E Stat Nonlin Soft Matter Phys, 66(6 Pt 1):061305, Dec 2002.

Radjai, Evesque, Bideau, and Roux. Stick-slip dynamics of a one-dimensional array of particles. Phys. Rev. E, 52(5):5555–5564, Nov 1995.

F. Radjai and S. Roux. Contact dynamics study of 2d granular media : Critical states and relevant internal variables. In H. Hinrichsen and D. E. Wolf, editors, The Physics of Granular Media, pages 165–186, Weinheim, 2004. Wiley-VCH.

F. Radjai, L. Brendel, and S. Roux. Nonsmoothness, indeterminacy, and friction in two- dimensional arrays of rigid particle. Phys. Rev. E, 54(1): 861, 1996a.

F. Radjai, J. Schäfer, S. Dippel, and D. Wolf. Collective friction of an array of particles: A crucial test for numerical algorithms. J. Phys. I France, 7:1053, 1997a.

F. Radjai, D. E. Wolf, S. Roux, M. Jean, and J. J. Moreau. Force networks in dense granular media. In R. P. Behringer and J. T. Jenkins, editors, Powders & Grains 97, pages 211–214. Balkema, Rotterdam, 1997b.

F. Radjai, D. E. Wolf, M. Jean, and J.J. Moreau. Bimodal character of stress transmission in granular packings. Phys. Rev. Letter, 80:61–64, 1998.

Farhang Radjai and Stphane Roux. Turbulentlike fluctuations in quasistatic flow of granular media. Phys Rev Lett, 89(6):064302, Aug 2002.

Farhang Radjai, Michel Jean, Jean-Jacques Moreau, and Stphane Roux. Force distributions in dense two-dimensional granular systems. Phys. Rev. Lett., 77(2):274–, July 1996b.

Farhang Radjai, Stephane Roux, and Jean Jacques Moreau. Contact forces in a granular packing. Chaos, 9(3):544–550, Sep 1999.

Mathieu Renouf, Frederic Dubois, and Pierre Alart. A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media. Journal of Computational and Applied Mathematics, 168 (1-2):375–382, July 2004.

V. Richefeu, F. Radjai, and M. S. El Youssoufi. Stress transmission in wet granular materials. Eur. Phys. J. E, Feb 2007.

Vincent Richefeu, Moulay Saïd El Youssoufi, and Farhang Radjai. Shear strength properties of wet granular materials. Phys. Rev. E, 73(5 Pt 1): 051304, May 2006.

Alexander Ries, Dietrich E Wolf, and Tams Unger. Shear zones in granular media: three-dimensional contact dynamics simulation. Phys Rev E Stat Nonlin Soft Matter Phys, 76(5 Pt 1):051301, Nov 2007.

G. Saussine, C. Cholet, P.E. Gautier, F. Dubois, C. Bohatier, and J.J. Moreau. Modelling ballast behaviour under dynamic loading. part 1: A 2d polygonal discrete element method approach. Computer Methods in Applied Mechanics and Engineering, 195(19-22):2841–2859, April 2006.

Leonardo E Silbert, Gary S Grest, and James W Landry. Statistics of the contact network in frictional and frictionless granular packings. Phys Rev E Stat Nonlin Soft Matter Phys, 66(6 Pt 1):061303, Dec 2002.

C. E. Smith and P.-P. Liu. Coefficients of restitution. J. Appl. Mech., 59:963, 1992.

L. Staron, J.-P. Vilotte, and F. Radjai. Preavalanche instabilities in a granular pile. Phys. Rev. Lett., 89:204302, 2002.

A. Taboada, K. J. Chang, F. Radjai, and F. Bouchette. Rheology, force transmission, and shear instabilities in frictional granular media from biaxial numerical test using the contact dynamics method. Journal Of Geophysical Research, 110:1–24, 2005.

C. Thornton. Computer simulation of impact fracture/fragmentation. In The first Nisshin Engineering Particle Technology International Seminar: Discrete Particle Simulations in Powder Technology, page 17, Osaka, Japan, 1993.

C. Thornton. Coefficient of restitution for collinear collisions of elastic-perfectly plastic spheres. Journal of Applied Mechanics, 64:383–386, 1997.

C. Thornton and K. K. Yin. Impact of elastic spheres with and without adhesion. Powder Technol., 65:153, 1991.

T. Unger, L. Brendel, D. E. Wolf, and J. Kertsz. Elastic behavior in contact dynamics of rigid particles. Phys Rev E Stat Nonlin Soft Matter Phys, 65(6 Pt 1):061305, Jun 2002.