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

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.


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

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 a ðqÞ P 0; ð1Þ where each label a = 1,2,. . . , j corresponds to a possible contact. Inequality holds as equality if this contact is effective. It is assumed in all the sequel that g a (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 # q of a time-interval [0, T] to R n . Since the standard differential equations of dynamics are of order 2, one also introduces the velocity function t # u 2 R n . In elementary practice qðtÞ ¼ qð0Þ þ R t 0 uðsÞ ds 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 · 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 _ 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Þ _ u ¼ F ðq; u; tÞ þ r; ð2Þ where the n-vector r is made of the generalized components of the-unknown-contact forces, while the nvector 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 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 R n -valued measures called a measure differential equation [14][15][16], Here dt denotes the Lebesgue measure on [0, T] (in fact it equals the differential measure of the function t # t) and ds the differential measure of the 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 pd(t c ), where d(t c ) denotes the Dirac measure at point t c while p in 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 ¼ _ u dt and ds = r dt.

Time-stepping approximation
One divides the interval [0, T] into subintervals, usually of the same length h. Let ]t i , t f ] (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 f ]. At mid-time t m ¼ t i þ 1 2 h, the value of q is estimated as q m ¼ q i þ 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 where u 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 f ], 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 þ 1 2 hu f in the simple case where q results from u by integration.
As a refinement on what precedes, a ''h-method'' is also used [17]: a real number h is chosen in interval [0.5, 1] and one takes t m = t i + (1 À h)h, q m = q i + (1 À h)u i h, q f = q m + hu f h.

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 Á V > 0 (one starts with V = G A G B /kG A G B k). By the ''shadow'' of a polygon, one means its orthogonal projection onto the direction of V, hence the shadow overlap If this is found 6 0, one concludes that A and B do not overlap. Otherwise, one makes V rotate in the direction of decreasing 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.

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 a, 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 n a is directed toward the candidate body; it is completed by a tangent unit vector t a (two such vectors in 3D situations) to make the local frame associated with contact a.
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 R n , the velocity vectors, in physical space 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 where H wa (q) denotes a linear mapping from R n to E 3 . In view of the definition of generalized components, the contact force R a exerted by the antagonist body upon the candidate one contributes in the element r of Eq. (2) by the n-vector where the linear mapping H a (q) from E 3 to R n equals the transpose of H wa (q). The same H a (q) applies to a possible contact percussion. When two-dimensional models are concerned as in the sequel of this paper, E 3 should naturally be replaced by E 2 .
For brevity in what follows, the label a is dropped.

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, 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 astrained 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-differential equations. The common approach of such systems consists in applying time-derivation to the non-differential 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 P 0 in the further motion implies that _ 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 AE 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 P 0. Furthermore, if at instant t one has g = 0 and _ g þ > 0, then g becomes >0 on a nonzero time interval 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 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 continuous way, so that in these questions one requires of the complementarity (8) to hold more generally if g(t) 6 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]: Assume that gð0Þ P 0 and that (8) holds for almost every t 2 ½0; T . Then gðtÞ P 0 for every t 2 ½0; T . Remark.
The same reasoning as above shows that if a motion verifies g P 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 6 0. Consequently, if the velocity U(t) exists in the standard two-side sense, the mere kinematics of noninterpenetration requires U N (t) = 0. 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 where the friction coefficient l > 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) For 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].

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 a and H Ã a 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 a 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 (a 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 a , U a f , R a , with a ranging in J(q m ), such that and that R a ; U a f verify for each a a contact law of the form (11), generically written as law a ðq m ; U a f ; R a Þ ¼ true.
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 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 a 2 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 a, 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 law a 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 atypic 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 a . 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 h-method variant mentioned in Section 2.2.

Collisions
In the above writing, interactions were described through the contact forces R a . 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 5 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.

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 where e ¼ LÀL 0 L 0 À e 0 with L 0 denoting the length of the spring in the reference configuration, L this length in the current configuration, e 0 a prestrain, k N the stiffness of the string and m 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.

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 q = 2367 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 N/m and tentative viscosity m = 80 N/m s À1 until they reach an equilibrium state under gravity.

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 sinusoï dal force applied to the sleeper: with A = À1000, B = À500, f = 20 Hz and / = Àp.

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

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 · 10 À4 to 1 · 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 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 · 10 À4 . The maximal interpenetration, sometimes significant, affects a small number of grains, and is generally located on a contact face-face between polygons.

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 computation is presented. Secondly the influence of the sub-ballast layer stiffness is numerically explored.

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ð2pft À /Þ, with A = À1000 N, B = À500 N and / = Àp, 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 computation has been made with a time step of 1 · 10 À4 s, so that 5 · 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 m of the sublayer have been experimentally determined by applying an increasing force to the sleeper directly reposing on the elastomer sublayer. From the displacement measured the values k n = 4100 N/m and m = 80 N/ 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.

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 P 4100 N/m (the reference value in the experiment), m = 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   15 has been evaluated from the evolution of the eight prepared samples for three different stiffness. We represent 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.

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 application 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 cyclic loading: mean value 1000 N, amplitude of the sinusoidal variation 500 N, frequency 20 Hz. The parameters of the sublayer are k N = 4100 N/m, m = 80 N/m s À1 .
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 identical 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 computations can be compared to quasi-static experiment results (Fig. 13) which have been made with a frequency 16 of 2 Hz. The numerical and experimental sample have the same properties and this comparison confirms the pertinence of the numerical computations presented.

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 · 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].