# Modelling ballast behaviour under dynamic loading.
## Part 1: A 2D polygonal discrete element method approach

Gilles Saussine, C. Cholet, P.E. Gautier, Frédéric Dubois, Claude Bohatier, Jean Jacques Moreau

► To cite this version:

Gilles Saussine, C. Cholet, P.E. Gautier, Frédéric Dubois, Claude Bohatier, et al. Modelling ballast behaviour under dynamic loading. Part 1: A 2D polygonal discrete element method approach. Computer Methods in Applied Mechanics and Engineering, 2006, 195 (19-22), pp.2841-2859. 10.1016/j.cma.2005.07.006. hal-00580793

HAL Id: hal-00580793

https://hal.science/hal-00580793v1

Submitted on 20 Dec 2016

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Modelling ballast behaviour under dynamic loading. Part 1:
A 2D polygonal discrete element method approach

G. Saussine<sup>a,*</sup>, C. Cholet<sup>a</sup>, P.E. Gautier<sup>a</sup>, F. Dubois<sup>b</sup>,
C. Bohatier<sup>b</sup>, J.J. Moreau<sup>b</sup>

Discrete element simulation provides some insight into the alteration of railway ballast after repeated train passings. The present Part 1 is devoted to a 2D model of this granular layer interposed between the deformable ground and the rail sleeper, to which a large number of loading cycles is applied. Ballast grains are modelled as indeformable polygonal solids. A detailed account of the application to this frictional dynamical problem of the Non-Smooth Contact Dynamics numerical method is given. Validation is obtained through comparison with physical experiments performed on assemblies of prismatic mineral grains. Numerical results on the settlement of a track submitted to 20,000 loading cycles or more are presented.

*Keywords:* Discrete element; Cyclic loading; Ballasted track; Contact detection

## 1. Introduction

Ballast is widely used as a constituent of railway tracks because of its mechanical properties and its flexibility in construction and maintenance [1]. This granular material, produced through rock crushing and sifting, achieves the transmission to the platform of static and dynamic efforts induced by the running of trains [2]. Railway companies have to preserve the track quality required by traffic safety and passenger

<footnote><sup>a</sup> SNCF, Research and Technology Department, 45 rue de Londres, 75379 Paris Cedex 08, France</footnote>

<footnote><sup>b</sup> LMGC, UM2, cc048, Place Eugène Bataillon, 34095 Montpellier Cedex 05, France</footnote>
comfort. Following up the track geometry is one of the principal inspection routines since, with time, some geometric *track defects* appear, partly due to the ballast settlement in its interaction with the deformable platform. The problem arises from the non-homogeneity of this settlement after a great number of loading cycles.

The increase of circulation speeds to more than 300 km/h makes still more important to analyse this deterioration process. For a more effective maintenance, one has to identify the factors inducing the appearance of defects.

Since the thickness of the ballast layer is only of the order of 10 grain diameters, *discrete element methods* offer a precious alternative to the continuous model or discrete/continuous approach [3–5] commonly used to study the response of the track to dynamic loading. By calculating the motion of all the grains interposed between the sleepers and the underlying deformable ground, they make it possible to take into account the rearrangements which may occur in the assembly.

Various strategies are available for computing the dynamics of body collections in which the contacts are not necessarily permanent (they constitute unilateral constraints) and affected with dry friction. Similar problems are met in various questions of Civil Engineering, Geodynamics, the Mechanics of Fractured materials and Robotics. The numerical simulation of the interaction of cells in living tissues or that of the constituents of colloids may take inspiration from the ideas put forward in the above domains. Surveys of the numerical techniques used in the treatment of such problems of “Non-Smooth Dynamics” may be found in [6,7].

*Event-driven* methods rest on the division of time into a priori unknown subintervals, on which the list of effective contacts remains unchanged as well as their respective frictional status (sliding or not). On each interval, the same integration techniques as in traditional bilateral dynamics are applied, while watching the evolution of some indicators marking the possible end of the interval. A transition problem then has to be solved in order to determine the contact conditions holding in the consequent interval. The complexity of this transition was discovered by Delassus [8] who showed that in systems with several unilateral contacts, the contacts which break are not necessarily those from which some unfeasibility signals emanate. The event-driven strategy, commonly used in the dynamics of machines, is inapplicable to dense granular collections, due to the large number of bodies and contacts.

In contrast *time-stepping* techniques rest on the use of time-meshes independent of the detected events, generally of uniform length. The most popular ones derive from the work of Cundall [9] or imitate the simulations tools applied in Molecular Dynamics [10], hence the initials MD used to refer to them. They consist in approximating the mechanical constraints of the non-interpenetrability of body pairs by some close-range repulsion laws; the abrupt features of dry friction are similarly smoothed. The resulting dynamics is governed by differential equations tractable by standard numerical integration techniques, practically of the explicit type. The drawback is that, for the sake of precision, very steep repulsion laws have to be used, generating stiff differential equations which require very short step-length for their integration, possibly resorting to the introduction of artificial damping and artificial inertia to secure numerical stability.

A different approach, called (*Non-Smooth*) *Contact Dynamics* method (abbr. NSCD), is invoked in this paper, facing unilaterality and dry friction without smoothing approximation (for its theoretical foundation, see e.g. [11]). It has been implemented in particular in the LMGC90 platform [12] which is expected to become, after some further exchanges between development and tests, a multipurpose set of modelling software involving various sorts of bodies, possibly deformable.

The aim of this paper is to expose some details of a numerical simulation tool derived from LMGC90, specially devoted to the study of ballast behaviour in 2D and 3D under dynamic load. The present first part is restricted to 2D polygonal modelling of ballast grains.

A comparison of some results with experimental measurements obtained at the Laboratoire Central des Ponts et Chaussées, Paris, on a 2D laboratory model of ballast subject to a few thousands of loading cycles, validates the method [13].
The sensitivity of the results to the code parameters is studied in regard to computation cost and precision. Finally we present the first results of investigations concerning the ballast behaviour under repeated cycles of dynamic loading.

## 2. The analytic framework

### 2.1. Parametrization of the system and dynamic equations

If the restrictions imposed to the considered body collection by mutual non-interpenetrability and by external confining boundaries are first disregarded, the possible positions of this system are parametrized through a $n$-vector $q$. The elements of $q$ are the position parameters of the respective bodies when assumed perfectly rigid, but may also include variables describing the possible deformability of some of them, e.g. the node coordinates of a finite-element mesh or the state parameters of a model of the deformable ground.

After that, the geometric effect of non-interpenetrability and confinement is expressed by a finite set of inequalities

$$g_x(q) \ge 0, \quad (1)$$

where each label $\alpha = 1, 2, \dots, \kappa$ corresponds to a possible contact. Inequality holds as equality if this contact is effective. It is assumed in all the sequel that $g_x(q)$ equals the measure of the gap counted along a common normal to the concerned material surfaces, viewed as negative in the case of overlap.

The evolution problem consists in determining the *position function* $t \mapsto q$ of a time-interval $[0, T]$ to $\mathbf{R}^n$. Since the standard differential equations of dynamics are of order 2, one also introduces the *velocity function* $t \mapsto u \in \mathbf{R}^n$. In elementary practice $q(t) = q(0) + \int_0^t u(\tau) \, d\tau$ but, if some elements of $q$ consist of the six position parameters of a three-dimensional rigid body $B$, it is classical to adopt as elements of $u$ the absolute Cartesian components of the velocity of the mass-center $G_B$ and the components of the spin vector of $B$ along principal axes of inertia at point $G_B$. In the latter case, $q(t)$ results from $u$ through slightly more complicated kinematic formulas than the above integration, with the considerable advantage that the contribution of $B$ in the inertia matrix $M(q)$ of the system reduces to a $q$-constant $6 \times 6$ diagonal block. Concomitantly, three of the generalized components of the forces acting on $B$ will consist of the components along the same principal axes of the moment of these forces about $G_B$.

As far as the derivative $\dot{u}$ of the velocity function exists, the dynamics of the system is governed by a differential equation (precisely an integro-differential equation, since $q$ is related to $u$ in integral way) of the form

$$M(q)\dot{u} = F(q, u, t) + r, \quad (2)$$

where the $n$-vector $r$ is made of the generalized components of the—unknown—contact forces, while the $n$-vector $F(q, u, t)$ encompasses the generalized components of all other efforts and also the velocity-dependent terms commonly referred to as *centrifugal* and *gyroscopic*. This dynamic equation has to be combined with some phenomenologic information connecting the unknown contact forces with the local kinematics at contact points. It will be shown in Section 2.4 how this information may be organized in a way which secures (1).

This differential equation does not govern what happens at singular instants such as that of *collisions* or other non-smooth phenomena. This defect is overcome by extending (2) to the setting of non-smooth dynamics. To that end the velocity function $u$, possibly discontinuous, is assumed to have *bounded variation* on the interval $[0, T]$. This entails the existence of an $\mathbf{R}^n$-valued measure on $[0, T]$, the *differential measure* (or Stieltjes measure) of $u$, noted $du$. The function $q$ is still continuous and one reformulates dynamics in terms of an equality of $\mathbf{R}^n$-valued measures called a *measure differential equation* [14–16],

$$M(q) \, du = F(q, u, t) \, dt + ds. \quad (3)$$
Here $dt$ denotes the Lebesgue measure on $[0, T]$ (in fact it equals the differential measure of the function $t \mapsto t$) and $ds$ the differential measure of the $\mathbf{R}^n$-valued function $s$ expressing the cumulated components of the contact impulsions experienced by the system from the initial instant. If a collision occurs in the system at some instant $t_c$, the measure $ds$ possesses an atom of the form $p\delta(t_c)$, where $\delta(t_c)$ denotes the Dirac measure at point $t_c$ while $p$ in $\mathbf{R}^n$ consists of the generalized components of the set of the contact percussions arising in the system at this instant.

In the smooth case, where velocities and contact efforts are sufficiently regular, the classical form (2) is recovered with $du = \dot{u}dt$ and $ds = rdt$.

## 2.2. Time-stepping approximation

One divides the interval $[0, T]$ into subintervals, usually of the same length $h$. Let $]t_i, t_{i+1}]$ (i as “initial”, $f$ as “final”) be one of them. From approximants $q_i$ of $q(t_i)$ and $u_i$ of $u(t_i)$ delivered by the antecedent step, one has to calculate approximants $q_f$ of $q(t_f)$ and $u_f$ of $u(t_f)$. One multiplies both members of (3) by the matrix $M^{-1}(q)$ (assumed to be a continuous function of $q$) and integrates over $]t_i, t_{i+1}]$. At mid-time $t_m = t_i + \frac{1}{2}h$, the value of $q$ is estimated as $q_m = q_i + \frac{1}{2}hu_i$ and these quantities are invoked to approximate the integral of the continuous function $M^{-1}F$ by $M^{-1}(q_m)F(q_m, u_i, t_m)h$. This yields

$$u_f = u_{\text{free}} + M^{-1}(q_m)c, \quad (4)$$

where $u_{\text{free}} = u_i + M^{-1}Fh$ is an approximant of the value that $u$ would take at instant $t_f$ in the absence of impenetrability constraints. The $n$-vector $c$ consists of the cumulated components of the contact impulsions arising in the system during $]t_i, t_{i+1}]$, possibly percussional. Calculating $u_f$ requires to combine this equation with relations conveying some information about the phenomenology of contacts, in order to connect $c$ with geometric and kinematic data. A typical feature of the Contact Dynamics strategy described in Section 2.7 is that the value of $u$ involved in such contact laws is the unknown final velocity $u_f$, so that the resulting numerical scheme is of “implicit” type.

After calculating $u_f$, there only remains to update $q$, for instance $q_f = q_m + \frac{1}{2}hu_f$ in the simple case where $q$ results from $u$ by integration.

As a refinement on what precedes, a “$\theta$-method” is also used [17]: a real number $\theta$ is chosen in interval $[0.5, 1]$ and one takes $t_m = t_i + (1 - \theta)h$, $q_m = q_i + (1 - \theta)u_ih$, $q_f = q_m + \theta u_f h$.

## 2.3. Contact detection

The implementation of the above time-stepping scheme requires that the contacts taken into account in the considered step are identified and geometrically described. In all numerical methods, the detection of contact between two bodies actually consists in observing some—hopefully small—overlap of the portions of space they occupy. It is in such a sense that “contact” will be understood in the sequel. Determining numerically whether two polyhedral bodies in 3D settings, or two polygonal ones in 2D settings, overlap is the subject of a recent book [18]. The treatment of the mechanical interaction requires additionally the identification of a plane (a line in 2D settings) mimicking what, in the case of ideal geometric contact, would be a common tangent plane. Of course, contact may take place through a larger contact zone than a single point.

In the 2D simulations which are the subject of the present paper, the detection of contact between two convex polygonal bodies has been implemented through the Shadow Overlap method devised by one of us, with reliability and robustness tested in several years of previous applications to diverse situations of granular mechanics (see e.g. [19]). Here is the principle.

Let $G_A$ and $G_B$ denote the centers of mass of the convex polygonal bodies $A$ and $B$ and let $V$ be a unit vector such that $G_A G_B \cdot V > 0$ (one starts with $V = G_A G_B / \|G_A G_B\|$). By the “shadow” of a polygon, one means its orthogonal projection onto the direction of $V$, hence the shadow overlap
$$ \mathrm{shov}(\mathbf{V}) = \max\{\mathbf{a} \cdot \mathbf{V} : \mathbf{a} \text{ vertex of } A\} - \min\{\mathbf{b} \cdot \mathbf{V} : \mathbf{b} \text{ vertex of } B\}. $$

If this is found $\le 0$, one concludes that $A$ and $B$ do not overlap. Otherwise, one makes $V$ rotate in the direction of decreasing $\mathrm{shov}(V)$ until producing a change of critical vertex (so are called vertices achieving the max or min above). At this stage, one concludes either to no overlap as above or to the existence of a corner of one polygon penetrating the other. In the latter case, one makes $V$ rotate once again in order to possibly detect a second intruding corner.

If a single corner is found crossing an edge of the partner polygon, the direction of this edge is viewed as the tangent direction. By orthogonally projecting the intruding vertex onto the said edge, one determines the penetration depth or violation of the impenetrability constraint, while the nominal contact point is chosen at half this distance.

In case of double intrusion, the common tangent line is fixed out of some averaging procedure. A segment of this line is identified as the contact segment. In all the numerical experiments reported in this paper the dynamics of such an edge-to-edge situation is treated as involving only two contact points, the extremities of the contact segment. This amounts to imagine the polygon edges affected by an infinitely slight concavity.

Though the above detection procedure to be executed at each time-step is fairly rapid, one wishes to avoid testing all the pairs of grains in the collection. In Molecular Dynamics, where time-steps are commonly hundred times shorter than those used in Contact Dynamics for treating the same problems, it is of crucial importance to save computation effort in the phase of detection. The sophisticated methods of predetection devised in the Molecular Dynamics context [10] are also useful here: the replacement of each grain by a rectangular bounding box, the covering of the numerical domain by rectangular tiles and the book-keeping of the locations of the respective bounding boxes relatively to these tiles.

## 2.4. Global and local levels

The above contact-detection, performed in a configuration $q$ of the system, yields a finite collection of punctual contacts, indexed by Greek labels as in (1). At the contact labelled $\alpha$, the two contacting bodies are conventionally called the candidate contactor and the antagonist one; this distinction is extended by continuity to neighbouring values of $q$ for which the two bodies are not necessarily in contact. The common normal unit vector $\mathbf{n}^\alpha$ is directed toward the candidate body; it is completed by a tangent unit vector $\mathbf{t}^\alpha$ (two such vectors in 3D situations) to make the local frame associated with contact $\alpha$.

In configuration $q$, the kinematic analysis of the way the parametrization has been constructed allows one to calculate, for every imagined value of $u$ in $\mathbf{R}^n$, the velocity vectors, in physical space $\mathbf{E}^3$, of all particles of the system, relatively to the selected reference frame. The relative velocity of the candidate body with regard to the antagonist one equals the vector difference of the velocities of the respective contacting particles, say

$$ U^{\alpha} = H^{\star\alpha}(q)u, \qquad (5) $$

where $H^{\star\alpha}(q)$ denotes a linear mapping from $\mathbf{R}^n$ to $\mathbf{E}^3$.

In view of the definition of generalized components, the contact force $R^\alpha$ exerted by the antagonist body upon the candidate one contributes in the element $r$ of Eq. (2) by the $n$-vector

$$ r^\alpha = H^\alpha(q)R^\alpha, \qquad (6) $$

where the linear mapping $H^\alpha(q)$ from $\mathbf{E}^3$ to $\mathbf{R}^n$ equals the transpose of $H^{\star\alpha}(q)$. The same $H^\alpha(q)$ applies to a possible contact percussion.

When two-dimensional models are concerned as in the sequel of this paper, $\mathbf{E}^3$ should naturally be replaced by $\mathbf{E}^2$.

For brevity in what follows, the label $\alpha$ is dropped.
## 2.5. Unilaterality of interaction

The geometric condition of impenetrability (1) has to be complemented by informations of mechanical nature. First it is put that grain interaction strictly consists in contact with no distant effect, i.e. $R=0$ as soon as the gap $g$ is strictly positive. Secondly it is assumed that no adhesion (gluing effect) takes place, i.e. the normal component $R_N$ is always non-negative. These two statements may be condensed with (1) so as to yield the gap complementarity condition, also called Signorini's condition,

$$g \ge 0, \quad R_N \ge 0, \quad gR_N = 0. \tag{7}$$

At the stage of devising numerical approximations, the time-discretizations of (2) or (3) proves difficult to combine with conditions involving $g$, the risk being to generate unstable algorithms as it frequently happens in MD methods. A similar difficulty is met in the numerical simulation of multibody systems atrained to traditional “bilateral” linkages, a standard situation in the computer-aided conception of mechanisms [20]. The geometric effect of bilateral linkages is then expressed by ordinary, i.e. non-differential equations, to be combined with the differential equations of dynamics yielding what is called a system of algebraic-differen- tial equations. The common approach of such systems consists in applying time-derivation to the non-dif- ferential equations so as to replace them by differential ones.

A unilateral analog to this treatment of algebraic-differential equations rests on the following remark [21,14]: if at some instant $t$ one has $g=0$, the preservation of condition $g \ge 0$ in the further motion implies that $\dot{g}^+$, the time derivative of $g$ on the right of $t$, is non-negative. Classically, $g$ being the normal gap between two bodies, its time-derivative equals the normal component $U \cdot n$, noted $U_N$, with $U$ and $n$ defined in the foregoing. The same equality holds for one-side time derivatives, so that the right-side velocity $U^+$ verifies $U_N^+ \ge 0$. Furthermore, if at instant $t$ one has $g=0$ and $\dot{g}^+ > 0$, then $g$ becomes $>0$ on a nonzero time inter- val on the right of $t$, hence $R$ vanishes on this interval; so does its limit on the right of $t$ noted $R^+$ and in particular the normal component $R_N^+$ of this limit. Finally, if $R_N^+$ exists and is strictly positive, then $R_N > 0$ on some nonzero interval on the right of $t$: on this interval, contact necessarily holds i.e. $g=0$, hence $\dot{g}^+ = U_N^+ = 0$.

This discussion is summarized in the velocity complementarity property

$$\text{if } g(t) = 0 \quad \text{then } U_N^+ \ge 0, \quad R_N^+ \ge 0, \quad U_N^+ R_N^+ = 0. \tag{8}$$

In computation, as well as in the mathematical investigation of the existence of solutions, an amount of violation of the unilateral constraint is tolerated and the definitions of $g$, $n$, etc. are extended in a contin- uous way, so that in these questions one requires of the complementarity (8) to hold more generally if $g(t) \le 0$.

In contrast, if $g(t) > 0$, one has $R^+ = 0$ and no inequality comes to restrain $U^+$.

A converse to (8) is established in [14,22]:

$$\begin{align}
& \text{Assume that } g(0) \ge 0 \nonumber \\
& \text{and that (8) holds for almost every } t \in [0, T]. \tag{9} \\
& \text{Then } g(t) \ge 0 \text{ for every } t \in [0, T]. \nonumber
\end{align}$$

**Remark.** The same reasoning as above shows that if a motion verifies $g \ge 0$ over a time interval and if $g(t) = 0$ at some $t$ interior to this interval, then the left-side velocity at this instant verifies $U_N^- \le 0$. Consequently, if the velocity $U(t)$ exists in the standard two-side sense, the mere kinematics of non- interpenetration requires $U_N(t) = 0$.
## 2.6. Coulomb friction

If contact is assumed frictionless, there only is to complement (8) with the equality $R = R_N n$ to obtain, in conjunction with (2), (5) and (6), a formulation of the dynamical problem.

As for dry friction, its simplest model obeys the law of Coulomb which, among other possible formulations, is classically written in terms of relations involving the normal and tangential components of $U$ and $R$

$$ \bullet \quad \|R_T\| \le \mu R_N, \qquad (10) $$

$$ \bullet \quad \text{if } U_T \ne 0 \text{ then } R_T = -\mu R_N U_T / \|U_T\|, $$

where the friction coefficient $\mu > 0$ may take different values at different contact points.

This formulation refers to the traditional situation of persistent contact, with $U$ denoting the velocity in the elementary two-side sense, so that $U_N = 0$ as precedently observed. But, in Unilateral Dynamics one should at every instant be ready to face the breaking of contact; consequently the information regarding the contact phenomenon should, in our views, be stated in terms of $U^+$, the right-side velocity possibly qualified as the “prospective” velocity (while $U^-$ could be called the “historic” one). For this reason the numerical approximation technique described in 2.7 below will be developed under the following synthetic reformulation of (8) and (10)

$$ \text{For } g > 0, \quad R^+ = 0. $$

$$ \text{For } g = 0: $$

$$ \bullet \quad \|R_T^+\| \le \mu R_N^+ \quad (\text{implying } R_N^+ \ge 0 \text{ provided } \mu > 0), \qquad (11) $$

$$ \bullet \quad U_N^+ \ge 0 \quad \text{and} \quad U_N^+ R_N^+ = 0, $$

$$ \bullet \quad \text{if } U_T^+ \ne 0 \text{ then } R_T^+ = -\mu R_N^+ U_T^+ / \|U_T^+\|. $$

This $q$-dependent relationship between $U^+$ and $R^+$ makes an example of what may generally be called a *contact law* [11]. In spite of its disparate look, it rests on a very consistent convex analytic foundation [21,23]. For other possible contact laws, see e.g. [24].

## 2.7. Implementation

One turns back to the discretized form (4) of the equation of Dynamics. The detection procedure, executed at $q = q_m$, yields the list $J(q_m)$ of the contacts to be treated as active during the time-step. For each of them, relations (5) and (6) are invoked with $q = q_m$, so that the $q$-continuous matrices $H_\alpha$ and $H_\alpha^*$ are treated as constant over the time-step. These relations connect quantities which, in view of Section 2.4, may respectively be classified as *global* (elements of $R^n$) and *local* (elements of $E^3$ or $E^2$ associated with each contact).

It is naturally at the local level that the phenomenologic information regarding the physics of contact is formulated in the terms of contact laws. In the present context the contact laws are that of Coulomb friction, as explicited in (11) for all $\alpha$ in $J(q_m)$, but when coming to discretization the question arises of which value of $U$ should be inserted there. The directing idea of Contact Dynamics is that the (unknown) final value $U_f = H^* u_f$ ($\alpha$ omitted here) mimicks the right-limit $U^+$. Similarly, the value of $R$ involved in the discrete relations is identified with $R^+$. Therefore, the discretization scheme one is constructing will be of implicit type.

Summing up, at each time-step there is to solve the *core problem* which consists in finding values of $u_f$, $r^\alpha$, $U_f^\alpha$, $R^\alpha$, with $\alpha$ ranging in $J(q_m)$, such that
$$u_f = u_{\text{free}} + M^{-1} \sum_{\alpha} h r^{\alpha}, \qquad (12)$$

$$r^{\alpha} = H^{\alpha}(q_m)R^{\alpha}, \qquad (13)$$

$$U_{\text{f}}^{\alpha} = H^{\star\alpha}(q_m)u_{\text{f}} \qquad (14)$$

and that $R^\alpha$, $U_f^\alpha$ verify for each $\alpha$ a contact law of the form (11), generically written as

$$\text{law}_{\alpha}(q_m, U_{\text{f}}^{\alpha}, R^{\alpha}) = \text{true.} \qquad (15)$$

Various strategies may be applied when handling this problem [25,7,26], either privileging the global or the local levels. For instance, one may use the first three equations to obtain

$$U_{\text{f}}^{\alpha} = H^{\star\alpha}u_{\text{free}} + h \sum_{\beta} H^{\star\alpha}M^{-1}H^{\beta}R^{\beta} \qquad (16)$$

an expression to be substituted into (15).

The numerical results presented in this paper were obtained through an iterative solver called “nonlinear block-Gauss-Seidel”. It consists in reviewing the values of $\alpha \in J(q_m)$ cyclically again and again, each time solving a single contact problem with the other contact forces treated as known and consequently updating interactions, until a certain convergence criterion is fulfilled. The criterion is linked with the iterative structure of the algorithm: for each $\alpha$, solving the single contact problem yields the corrective quantities to be added to the local variables in order to comply with (15). This yields an assessment of the precision at which $\text{law}_\alpha$ was satisfied. Diverse ways are open for inserting these local assessments into the management of the iterative procedure. For instance, one may impose on all imperfections to come below a certain threshold (strict Gauss-Seidel policy); alternatively one may trigger the end of iterations after inspecting some average involving all contacts together (non-strict Gauss-Seidel policy), possibly with energetic meaning [24] (Energy criterion). Additionally, in the course of a calibrated computing session, one usually fix an upper bound for the number of iteration so as to avoid wasting computation time when atypical numerical configurations occasionally occur. Fixing also a lower bound for this number has been found to have a favorable effect on the consistency of the time-step succession.

Gauss-Seidel iterations have to be launched from some initial guess of the contact forces $R^\alpha$. One may choose zero but, in the dynamics of dense assemblies like the ballast layer, where the set $J(q_m)$ exhibits only moderate changes from one step to the next, convergence is greatly accelerated by taking as initial guess the contact forces found at the antecedent step for the contacts already active.

All this is compatible with the $\theta$-method variant mentioned in Section 2.2.

## 2.8. Collisions

In the above writing, interactions were described through the contact forces $R^\alpha$. However, when two bodies modelled as perfectly rigid hit each other at some instant $t_c$, their interaction cannot be described in terms of force anymore but in terms of percussion. If they belong to clusters of bodies in contact at instant $t_c$, percussions should also be expected at all contact points.

In the course of a time-stepping computation, collisions are revealed by the presence in the set $J(q_m)$ of contacts which were not considered at the antecedent step. This does not prevent the program from running; there remains to make precise the meaning of the obtained results. The calculation made amounts to admit that, at each examined contact, the contact law (11) is applicable to relate a possible contact percussion $P$ to the right-side relative velocity $U^+$. In view of the fourth line of (11), as soon as $P \neq 0$ the law implies $U_N^+ = 0$, i.e. the collision is treated as completely inelastic.

An algorithmic variant entailing no extra computing cost allows one to model more general situations. It consists in inserting in (11) at the place of $U^+$ some weighted mean of $U^+$ and $U^-$, possibly using different
weights for normal and tangential components. In the simple case of a binary collision, i.e. occurring between two bodies free of any other contact, this procedure is found equivalent to the classical use of (normal and tangential) coefficients of restitution. But, due to the algorithm involving all contacts together, the method still yields plausible results in more general situations [16]. This success should not make one forget that the concept of coefficient of restitution provides only a very rough description of the complex phenomena occurring during the very short time of the collision such as local deformations, elastic propagation of disturbances, etc. Modelling collisions realistically is still an object of research [27,28].

It is a common observation that, even though high restitution could be measured in binary collisions experiments (with steel beads for instance), the impact at moderate velocity of an object against a dense pack produces very little bounces. For this reason, the campaign of numerical experiments reported in this paper has been conducted primarily with zero restitution coefficients.

## 2.9. Deformable layer

The deformability of the track sublayer is essential in the ballast evolution. In fact the whole dynamical behaviour of the structure will depend on the sublayer properties: stiffness, damping, etc. For numerical simulation, one possibility is to discretize the sublayer into two-dimensional finite elements. Since this is time consuming, the following one-dimensional model was used.

One replaces the sublayer by a set of small rigid bodies supported by visco-elastic springs. The visco-elastic law relating the reaction force to the deformation is chosen as

$$R_N = -k_N \varepsilon + v \dot{\varepsilon}, \qquad (17)$$

where $\varepsilon = \frac{L-L_0}{L_0}$ - $\varepsilon_0$ with $L_0$ denoting the length of the spring in the reference configuration, L this length in the current configuration, $\varepsilon_0$ a prestrain, $k_N$ the stiffness of the string and v the viscosity of the “dashpot”. This law involving visco-elasticity may be called *derived Signorini law*, in the sense that some affine change of variables could be used to reduce it to the standard Signorini law.

## 3. Numerical experiments

The object of this part is to assess the pertinence of the numerical method in the case of repeated loading cycles applied to a granular sample. The challenge is to adjust the algorithm parameters in order to compute a significant number of cycles with a reasonable compromise between precision and computer time.

Ballast settlement in real conditions equals a few millimeters after a few hundred thousands of train runnings. The discrete element simulation is expected to yield information about the parameters influencing the settlement amplitude. It is of primary importance to keep under control the numerical errors entailing some violation of interpenetrability constraints. A significant interpenetration of grains would clearly distort the evaluation of settlement.

The whole computations have been undertaken on samples which have the following properties:

* the grains have a pentagonal shape,

* the size distribution of pentagons approximates, at scale one third, the grading of the real ballast. The sample are composed of 361 grains of 1 cm diameters, 242 of 1.5 cm and 121 of 2 cm,

* The volumetric mass is $\rho = 2367 \text{ kg/m}^3$, the friction coefficient is chosen equal to 0.5 and the restitution coefficient to 0 because we consider dense granular packing,

Grains are let to fall down under gravity onto the sublayer, with stiffness $k = 4100 \text{ N/m}$ and tentative viscosity $v = 80 \text{ N/m s}^{-1}$ until they reach an equilibrium state under gravity.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Sample snapshot.</figcaption></figure>

## 3.1. A numerical sample

Comparison with the exploratory experiment presented in [13] will be used as a test of the numerical method. This was an experiment in *Schneebeli* style: pentagonal prisms of uniform length made of high-performance concrete are stacked parallel. The evolution of the prisms is recorded in photographs of one of the end sections of the pile. This makes a physical realisation of a two-dimensional granular material (Fig. 1).

After this preparation, one uses the sample to evaluate the influence of convergence criteria and the influence of step-length on the quality of computation. The tests consist in submitting this sample to vertical loading cycles with a sinusoidal force applied to the sleeper:

$$F = A + B \cos(2\pi ft - \phi)$$

with $A = -1000$, $B = -500$, $f = 20$ Hz and $\phi = -\pi$.

## 3.2. Influence of the convergence criterion

Since the main objective of our simulations is the prediction of ballast settlement, a particular attention should be paid to the study of numerical error coming from the handling of the non-interpenetrability of ballast grains.

So careful preliminary investigations is needed to assess the amplitude of the geometric error resulting from overlap so as to make sure its order of magnitude is inferior to that of the settlement amplitude.

The adjustment of parameters for minimal interpenetration and reasonable computational cost have been made on a representative sample. Computations have been performed on a SUN-Enterprise 450 (450 MHz) unix station (Fig. 2).

### Statistics with energy criterion

<table><tr><td>Norm</td><td>1.66 × 10<sup>-3</sup></td><td>1.66 × 10<sup>-4</sup></td><td>1.66 × 10<sup>-5</sup></td></tr><tr><td>Computation time</td><td>3H57</td><td>7H45</td><td>33H11</td></tr><tr><td>Average number of iteration</td><td>52.8</td><td>186</td><td>1011.8</td></tr></table>

### Statistics with strict Gauss-Seidel

<table><tr><td>Norm</td><td>1 × 10<sup>-4</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td>Computation time</td><td>53H30</td><td>63H53</td></tr><tr><td>Average number of iteration</td><td>2480.7</td><td>2675.9</td></tr></table>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Evolution of cumulated interpenetration (100 loading cycles, 20 Hz).</figcaption></figure>

<table><caption>Statistics for non-strict Gauss-Seidel</caption><tbody><tr><th scope="row">Norm</th><td>1 × 10<sup>-5</sup></td><td>1 × 10<sup>-6</sup></td><td>1 × 10<sup>-7</sup></td></tr><tr><th scope="row">Computation time</th><td>10H58</td><td>15H10</td><td>17H57</td></tr><tr><th scope="row">Average number of iteration</th><td>375.4</td><td>514.9</td><td>664.2</td></tr></tbody></table>

From this results one can do some remarks:

* one has to find a compromise solution between numerical precision and computing time. A too weak criterion leads to a catastrophic solution with an accumulation of interpenetration,

* the different convergence policies do not lead to the same behaviour of the algorithm. The strict policy implies for each time step a large number of iterations. A severe criterion is need to avoid numerical error accumulation,

* it seems that to impose a minimum number of iteration increases the accuracy of the computation. The choice of this minimum number of iteration depends on the load, the influence of strong and weak contact networks.

### 3.3. Influence of step-length

This parametric study was carried out using a sample reproducing an experiment of LCPC [13], but raising the frequency of loading to 20 Hz. For two loading inclinations, 0° and 5°, we have analysed the evolution of the average and maximal interpenetration for 10 different step lengths.

During the loadings, the average interpenetration evolves slightly (Fig. 4). The step-length has a significant influence on the magnitude which is multiplied by 8 if the step-length goes from $1 \times 10^{-4}$ to $1 \times 10^{-3}$. The shorter length provides a more stable evolution of the maximal interpenetration and keeps it smaller.

Maximal interpenetration (Fig. 3) does not influence the evolution of the average interpenetration in the whole of computations. One observes the same trend by applying cyclic loading with an angle of 5°. During the first cycles, one notices a reduction in the average interpenetration which can be justified by the shifting
<figure><img src="image_4.png" /><figcaption>Fig. 3. Evolution of maximal interpenetration.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Evolution of average interpenetration.</figcaption></figure>

of the sleeper on the underlaying grains, as well as a reduction of contact number. Then we find similar behaviour as previously noticed.
The evolution of interpenetration in the particular configuration of our sample, can be controlled
by using a sufficiently small step-length, about $2 \times 10^{-4}$. The maximal interpenetration, sometimes
significant, affects a small number of grains, and is generally located on a contact face-face between
polygons.

**4. Application to track settlement**

The object of this part is to test the accuracy of NSCD discrete element computation in reproducing and
analysing the settlement process. First a comparison between laboratory experiment and numerical compu-
tation is presented. Secondly the influence of the sub-ballast layer stiffness is numerically explored.

4.1. *Comparison with a reference experiment*

A reference experiment has been developed by Combe [13] in order to test the ability of discrete element computation to reproduce the behaviour of a granular medium during repeated loading cycles.

The setting is that of *Schneebeli materials*, often used in Civil Engineering to realize two-dimensional
models of soils. Such materials consist of cylindric or prismatic rods stacked parallel and observed laterally.
In the reported experiment, prisms of pentagonal section are used to mimick the angularity of ballast
grains. The prisms are made of high performance cement with sufficient resistance to endure tests without
damage. Prisms of three different sizes (diameters of 10, 15 and 20 mm) are randomly mixed to form a layer
about 10 cm thick, ending with a free bank (Fig. 5). An aluminium sleeper is put on top; bars of square
section have been fixed to its lower face in order to improve its resistance to lateral sliding. An elastomer
sublayer stands on the frame bottom to model the sub-ballast soil.

Vertical and lateral forces are exerted on the sleeper using a pneumatic jack. A sinusoidal loading, $F = A + B \cos(2\pi ft - \phi)$, with $A = -1000$ N, $B = -500$ N and $\phi = -\pi$, is applied at frequencies of 0.25 Hz to 1 Hz with different inclinations (0–15°). The general test arrangement is shown in Fig. 5. Lateral views of the sleeper and the prisms are captured using a digital camera at regular intervals.

A whole sample was digitized in its initial configuration and was used to compare LMGC90 simulations
with experiments. The digital image did not enable to position the elastomer layer: thus, the first stage was
to fix the position of the polygons that model the elastomer underlayer. The second computation stage was
to let the numerical sample stabilize under gravity. Thereafter, loading cycles were applied. This computa-
tion has been made with a time step of $1 \times 10^{-4}$ s, so that $5 \times 10^{6}$ steps were needed to compute 1000 loading
cycles. It took about three weeks on a Pentium 4 (2.5 GHz) Unix station. The parameters $k_n$ and $v$ of the
sublayer have been experimentally determined by applying an increasing force to the sleeper directly repos-
ing on the elastomer sublayer. From the displacement measured the values $k_n = 4100$ N/m and $v = 80$ N/

<figure><img src="image_6.png" /><figcaption>Fig. 5. Experiments and prisms of pentagonal section.</figcaption></figure>
$m s^{-1}$ were adopted. The friction coefficient between concrete prisms has been estimated from experiment to equal 0.5 and the restitution coefficient is fixed at 0.

The simulation shows a good agreement with the experimental settlement (see Fig. 6) and highlights some outward migration of the grains lying beneath the sleeper. Apart from this zone the grains are only affected with small displacements (Figs. 8 and 9).

The same motions were also observed in the physical experiment and underline the complexity of the settlement phenomenon which depends on the applied loading but also on the evolution of the *contact chains*. In fact the most loaded particles are located directly under the sleeper.

The average grain velocity in the sample decreases as the settlement goes on (Fig. 7).

The results obtained are not in sufficient number to draw a final conclusion but they make one confident in the ability of discrete element method to reproduce the observed phenomena.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Settlement comparison.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Average velocity.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Displacements in the sample.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Displacements in front of the sleeper.</figcaption></figure>

## 4.2. Influence of the sublayer on ballast settlement

The mechanisms involved during settlement are difficult to grasp; furthermore the mechanical or geometrical parameters are ill known. Our aim is to study the sensitiveness of the ballast settlement to the characteristics of track components; first, we are interested in the role of the underlayer, specifically its stiffness.

Eight samples of same characteristics were built but with a different geometrical piling up. This statistical approach of the problem makes it possible to limit the impact of the variations noticed during the experimental phase.

We tested three spring stiffness: $k_N \ge 4100$ N/m (the reference value in the experiment), $v = 80$ N/m s$^{-1}$, then increased by a factor 100 then 1000. The only difference between these samples is the grain dispositions before deposit under gravity. All the computations have been made on a Compaq EV6 (667 MHz) Unix station; computating 3000 loading cycles takes about 3 days (Fig. 10).

Three thousand loading cycles are applied at the top of the sample: the average charge applied is of 1000 N, the signal is sinusoidal, the load amplitude is 500 N and the frequency 20 Hz. The settlement

<figure><img src="image_11.png" /><figcaption>Fig. 10. Example of sample.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Evolution of settlement.</figcaption></figure>

has been evaluated from the evolution of the eight prepared samples for three different stiffness. We repre-
sent the average settlement computed on the same graph (Fig. 11).

As intuitively expected, the more the underlayer is flexible the larger is the settlement (Fig. 11).

At grain scale, we observe that during a loading cycle, the decrease of load magnitude triggers a “breath” phenomenon in the granular medium: the grains can circulate under the sleeper or rearrange themselves because of lower force intensity. High stiffness of the underlayer implies a strong contact force network, more compactness of the pack and less recirculation. Numerical simulation allows one to watch the evolution of the contact force network (the well known “force chains”), giving an orientation for future investigation of the intimate mechanism of settlement.

4.3. Exerting a great number of loading cycles

Ballast settlement is the result of several million of loading cycles, each corresponding to the passing of an axle. Various experiments performed by SNCF reveal two periods in the life of a track. When plotting the magnitude of settlement against time, one observe an initial period with rapid and fluctuating changes, followed by a period of slower and more regular evolution.

The use of numerical simulation for diagnosis and prediction of settlement requires a precise knowledge
of the used model. Presently, the objective is to investigate the behaviour of our model during the applica-
tion of a significant number of cycles.

To this end, we prepared three samples characterized by similar distributions of grain sizes but different
configurations before the deposition under gravity. Following this preparation, we applied a sinusoidal cyc-
lic loading: mean value 1000 N, amplitude of the sinusoidal variation 500 N, frequency 20 Hz. The param-
eters of the sublayer are kN = 4100 N/m, v = 80 N/m s⁻¹.

The settlement evolution for the three samples is different. We observes a brake of slope in the evolution
of the settlement beyond 13,000 cycles that can be explain with a phenomenon of reorganization of the
grains under the effect of the loading (Fig. 12).

After this change, the evolution of settlement is more regular. It should be noticed that in spite of iden-
tical properties, same grading, same density, the final settlement is different for the three samples at the end
of the application of loading cycles. The area where the force network is the strongest remains unchanged
[13], located right beneath the sleeper.

Moreover it is important to notice that the settlement after 20,000 loading cycles of numerical compu-
tations can be compared to quasi-static experiment results (Fig. 13) which have been made with a frequency
<figure><img src="image_13.png" /><figcaption>Fig. 12. Settlement for three different samples.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Comparison with experimental results.</figcaption></figure>

of 2 Hz. The numerical and experimental sample have the same properties and this comparison confirms
the pertinence of the numerical computations presented.

**5. Conclusion**

In numerous studies, ballast is modelled as a continuous medium for which some empirical constitutive laws should be devised, with parameters adjusted through experimental identification. Finite element methods are then applied to solve boundary value problems. Such a modelling should be relevant when considering the whole track system or the sub-ballast structure. To grasp the alteration of the ballast layer in the
medium-to-long term under cyclic loading, we developed an approach at grain scale through a discrete element method based on the Non-Smooth Contact Dynamics numerical strategy.

This required a parametric study in order to determine the way of conducting calculations for a granular medium submitted to loading cycles. The use of a sufficiently small step-length about $2 \times 10^{-4}$ s and an average number of Gauss-Seidel iterations of about 400–500 ensures that the numerical errors of interpenetration do not significantly distort the evaluation of settlement. Available computation speed makes it possible to carry out approximately 1000 cycles of loadings in 24 h with a frequency of 20 Hz.

Comparison between experiment and simulation highlights the effectiveness of this method in reproducing a long process and the evolution of a granular medium subjected to cyclic loading. The experimentally observed settlement and the simulated one are of the same order and the variability observed in the experiment was similarly found in simulations. These first results make it possible to carry out a parametric study to put in evidence the role of the deformable underlayer in the settlement amplitude: the higher is the stiffness the less important is the settlement. Computation on a significant number of cycles, beyond 20,000 cycles highlights the existence of two successive phases in the evolution of settlement, a fact correlated by experimental tests.

The whole of these results demonstrate the feasibility of the discrete element methods in studying the settlement phenomenon. In spite of numerical difficulties it is possible to evaluate the settlement arising in the evolution of a collection of rigid bodies with angular shape and to study micromechanic quantities not currently measurable by experiment. The simulations made show that the low number of grains under the sleeper entails that such traditional average quantities as the stress tensor are insufficient to explain the settlement phenomena. It will be necessary to use tools of statistical physics and to investigate this problem, and to enrich this model by using three-dimensional discrete element method. Currently in progress is a three-dimensional polyhedral approach with the Non-Smooth Contact Dynamics resolution method allowing one to handle 30,000 digitized ballast grains with a view to study for example the lateral resistance of the track [29].

## Acknowledgements

We thank M.M.J. Lux and Ph. Montier from RFF who provide the funding for this research. The development of software devoted to the modelling of ballast by discrete element method has been undertaken by the The Research and Technology Department of the SNCF in collaboration with the Laboratoire de Mécanique et Génie Civil from University of Montpellier.

## References

[1] Direction de la recherche et le technologie. Rail et Recherche, vol. 24, 45, rue de Londres, 75379 Paris Cedex 08, 2002, SNCF.

[2] X. Oviedo, Etude du comportement du ballast par un modèle micromécanique, Ph.D. thesis, LCPC, 2001.

[3] A.S.J. Suiker, A.V. Metrikine, R. De Borst, Dynamic behaviour of layer of discrete particles, part 1: analysis of body waves and eigenmodes, J. Sound Vibr. 240 (May) (2000) 1–18.

[4] A.S.J. Suiker, A.V. Metrikine, R. De Borst, Dynamic behaviour of layer of discrete particles, part 2: response to a uniformly moving, harmonically vibrating load, J. Sound Vibr. 240 (May) (2000) 19–39.

[5] K.Y Zhai, W.M. Wang, J.H. Lin, Modelling and experiment of railway ballast vibrations, J. Sound Vibr. (2003).

[6] B. Brogliato, A.A. ten Dam, L. Paoli, F. Génot, M. Abadie, Numerical simulation of finite dimensional multibody nonsmooth mechanicals systems, Appl. Mech. Rev. 55 (2) (2002) 107–150.

[7] M. Jean, Simulation numérique des problèmes de contact avec frottement, Matériaux et techniques, 1993, pp. 1–3.

[8] E. Delassus, Mémoire sur la théorie des liaisons finies unilatérales, Ann. Sci. Ecole Normale Sup. 34 (1917) 95–179.

[9] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1) (1979) 47–65.

[10] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science Publications, 1987.
[11] J.J. Moreau, An introduction to unilateral dynamics, in: M. Frémond, F. Maceri (Eds.), Novel Approaches in Civil Engineering, Lecture Notes in Applied and Computational Mechanics, vol. 14, Springer-Verlag, 2004, pp. 1–46.

[12] F. Dubois, M. Jean, LMGC90 une plateforme de développement dédiée à la modélisation de problèmes d'interaction, in: M. Potier-Ferry, M. Bonnet, A. Bignonnet (Eds.), Sixième Colloque National en Calcul des Structures, Giens, 2003, Ecole Polytechnique, pp. 111–118.

[13] C. Cholet, G. Saussine, P.E. Gautier, F. Dubois, C. Bohatier, G. Combe, K. Sab, Application of discrete element methods to the modelling of ballasted track, in: World Congress on Railway Research (WCRR), 2003.

[14] J.J. Moreau, Unilateral contact and dry friction in finite freedom dynamics, in: J.J. Moreau, P.D. Panagiotopoulos (Eds.), Nonsmooth Mechanics and Applications, CISM Courses and Lectures, vol. 302, Springer-Verlag, 1988, pp. 1–82.

[15] J.J. Moreau, Bounded variation in time, in: J.J. Moreau, P.D. Panagiotopoulos, G. Strang (Eds.), Topics in Nonsmooth Mechanics, Birkhauser, 1988, pp. 1–74.

[16] M. Jean, J.J. Moreau, Unilaterality and dry friction in the dynamics of rigid bodies collection, in: A. Curnier (Ed.), Contact Mechanics International Symposium, Presses Polytechniques et Universitaires Romanes, 1992, pp. 31–48.

[17] D. Vola, E. Pratt, M. Jean, M. Raous, Consistent time discretization for a dynamical frictional contact problem and complementarity techniques, Rev. Eur. éléments finis 7 (1–3) (1998) 149–162.

[18] G. Van Der Bergen, Collision Detection in Interactive 3D Environments, Elsevier, 2003.

[19] C. Nouguier, Simulation des interactions outils-sol; application aux outils de traitement des sols, Ph.D. thesis, Université Montpellier 2, Montpellier, France, 1998.

[20] E.J. Haug, Computer Aided Kinematics and Dynamics, vol. 1, Allyn and Bacon, Boston, 1989.

[21] J.J. Moreau, Une formulation du contact à frottement sec; application au calcul numérique, Comp. Rend. Acad. Sci., Paris Série II 302 (1986) 799–801.

[22] J.J. Moreau, Some basics of unilateral dynamics, in: F. Pfeiffer, C. Glocker (Eds.), Unilateral Multibody Contacts, Kluwer, Dordrecht, 1999, pp. 1–14.

[23] G. De Saxcé, Z.Q. Feng, New inequation and functional for contact with friction, J. Mech. Struct. Mach. 19 (1991) 301–325.

[24] B. Cambou, M. Jean, Micromécanique des matériaux granulaires, Hermes Sciences, 2001.

[25] M. Renouf, F. Dubois, P. Alart, A parallel version of the non-smooth contact dynamics algorithm applied to the simulation of granular media, J. Comput. Appl. Math. 168 (2004) 375–382.

[26] P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to Newton-like methods, Comput. Methods Appl. Mech. Engrg. 92 (3) (1991) 353–375.

[27] W.J. Stronge, Impact Mechanics, Cambridge University Press, 2000.

[28] D. Stoianovici, Y. Hurmuzlu, A critical study of the applicability of rigid body collisions theory, ASME J. Appl. Mech. 63 (1996) 307–316.

[29] G. Saussine, C. Cholet, P.E. Gautier, F. Dubois, C. Bohatier, J.J. Moreau, Modelling ballast under cyclic loading using discrete element method, in: Proceedings of the Cyclic Behaviour of Soils and Liquefaction Phenomena, April 2004.