Noura Ouhbi
email: ouhbi.noura@gmail.com

Charles Voivret

Guillaume Perrin

Jean-Noël Roux
email: jean-noel.roux@ifsttar.fr

J.-N Roux

3D Particle Shape Modelling and Optimization through Proper Orthogonal Decomposition Application to Railway Ballast

Keywords: Railway Ballast, Shape, Characterization, Modelling, Proper Orthogonal Decomposition, Generation, Stochastic

Based on Proper Orthogonal Decomposition (POD), a new method is presented in order to statistically characterize arbitrary particle shapes using an optimal choice of shape functions identified on a set of 1000 digitized railway ballast particles obtained through 3D Scan. The coefficients of the POD expansion enable a description of ballast grains with varying levels of accuracy. On exploiting the knowledge of their statistical distribution we are able, implementing an appropriate Multivariate Kernel Density Estimation (Multivariate KDE) method, to generate irregular particles with similar morphological features. The description and generation methods are validated by comparing statistical distributions of basic characteristics: surface area, volume, average radius, elongation, flatness, and aspect ratio. Using suitable geometric descriptors defining local curvatures, we identify which surface points might be regarded as forming faces. This shows that the proposed particle generation method is well suited for irregularly shaped granular materials, as a first geometric definition step, before numerical simulations of their collective mechanical properties are carried out by a Discrete Element code dealing with

Introduction

Understanding the mechanical behaviour of the ballasted track [START_REF] Selig | Track geotechnology and substructure management[END_REF] is a major issue for railway companies. As with other granular materials, numerical simulations by the Discrete Element Method (DEM) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Radjaï | Discrete-element modeling of granular materials[END_REF] are a particularly valuable tool, thanks to its unique ability to control independently the micromechanical parameters (e.g. friction or restitution coefficients), and to the insights it provides on granular materials at different scales.

An important question, though, is the modeling of grain shape and of its influence on mechanical properties. As shown in Fig. 1, the ballast grains, with sizes of the order of a few centimeters, assume rather different, irregular shapes. They typically exhibit angular shapes and rough surfaces composed of large quasi planar faces limited by sharp edges. The issue of how to characterize the size and shape of solid particles is currently being investigated in many different fields of applied physical sciences and engineering, including chemical, atmospheric, pharmaceutical and soil sciences [START_REF] Riley | Quantitative shape measurements of distal volcanic ash[END_REF][START_REF] Vallebuona | A procedure to estimate weight particle distributions from area measurements[END_REF][START_REF] Taylor | Some properties of irregular 3 D particles[END_REF][START_REF] Matzl | Stereological measurement of the specific surface area of seasonal snow types: Comparison to other methods, and implications for mm-scale vertical profiling[END_REF][START_REF] Bagheri | Dedicated vertical wind tunnel for the study of sedimentation of non-spherical particles[END_REF]. A number of numerical studies [START_REF] Lim | Mechanics of railway ballast behaviour[END_REF][START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF][START_REF] Azéma | Etude numérique des matériaux granulaires à grains polyédriques: rhéologie quasi-statique, dynamique vibratoire, application au procédé de bourrage du ballast[END_REF][START_REF] Azéma | Packings of irregular polyhedral particles: Strength, structure, and effects of angularity[END_REF][START_REF] Quezada | Penetration test in coarse granular material using contact dynamics method[END_REF][START_REF] Saint-Cyr | Particle shape dependence in 2D granular media[END_REF][START_REF] Yang | Exploring the relationship between critical state and particle shape for granular materials[END_REF][START_REF] Azéma | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF], as well as laboratory experiments [START_REF] Quezada | Penetration test in coarse granular material using contact dynamics method[END_REF], have shown the strong influence of grain shape on granular material behaviour, in both quasistatic and inertial flow regimes. The principal goal of this paper is to provide a quantitative method to model the shape of ballast grains, a necessary step in a research effort based on DEM, towards a quantitative and relevant evaluation of the sensitivity of ballast performance to shape variability. In this paper, we propose a novel statistical method to achieve an optimal and exhaustive particle shape characterization. The input quality is rendered relevant by the millimetric accuracy of 3D scans of real ballast grains. Our approach enables a DEM modelling of ballast grains as polyhedra with, in theory, an arbitrarily high accuracy.

In order to characterize particle shape, various shape descriptors of different levels of complexity have been introduced. As for railway ballast, we can find several shape indicators in the European and French norms (BS EN 13450, 2003). They are categorized in 1D, 2D and 3D descriptors, and are based on simple formulas combining some basic properties such as surface area, volume, or diameter. Such simple descriptors (e.g., elongation or roundness), while providing global shape information, fail to capture the complex morphology of ballast grains (faces,edges, angularity,etc.). More advanced shape descriptors have been introduced in the last years in order to both quantify and represent the realistic geometric characteristics of granular particles before incorporating them in DEM simulations [START_REF] Blott | Particle shape: a review and new methods of characterization and classification[END_REF][START_REF] Williams | Particle shape characterisation and its application to discrete element modelling[END_REF][START_REF] Petrak | Twodimensional particle shape analysis from chord measurements to increase accuracy of particle shape determination[END_REF][START_REF] Bagheri | On the characterization of size and shape of irregular particles[END_REF][START_REF] Jia | Advances in shape measurement in the digital world[END_REF][START_REF] Zhang | Particle shape characterisation and classification using automated microscopy and shape descriptors in batch manufacture of particulate solids[END_REF][START_REF] Yuan | Shape classification of wear particles by image boundary analysis using machine learning algorithms[END_REF][START_REF] Scheibelhofer | Comparing particle size distributions of an arbitrary shape[END_REF].

Most of these approaches are based on computing statistical properties of different shape features, and using them to describe the shape of a set of particles. Sophisticated three-dimensional shape characterizations recently used include 3D extensions of Fourier descriptors [START_REF] Mollon | Generating realistic 3d sand particles using Fourier descriptors[END_REF][START_REF] Mollon | 3d generation of realistic granular samples based on random fields theory and fourier shape descriptors[END_REF], or spherical harmonics (SH) series (which appropriately generalize Fourier series in a spherical topology) [START_REF] Garboczi | Three-dimensional mathematical analysis of particle shape using x-ray tomography and spherical harmonics: Application to aggregates used in concrete[END_REF][START_REF] Grigoriu | Spherical harmonic-based random fields for aggregates used in concrete[END_REF][START_REF] Venkatraman | Potential for protein surface shape analysis using spherical harmonics and 3d zernike descriptors[END_REF][START_REF] Zhou | Random generation of natural sand assembly using micro x-ray tomography and spherical harmonics[END_REF][START_REF] Zhou | Micromorphology characterization and reconstruction of sand particles using micro x-ray tomography and spherical harmonics[END_REF]. They were implemented for shape characterization as well as generation of grain shapes. As explained further in this paper, such approaches, however powerful they may be, prove not convenient for highly angular shapes, as in the case of ballast grains.

Thus, the method we propose is based on the Proper Orthogonal Decomposition (POD) technique [32]. Historically introduced in studies in turbulence as an efficient method to identify and extract coherent structures of a flow [START_REF] Adrian | Stochastic estimation of the structure of turbulent fields[END_REF], POD can be regarded as a way to substitute the Fourier decomposition [START_REF] Lasser | Introduction to Fourier series[END_REF] when the flow directions cannot be assumed uniform or periodic. POD is also known as the Karhunen-Loève Decomposition [START_REF] Karhunen | Zur spektraltheorie stochastischer prozesse[END_REF][START_REF] Loève | Probability Theory; Foundations, Random Sequences[END_REF] or Hotelling analysis [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF], and may be viewed as a variant of Principal Component Analysis [START_REF] Pearson | Liii. on lines and planes of closest fit to systems of points in space[END_REF], or of Singular Value Decomposition [START_REF] Baker | Singular value decomposition tutorial[END_REF]. It is being used through a broad range of applications (e.g., image processing, signal analysis, data compression and optimal control [START_REF] Yin | An efficient data compression model based on spatial clustering and principal component analysis in wireless sensor networks[END_REF]), for which the common objective is to extract the dominant features out of a set of data, and to build a reduced model. As shown in [START_REF] Zhang | Particle shape characterisation and classification using automated microscopy and shape descriptors in batch manufacture of particulate solids[END_REF], this method is also efficient for particle shape classification.

Sec. 2 presents the POD method and identifies its ingredients on a set of scanned ballast grains. Using statistical distributions provided by POD, the numerical generation of sets of equivalent grains is made possible as described in Sec. 3, on resorting to a stochastic method known as the Multivariate Kernel Density Estimation (Multivariate KDE). The method is validated by comparing statistical distributions of some basic shape descriptors (e.g., surface area, elongation or flatness) between the original grain and the numerically generated grain sets. More advanced geometric descriptors are discussed in Sec. 3.3, such as curvature and large faces (identified as low curvature subsets of the grain surfaces), allowing both more physical characterizations and more stringent tests of the modeling of grain morphology. Finally, Sec. 4 reports on a preliminary numerical simulation of some mechanical properties of the model material.

Grain shape characterization

Data Acquisition and Pre-processing

We start from a dataset (provided by SNCF) with n = 1000 point clouds obtained by 3D digitization (3D scan) of n real ballast grain surfaces, each one comprising about 300000 points. Fig. 2 shows 2 examples of these point clouds. To increase POD efficiency we have to maximize the resemblance between the grains. After cleaning the point clouds from possible outliers, the centre of mass of all grains is moved to the origin of coordinates and the grains are rotated so that their principal axes of inertia coincide with a common fixed frame. The resulting grains (see Fig. 3), are represented by the values of their radii, computed, from their centre of mass to their surface, along d = 800 isotropically distributed directions, that are common to all the grains (using the point repulsion algorithm defined in Ref. [START_REF] Albrecht | Curvature-based analysis of point clouds[END_REF]). This relies on the assumption, reasonable for ballast, of star-like grain shapes and uniquely defined radii in each direction. Denoting as r ij the radius of grain i along direction j, we thus define a matrix A, each line of which containing the d radii of one of the n grains, along each one of the d directions. Each one of the d columns of A contains the n radii of the grains as measured in one direction.

We then define the average grain by its radii:

r j = 1 n n k=1 r kj (1 ≤ j ≤ d). (1)
The average grain, as shown in Fig. 4, is close to an ellipsoid with large (L), intermediate (I) and small (S) lengths along principal axes given by: L 29.6 mm; I 25.7 mm; S 21.6 mm.

(2)

Fig. 4 Average grain shape.

POD applied to particle shape

To organize morphological information more efficiently, it is natural to work on the difference between each grain and the average grain, using

∆r ij = r ij -r j . (3)
The radius differences written in (3) define matrix A, and its line vectors ∆r i for 1 ≤ i ≤ n, by subtracting vector r from each line of A. Let us, for any matrix a denote its Frobenius norm as ||a||: ||a|| 2 = i j a 2 ij . ||r|| being the norm of d-dimensional vector r, one has by definition:

||A|| 2 = n||r|| 2 + ||A|| 2 , (4)
and thus the importance of the variation from grain to grain, relatively to the average shape, may be expressed by the following ratio:

δ 1 = ||A|| n 1/2 ||r|| 20% (5)
Applying the POD to the data contained in A consists in carrying out a singular value decomposition of this matrix, i.e.,

A = U • Σ • T V (6)
in which U and V are orthogonal square matrices, of respective order n and d, and Σ is a n × d matrix with all elements equal to zero except the first r ones on the diagonal, σ k , 1 ≤ k ≤ r, which are positive and classified in decreasing order, r being the rank of A.

Equivalently, matrix T A•A, which is square symmetric of order d and contains the correlation coefficients of the column vectors of A, has to be diagonalised. Its d ≥ r eigenvalues λ j are λ j = σ 2 j for 1 ≤ j ≤ r, λ j = 0 for j ≥ r + 1, and, written in decreasing order from the upper left corner, define diagonal matrix D. Matrix V, as defined in [START_REF] Taylor | Some properties of irregular 3 D particles[END_REF], contains the coordinates of the corresponding (normalized) eigenvectors, such that

T A • A = V • D • T V. (7)
As to matrix U, it is easily deduced from equality

U • Σ = A • V, (8)
on exploiting the simple diagonal structure of Σ.

The main interest of the POD in the context of data modeling is its ability to define a set of approximations A (k) for A, each one being optimal in the sense that A (k) minimizes ||A -a|| 2 among all matrices a of rank k.

A (k) is obtained on setting all singular values σ l to zero for l ≥ k + 1 in SVD decomposition [START_REF] Taylor | Some properties of irregular 3 D particles[END_REF], thereby replacing Σ by Σ (k) , containing only k non-vanishing diagonal elements:

A (k) = U • Σ (k) • T V. (9)
A relative error in approximating A by A (k) might be defined as:

(k) = ||A -A (k) || 2 ||A|| 2 = r l=k+1 λ l r l=1 λ l . ( 10
)
By construction, it decreases as k increases and is equal to 0 for k = r (with r = d in the present case, because 800 = d < n = 1000). Conversely, the difference

E(k) = 1 -(k)
, to which we refer as the information content of the approximate morphological description at order k, grows from zero to one. Returning to the vectors (or sets of radius differences) ∆r i , with j coordinates that describe the shapes of the n grains, relation 6 expresses them as linear combinations, with coefficients U il σ l , of a set of basis vectors, hereafter referred to as shape functions or modes, S l , defined by the columns of matrix V, S lj = V jl (the eigenvectors of T A • A). Upon truncating the decomposition at order k one only uses the k first functions, and characterise the shape of the grains by the values taken by k coefficients.

Results and discussion

Fig. 5 shows the error defined in [START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF] as a function of the number k of basis functions. This error measures a relative shape deviation, averaged over all grains, between their real and their approximated shapes (truncated at the order k). Fig. 5 shows that the method is successful in extracting the important information: keeping only the 20 first modes instead of 800 (2.5% of the total number), the error remains below 10%. This performance can be increased as we have larger sets of particles. For more details, see [START_REF] Liu | Spherical harmonic-based random fields based on real particle 3d data: improved numerical algorithm and quantitative comparison to real particles[END_REF].

It is interesting to test the sensitivity of the set of basis functions to the initial number of scanned grains. We reproduced the procedure using randomly chosen subsets of the grain population, and compared the resulting covariance matrices using the method of Ref. [START_REF] Garcia | A simple procedure for the comparison of covariance matrices[END_REF]. The outcome of this study [START_REF] Ouhbi | An accurate description of grain shape for the numerical study of the mechanical behaviour of railway ballast[END_REF] is that the differences become negligible as soon as a population of 400 ballast grains is used in the statistical analysis.

In order to prove the efficiency of our method, we compare its results with those of the well known method shape characterisation method relying on an expansion in Spherical Harmonics (SH). Using the same number of coefficients, N c = 800, corresponding to 100% of the information within the POD approach, we reconstructed our dataset with different shape precision levels and compared with POD results. Fig. 6 shows an example of a grain reconstructed with 800 coefficients using POD and using SH, and one can notice that for a complex shaped particle as ballast, a larger SH basis is required: SH expansions with 800 coefficients fail to construct a grain similar to the scanned one. On trying different basis sizes, we found that the complete precision (100% of the shape information) is approached with 3000-3500 SH coefficients. Fig. 7, plotting SH and POD errors versus the number of coefficients, clearly shows that the POD approach is, by far, the more efficient one. As in the SH analysis, the POD modes of increasing order tend to characterise details of smaller scale, from global shape to surface texture and roughness details. This is apparent in Fig. 8, showing the first mode, containing more than 50% of the morphological information (E(1) > 0.5), to have roughly ellipsoidal shape, while the following ones, with positive and negative lobes, are associated with regions of higher curvatures. The first mode (Fig. 8a) is quite similar to the average grain (Fig. 4), albeit slightly more elongated in shape ????????DE COMBIEN ???? Consequently, the coefficient of the first mode is, to a large extent, characteristic of the grain size, but also of its shape, which is correlated to it. Fig. 9 illustrates how the shape of one grain is gradually retrieved by an increasing number of modes used to reconstruct the shape. Table 1 lists the number of shape functions necessary to achieve the different levels of accuracy illustrated on this figure, as well as in the measurements of characteristic morphological data in the following. Fig. 9 One grain, approximated with the different levels of approximations of Tab. 1, improving from (a) to (f), compared to exact shape (g).

To quantify the evolution of shape descriptors with shape accuracy, we computed surface areas S i , volumes Ω i , and average radii R i , 1 ≤ i ≤ n, defined as follows, from radius values interpolated along directions distributed over the unit sphere such that both angles θ and φ vary by constant steps, ∆θ and ∆φ. ????

S i = ∆θ∆φ d j=1 r 2 ij (θ j , φ j ) sin φ j Ω i = ∆θ∆φ 3 d j=1 r 3 ij (θ j , φ j ) sin φ j R i = ∆θ∆φ S i d j=1 r 3 ij (θ j , φ j ) sin φ j = 3Ω i S i . (11)
Other morphological featured are also investigated. Identifying three linear dimensions of a grain, as length (largest dimension L), width (intermediate dimension I) and thickness (smallest dimension S), as in Ref. [START_REF] Blott | Particle shape: a review and new methods of characterization and classification[END_REF] (which involves embedding the smallest projected surface of the grain onto a plane within a minimum rectangle), one defines the grain flatness as F = S/I, its elongation as E = I/L and its aspect ratio as a = EF = S/L. The shape of mode 1 is such that the larger grains (with large positive coefficient C 1 ) have larger flatness and aspect ratio than the samller ones (those with negative C 1 of larger absolute value).

Probability density functions (PDF) of all these descriptors computed for reconstructed particles with different levels of shape precision are compared with those of the scanned set. Fig. 10 plots the probability density functions (PDF) of S, Ω, R. The ability of the sequence of shape approximation to describe those distributions is more quantitatively appreciated on evaluating the Kolmogorov-Smirnov (KS) likelihood parameter (which quantifies the probability, assuming equidistribution of two sets of data, that the difference in their empirical distributions is as large as the one that is actually measured). Likelihood values for all 6 descriptors are very good, showing that the distributions of morphological features are well reproduced, even with more than 10% error on shape (Tab. 2). The POD modeling scheme thus enables the control of shape precision, but also significant data reduction, as we noticed the very good agreement of shape approximation with truncated information when it comes to simple shape descriptors. This raises the question of comparing descriptors that are representative of a complex shape (faces, edges, angularity etc.). In other words, if the method shows significant reduction, then we should be able to generate grains that are equivalent to the scanned set (in terms of shape), using only a reduced number of modes. According to the obtained KS results, we are able to generate equivalent grains using only 117 modes instead of 800 modes, since it represents 99% of all information. This will be the objective of Sec. 3.

Generating new grains

Procedure

We now exploit the POD procedure output statistics in order to generate distinct and representative 3D ballast point clouds. The distributions of some of the coefficients of the expansion of grain shapes as series of basic shape functions are shown in Fig. 11. By construction, these coefficients are uncorrelated [START_REF] Wikipedia | Correlation and dependence -wikipedia, the free encyclopedia[END_REF] and of zero mean, while their variances add up to ||A|| 2 /n. In order to generate grains with properties similar to the scanned ballast grains from which the basic data set was obtained, one should reproduce those distributions, as well as the dependencies between the different coefficients. From a stochastic point of view, generating statistically equivalent particles means to estimate the multivariate probability density functions (PDF) representing our experimental data. To do so, the parametric approaches propose to search this PDF in an algebraic class of distributions, which are parametrized by a relatively small number of quantities. However, when the dependence structure associated with the unknown PDF is complex, the definition of a relevant parametric class can become very difficult. For instance, assuming for the PDF of the grain coefficients simple representations (such as the multidimensional Gaussian distribution) can lead to the generation of non physical shapes, once the generated coefficients have been applied to the modes.

To circumvent this problem, nonparametric approaches, such as the multivariate Kernel Density Estimation, are generally preferred. This method consists in approximating the PDF of the coefficients as a sum of multidimensional Gaussian PDFs, which are centred at each available value of these coefficients. In that case, the dependence structure that is present in the available data is directly integrated in the modeling of the PDF. How-ever, to correctly concentrate the new generated coefficients on their regions of high probability, and therefore obtain realistic and representative grain shapes, the covariance of these Gaussian PDFs has to be carefully chosen. This is particularly important when the number of available data is not very high compared to the number of coefficients. More details about the method and the algorithm associated with this covariance optimization are presented and discussed in another publication [START_REF] Perrin | Data-driven kernel representations for sampling with an unknown block dependence structure under correlation constraints[END_REF]. The method enables generating as many grains as wanted in some seconds. To be able to compare correctly PDFs with the scanned set, 1000 new ballast grains are thus generated, using 100 shape functions -a limitation motivated by the observation made in Sec. 2 (see Fig. 5) that modes of higher order only contain about 1% of the total information.

3.2 Representativity: simple shape descriptors. Fig. 12 shows some of these simulated grains. Beyond their realistic visual appearance, their representativity is tested with the distribution of the morphologic parameters measured on real grains: volume, surface area, average radius, flatness, elongation, and aspect ratio. The distribution of the three latter quantities is shown in Fig. 13, showing good coincidence between original and simulated grain data.

As a more decisive check, the Kolmogorov-Smirnov test values are above 0.7 for all six quantities Ω, S, R, F , E and a.

Representativity: advanced shape characterization.

While distributions of simple shape parameters appear to validate our model, one essential motivation of the study was the modeling of angular and facetted grains. To assess the performance of the method in this respect, we first have to explain how to identify normal vectors, faces and edges on grains represented by point clouds.

Normal vectors and surface variation at vertices

We implement the method described in Ref. [START_REF] Albrecht | Curvature-based analysis of point clouds[END_REF]. After a selection of a local neighboorhood N p of each point p (9 neighbours per vertex using the KDtree algorithm from Scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]), an eigenanalysis of the covariance matrix of N p is carried out. This amounts to ap-plying PCA to matrix C defined as follows:

C =   p i1 -p ... p i k -p   T •   p i1 -p ... p i k -p   , i j ∈ N p ( 12
)
where p is the centroid of the neighbours p ij of the sample point p. This 3x3 matrix is used to solve the eigenvector problem:

C • v l = λ • v l , l ∈ {0, 1, 2}. ( 13
)
λ l measure the variation of p i , i ∈ N p , along the direction of the corresponding eigenvectors v l . Labelling eigenvalues such that λ 0 ≤ λ 1 ≤ λ 2 , the plane:

T (x) : (x -p) • v 0 = 0 ( 14
)
through p minimizes the sum of squared distances to the neighbours of p. Thus v 0 approximates the surface normal at p. Consequently, λ 0 quantitatively describes the variation along the surface normal. We define the surface variation σ m (p) at point p in a neighbourhood of size m by:

σ m (p) = λ 0 λ 0 + λ 1 + λ 2 (15)
A small value of the surface variation at p signals a well identified apparent tangential plane on the scale of a neighbourhood of m points about point p. Large values of σ m (p), conversely, correspond to large surface curvatures. Given a choice of m (in practice we set m = 9), to each vertex is attributed a normal vector and a value of surface variation. Fig. 14 visualizes the normal vectors computed at each one of the 800 points used to represent one ballast grain, while the surface variation values are shown in Fig. 15. The procedure for face detection is based on a clustering algorithm called DBSCAN [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. A clustering is defining groups of objects (points) sharing the same chosen property. In this case, we chose to define a cluster as a group of neighbouring points with normal vectors satisfying the following condition:

|| → n i - → n j || 2 = 3 k=1 (n i k -n j k ) 2 ≤ 2 (1 ≤ i, j ≤ d) (16)
The chosen value of has a significant impact on how the points will be distributed in clusters. An optimal choice, minimizing the sensitivity of the obtained set of clusters to variations of in some range about the set value, is identified on implementing the HDBSCAN algorithm [START_REF] Campello | Densitybased clustering based on hierarchical density estimates[END_REF]. The best value is = 0.116, which corresponds to a solid angle of 0.17 steradians.

The clustering procedure attributes cluster labels to the vertices of each grain. All vertices belonging to the same cluster abide by the criterion defined in Eq. 16. Some points are not labelled, for lack of neighbours sharing the same normal direction within the set tolerance. In order to interpret the geometrical meaning of these clusters we classify them as follows.

F Large faces, clusters with more than 50 vertices: see e.g. large faces 1, 2, 3 and 4 in Fig. 16.

sF Small faces, clusters with less than 50 vertices.

E Edge vertices, defined as non-labelled ones with a high surface variation (we chose the threshold of 5%,the standard deviation of the surface variation distribution of non-labelled points, whose mean is 0.5%)(Fig. 16).

R Remaining points (Fig. 16), which belong to none of the previous categories.

The result of such a clustering procedure enables the detection of faces and edges, as shown in Fig. 16. Fig. 17 also illustrates this classification, showing edges and remaining points U, i.e., all points that do not belong to large or small faces. The clustering method has shown its ability to detect faces and edges in the same way we would while looking at a ballast grain. To quantify the obtained results, Fig. 18 shows the distribution of the numbers of large faces observed on both scanned and numerically generated grains. The excellent agreement the distribution functions shows the success of the POD method to generate statistically representative ballast grains with the correct facetted, angular morphology. The method also shows its ability to detect faces and edges of a particle even with a high level of shape irregularity.

Next we turn, in a preliminary study, to the scanned and generated grains mechanical properties.

Meshing and Numerical Simulations

A necessary step to carry out DEM simulations of the ballast grains as modelled via the POD approach is to mesh them, so that the point cloud is transformed into a polyhedron with a set of faces and edges. We may then implement the DEM numerical platform LMGC90 [START_REF] Dubois | Lmgc90. 10e colloque national en calcul des structures[END_REF] to simulate mechanical tests.

The meshing procedure should be carried out in a way that modifies the least the morphology of the point cloud. Our practical procedure involves the three following steps. First, we resort to the α-shape method [START_REF] Fischer | Introduction to alpha shapes[END_REF], involving a Delaunay triangulation [START_REF] Lee | Two algorithms for constructing a delaunay triangulation[END_REF], thereby building a mesh on the grain surface, as a set of triangles, whose number is chosen to produce the least error on shape. The second step is a surface subdivision, for which the "Loop S3 method" [START_REF] Boyé | Least squares subdivision surfaces[END_REF] is applied, within the Meshlab software [START_REF] Cignoni | Meshlab: an open-source mesh processing tool[END_REF], resulting in a more regular mesh with many faces. Finally, the third operation reduces the number of faces to the preset value via the "edge collapse" method [START_REF] Kobbelt | A general framework for mesh decimation[END_REF].

Fig. 22 shows how these steps operate on one grain, transforming the initial point cloud into a polyhedron with a number of faces set to 100. Fig. 23 shows how the choice of a number of faces may affect the grain shape. In a first stage, the shapes are further simplified as convex polyhedra that correspond to the convex hull of our non-convex meshed polyhedra, which enables much faster computations in terms of contact detection. This entails a loss of precision on shape, which needs to be quantified. In case a meshing is done with a high number of faces (50000), and thus a better quality, the loss of information (in the sense of the POD error) is about 20% (which would be reduced to less than 5% with nonconvex polyhedra) (Fig. 24). In our case, we use only 100 faces, which induces a larger error.

The evolution of the error restricts the relevance of simulations with large numbers of modes. Here we simulate an assembling procedure via one-dimensional compression between parallel planes with periodic lateral boundary conditions, without gravity. It is carried out using 300 grains, with 100 faces and 3 different shape precision levels after convex meshing: 30% (1 mode), 80% (11 modes) and 100% of information (800 modes). The sample is first generated by depositing the spheres circumscribed to the polyhedron (Fig. 25) in a cell with a constant square section (S = 30 cm ×30 cm). Three values of the inter-granular friction coefficient µ are compared: 0, 0.3 and 0.8. As a reference, a similar simulation is carried out with meshed spheres with 100 faces (and µ=0.3). Contacts with the lower and upper plane walls are frictionless. A constant force of F 0 = 500 N is applied to the upper mobile plane until an equilibrium state is approached, in which the contacts resist further compression (Fig. 25). Figs. 26 to 29 show how the solid fraction and the coordination numbers corresponding to simple (vertex-face or edge-edge), double (edge-face) and triple (face-face) contacts respectively denoted as z 1 , z 2 , and z 3 evolve in time towards their final equilibrium values. We note that the force supported by the bottom fixed plane converges to F 0 after 3.8s. The ratio K= Ec σD 3 of the kinetic energy per grain E c , to the characteristic energy scale σD 3 (σ denoting stress F0 S , and D the average grain diameter) gradually decreases to 4.10 -6 in 10 s. Such time scales are related to the initial state density and to grain and upper wall inertia. The final, equilibrium values of solid fraction ρ, of coordination numbers z 1 , z 2 , z 3 are given in Tab. 3, along with those of x 0 , the proportion of "rattlers" (grains carrying no force), for different numbers of modes and friction coefficients.

A slight increase of solid fraction ρ with the number of modes is observed, while z 1 decreases and z 2 and z 3 increase. This effect of a more accurate, and therefore more complex, shape description, could reflect some evolution in the accomodation of steric hindrance. This shape effect is stronger for low friction coefficients. A strong increase of rattler fraction x 0 should also be noted for growing µ, as in spherical bead assemblies. Fig. 30 plots the distribution of connectivity Z (the fractions of grains having exactly Z contacting neighbours). A significant difference with spheres can be noted, while different shape models yield very similar results: increasing the precision on shape modelling changes the type of contacts without affecting contact numbers. The main outcome of this preliminary mechanical study is thus the observed influence of the accuracy with which the grain shape is modeled on the collective properties of static granular systems assembled by a compression process and equilibrated under prescribed confining stresses. This influence is stronger for low friction coefficients, and tend to fade out for large µ values.

Conclusions.

We propose an innovative method to characterize particle shape, with ballast grains as an application. This method enables shape description with a controlled accuracy, and the generation of new sets of grains, that are morphologically equivalent to real ones, as regards simple geometric characteristics (volume, surface, elongation and aspect ratio) as well as suitably introduced more sophisticated descriptors, such as the number of large faces and the proportion of grain surface area they contain.

Some preliminary numerical simulations, carried out to quantify the impact of particle shape on the mechanical behaviour, reveal that a very significantly reduced number of parameters (11 shape functions) may be sufficient for our model, and also a dominance of friction over particle shape for high friction coefficients. Beyond the simulated assembling stage for a solid granular packing, we expect the shear response to exhibit a higher sensitivity to shape, as observed in [START_REF] Azéma | Packings of irregular polyhedral particles: Strength, structure, and effects of angularity[END_REF] for polyhedra with varying numbers of faces inscribed within the same sphere.

Fig. 1

 1 Fig. 1 (a) Railway Ballast. (b) Two ballast grains on a larger scale.

Fig. 2

 2 Fig. 2 Scanned ballast grains.

Fig. 3

 3 Fig. 3 Pre-processed grains: effects of reorientation (left), of description with 800 directions (right).

Fig. 5

 5 Fig. 5 Error (k), as defined in Eq. 10, versus number k of shape functions.

Fig. 6 A

 6 Fig. 6 A grain reconstructed using (SH) analysis (a) and POD procedure (b).

Fig. 7

 7 Fig. 7 Error (k), as defined in Eq. 10 for POD and SH, versus number N c of coefficients.

Fig. 8

 8 Fig. 8 First shape functions. From top to bottom: mode 1 to mode 4.)

Fig. 10

 10 Fig. 10 Convergence study for distributions of morphological properties (a) surface area, (b) volume, (c) average radius.

Fig. 11

 11 Fig. 11 Distributions of the first POD coefficients.

Fig. 12

 12 Fig. 12 Two generated ballast grains.

Fig. 13

 13 Fig. 13 Distributions of (a) flatness, (b) elongation and (c) aspect ratios in simulated grains, compared to the real ones.

Fig. 14

 14 Fig. 14 Normal vectors at each vertex of one grain.

Fig. 15

 15 Fig. 15 Surface variation at each vertex of one grain, with values encoded as dot size and color.

3. 3 . 2

 32 Clustering procedure, identification of faces and edges.

Fig. 16

 16 Fig. 16 Clustering of normal directions for one grain (different colours represent different clusters).

Fig. 17

 17 Fig. 17Edges E and remaining points U on the surface of one grain.

  Fig. 19 is a plot of distribution, among the population of grains, of the proportion ζ LF i of the d vertices belonging to large faces. We clearly observe a good agreement between scanned and generated distributions. The distributions of surface variation are shown in Figs. 20 and 21 for large faces and edges.

Fig. 18

 18 Fig. 18 Statistical distributions of numbers of large faces for scanned and generated grains.

Fig. 19

 19 Fig. 19 Statistical distributions of the proportion ζ LF for scanned and generated grains.

Fig. 20

 20 Fig. 20 Statistical distributions of surface variation of large faces for scanned and generated grains.

Fig. 21

 21 Fig. 21 Distribution of surface variation of edges for scanned and generated grains.

Fig. 23

 23 Fig. 23 One grain meshed, as in Fig. 22, using different numbers of faces: (a) 4000, (b) 400, and (c) 40.

Fig. 24

 24 Fig. 24 Meshing effect on POD error.

Fig. 22 Fig. 25

 2225 Fig. 22 Meshing steps: (a) point cloud; (b) surface construction; (c) surface subdivision; (d) mesh decimation.

Fig. 26

 26 Fig. 26 Solid fraction vs time (µ=0.3).

Fig. 27 z

 27 Fig. 27 z 1 vs time (µ=0.3).

Fig. 28 z

 28 Fig. 28 z 2 vs time (µ=0.3).

Fig. 29 z

 29 Fig. 29 z 3 vs time (µ=0.3).

Fig. 30

 30 Fig. 30 Connectivity distribution (µ=0.3).

Table 1

 1 Number of modes corresponding to different error levels. The error decreases to zero for 800 modes.

	Error	20% 15 % 10% 5%	1%	0.1%
	Mode number	11	16	22	40	117	350

Table 2

 2 Kolmogorov-Smirnov test results, at different levels of approximation k, for the distributions of volume Ω, surface area S, average radius R, flatness F , elongation E and aspect ratio a.

	k	11	16	22	40	117	350
	KS (Ω)	0.866 0.879 0.984 0.997 0.997 0.997
	KS (S)	0.975 0.984 0.997 0.999 0.999 0.999
	KS (R)	0.725 0.737 0.764 0.785 0.849 0.947
	KS (F )	0.926 0.943 0.965 0.973 0.987 0.998
	KS (E)	0.934 0.941 0.945 0.981 0.986 0.996
	KS (a)	0.887 0.888 0.919 0.928 0.965 0.983

Table 3

 3 Final state characteristics for different numbers of modes N and friction coefficients µ.

	N	µ	ρ	z 1	z 2	z 3	x 0
	11	0	0.690 4.360 3.766 0.320	0%
	800	0	0.697 3.877 0.997 0.427	0%
	1	0.3	0.590 3.046 0.764 0.213	7.6%
	11	0.3	0.594 2.893 1.440 0.280	7.3%
	800 0.3	0.602 2.874 1.593 0.332	5.1%
	1	0.8	0.566 1.893 0.860 0.080 26.6%
	11	0.8	0.567 1.886 0.980 0.080 26.3%
	800 0.8	0.567 1.884 0.981 0.094 22.4%

Acknowledgements This work has been supported by the Lines, Track & Environment Department in SNCF-Réseau, Engineering & Projects.