Imants Svalbe

Shekhar Chandra

Andrew Kingston

Jean-Pierre Guédon

Quantised Angular Momentum Vectors and Projection Angle Distributions for Discrete Radon Transformations

Keywords: Discrete projection, tomography, digital angles, finite Radon transforms

A quantum mechanics based method is presented to generate sets of digital angles that may be well suited to describe projections on discrete grids. The resulting angle sets are an alternative to those derived using the Farey fractions from number theory. The Farey angles arise naturally through the definitions of the Mojette and Finite Radon Transforms. Often a subset of the Farey angles needs to be selected when reconstructing images from a limited number of views. The digital angles that result from the quantisation of angular momentum (QAM) vectors may provide an alternative way to select angle subsets. This paper seeks first to identify the important properties of digital angles sets and second to demonstrate that the QAM vectors are indeed a candidate set that fulfils these requirements. Of particular note is the rare occurrence of degeneracy in the QAM angles, particularly for the half-integral angular momenta angle sets.

Introduction

The ultimate quality with which digital images can be reconstructed from projected views is highly sensitive to the selection of the viewing angles [START_REF] Servières | The mojette transform: Discrete angles for tomography[END_REF][START_REF] Subirats | Angular assessement of the mojette filtered back projection[END_REF][START_REF] Svalbe | Reconstruction of tomographic images using analog projections and the digital radon transform[END_REF]. Conventional CT view angles, Figure 1(a), are constrained by the configuration of the x-ray source and detectors. In contrast, digital image angles are constrained only by pixellation of the discrete array on which the image is to be reconstructed, Figure 1(b). Simply dividing an angle interval into equal or integral steps does not provide descriptive digital angle sets. This becomes even more critical for asymmetric digital images that have one or more elongated axes. True digital angle sets should satisfy the following properties:

1. Generate a set of O(N ) discrete angles for a symmetric N × N array (to balance O(N ) view angles with O(N ) projected elements in each view. For an asymmetric discrete array, far fewer than N angles would be needed). 2. These angles should be constructed in a way that accommodates the integer spacing of pixels. 3. A digital angle set should have the properties of being uniformly distributed over the range 0 • -180 • , or at least be locally uniform over a more limited angle range. 4. Be increasing; any new angle set for N ′ >N should contain all of the previously generated angles for N , with each resulting angle remaining unique.

A digital angle set that satisfies these criteria is derived from the Farey series [START_REF] Hardy | An Introduction to the Theory of Numbers[END_REF]. A different scheme to span analogue real and Fourier spaces to create digital angles was developed through the pseudo-polar Fourier transform scheme of [START_REF] Averbuch | The pseudo-polar fft and its applications[END_REF]. The Farey angles are important here as they arise intrinsically in the Mojette Transform [START_REF] Guédon | The mojette transform: the first ten years[END_REF] and the Finite Radon Transform (FRT) [START_REF] Matúš | Image representation via a finite radon transform[END_REF]. The discrete angle properties of these transforms are reviewed in Section 2.

This paper examines an alternative set of angles designed for use on a discrete grid that are derived from the spatial quantisation of the 3D angular momentum in quantum physics. We will show that the digital angle set derived from the QAM vectors satisfy the above criterion. The idea of using an analogue version of the quantum vector spaces has already been utilised by [8] to represent polar colour variables.

This "naturally" occurring quantum angle set has attractive properties that seem to be relatively unexplored. The quantised angular momenta (QAM) angle distribution has the property of being increasing and appears to be locally uniform, with a slow and smooth decrease in density at larger angles. The set of all possible QAM angles is almost, but not quite, unique. Understanding the formation and distribution of the small number of redundant QAM angles is important when choosing appropriate digital projection and reconstruction angles. This alternative set of quantum-based angles may improve the attainable quality of digital image reconstructions and help optimise the number of projected views required, particularly for the projection of asymmetric digital objects. This paper is organised as follows: a review of the Farey-based digital angles and their relationship to the Mojette and FRT transforms is presented in Section 2. Example sets of QAM angle vectors are compared to similar sets of Farey angles in Section 3. Section 4 describes the origin of integer and half-integer quantisation of angular momentum. Section 5 shows how digital angles can be obtained from these QAM vectors. The remainder of the paper demonstrates the properties of the QAM angles against the requirements 1-4 outlined above; Section 6 demonstrates the angle set has extremely low degeneracy and Section 7 explores the local density and range of the QAM angles.

Existing Digital Angle Schemes

An N × N square (or hexagonal) regular discrete array generates a natural set of projection angles [START_REF] Svalbe | Sampling properties of the discrete radon transform[END_REF][START_REF] Svalbe | On correcting the unevenness of angle distributions arising from integer ratios lying in restricted portions of the farey plane[END_REF]. These unique angles have tangents that are based on the ratios of relatively prime integers belonging to members of the set of Farey fractions F N ranging from 1/N up to 1/1. Taking the arctangent of these fractions produces a set of angles lying between 0 • and 45 2(a). The density of Farey angles is remarkably even and exhibits no degeneracy or replication of a n ya n g l e sa st h ei n t e g e rN increases to infinity. Each Farey angle is defined by a unique vector θ ab that links the origin (0, 0) to a co-prime pair of Cartesian coordinates (a, b). Farey angles ranging from 0 • to 180 • degrees are obtained using (a, b), (b, a), (-a, b)and( -b, a) as four-fold symmetric vectors oriented at θ ab , 90-θ ab , 90+θ ab and 180-θ ab respectively, as shown for F 20 in Figure 2(b). Any N × N image can be reconstructed exactly if N>1+max(|a i |, |b i |) for any set of projections taken at rational angles a i /b i . The Katz criterion [START_REF] Katz | Questions of Uniqueness and Resolution in Reconstruction from Projections[END_REF] ensures that ambiguous "ghost functions" do not exist in projections of the reconstructed image space (and hence that the reconstructed image is unique).

The Mojette transform developed in [START_REF] Guédon | The mojette transform: the first ten years[END_REF] is a generalisation of the FRT [START_REF] Matúš | Image representation via a finite radon transform[END_REF]. Many properties of the FRT have been investigated and applied by Kingston and Svalbe [START_REF] Kingston | Projective transforms on periodic discrete image arrays[END_REF]. The Mojette and FRT formalisms both map between digital images and digital projections. This is done exactly and invertibly, with no interpolation (hence preserving image sharpness) by a deliberate selection of grid-dependent digital view angles and projection paths.

The FRT restricts N to be prime. This endows the projections with the property of minimal information redundancy [START_REF] Svalbe | An image labeling mechanism using digital radon projections[END_REF][START_REF] Kingston | Geometric effects in redundant keys used to encrypt data transformed by finite discrete radon projections[END_REF] and enables the use of very simple projection and reconstruction algorithms. Each 2D square array of prime size p has a pre-determined set of p + 1 rational slopes that define the digital projection orientations. These orientations are a subset of the Farey series for F N , as shown in Figure 2(c). The subset of Farey fractions (a/b) that are selected at each prime size p has its own interesting behaviour, as discussed in [START_REF] Svalbe | Intertwined digital rays in discrete radon projections pooled over adjacent prime sized arrays[END_REF]. The FRT has robust, efficient real-space and Fourier-space reconstruction algorithms based on simple addition. It automatically satisfies the Katz reconstruction criterion. The Mojette Transform allows much more flexibility in the shape and size of the discrete array chosen to represent some discrete object. A set of Farey vectors, tailored to satisfy the Katz criterion, are selected for projection orientations based on the shape and size of the array. Some degree of information redundancy arises in this more general Mojette projection representation. However this redundancy may be exploited usefully in the design of very efficient data transmission, storage and encryption schemes [START_REF] Autrusseau | Mojette cryptomarking scheme for medical images[END_REF][START_REF] Normand | Controlled redundancy for image coding and high-speed transmission[END_REF]. The algorithms to reconstruct images from Mojette projections [START_REF] Guédon | The mojette transform: the first ten years[END_REF] are more complex than for the FRT, largely because of the increased level of redundancy.

Quantised Angular Momentum Directions

The 3D angular momentum vector in quantum physics also generates a "natural" set of discrete angles. The angular momentum vector (j) can only take on values that are integer or half-integer multiples of the reduced Planck's constant (h). When this vector is aligned with respect to an externally imposed reference direction (such as that of the total local magnetic field), the magnitude of the zprojection of j can change only in integer increments, Figure 3(a). This alignment constraint results ultimately from the quantisation of stable energy states for bound particles, and the sensitivity of the state energy to the orientation of the particle orbit.

In atoms and nuclei, the use of quantum mechanics is essential. There the observed magnitude of j ranges from 0 to 8 as electron or nucleon orbits are filled in the atoms and nuclei from hydrogen to uranium. The correspondence principle argues that, in the classical (large angular momentum) limit, the 3D vector j is equally free to take any alignment direction with respect to any z-axis. The spectra of vibrational modes which arise as allowed excitations in finite discrete lattices is another relevant physical example of where a set of discrete angle trajectories emerges naturally from the underlying matrix of the data structure [START_REF] Homolya | Density of states for vibrations of fractal drums[END_REF].

The QAM vectors provide a locally smooth and reasonably uniform global coverage of angle space, particularly over the interval 0 • -45 • , for large angular momentum values, Figure 3(c). At angles closer to 90 • , the density decreases steadily and this property may find application in the limited field of view case encountered when reconstructing SPECT data using cone-beams or for PET/CT tomosynthesis.

The size of the set of QAM angle vectors increases as the magnitude of the maximum angular momentum j is increased, adding 2j + 1 vectors when j -1 increments to j. In the main, the added vectors are new and do not occur at angles generated previously by smaller values of j. There are, however, some values of j that do result in degenerate angles.

The occurrence of such replicated angles in the QAM set is rare, even more so for the half-integer QAM case. In this paper, these relatively rare degenerate angles will be examined with a view to being able to predict and quantify any clumping in the local smoothness of the QAM angle set. Quite localised nonuniformities in the density of angles also occur for the Farey sets [START_REF] Svalbe | On correcting the unevenness of angle distributions arising from integer ratios lying in restricted portions of the farey plane[END_REF]. The analytic work of [START_REF] Boca | Distribution of lattice points visible from the origin[END_REF] has examined and modelled the details of those variations.

Quantised Angular Momentum

Classical angular momentum is a measure of the "turning moment" of a moving object about some axis: it is a vector quantity of magnitude proportional to the radius r of the object from the axis and to the linear momentum p = mv of the object, where v is the velocity of the particle which has mass m. Formally, L = r × p, with the direction of vector L being normal to the plane defined by the vectors r and p. In classical mechanics, the direction of L is free to take any direction in space and |L| is a continuous variable, given that values for r and p are unconstrained.

In quantum mechanics, the measurable momentum and location of all objects are subject to the Heisenberg Uncertainty Principle that reflects the "graininess" of space-time and the particles that can exist within it. Momentum and position cannot be simultaneously specified to a precision below the value of the reduced Planck's constant (h). This automatically imposes an uncertainty on the angular momentum, it being a product of position and linear momentum. This uncertainty means that the angular momentum vector can only be observed to change by amounts proportional to h.

If we take a quantised angular momentum vector L in 3D (x, y, z)spacewit h components L x , L y , L z , then the uncertainty in L means that if the direction of L z is fixed, the vector must "precess" (with unknown phase) around the zaxis so that L x and L y are uncertain, Figure 3(a). The uncertainty in L,a n d the expression for its magnitude, emerges naturally in quantum mechanics after separation of the radial and angular part in the solution to the Schrödinger equation when applied to any particle constrained by a potential well.

In quantum mechanics, all fundamental particles (like electrons) have an internal angular momentum (called spin S). This spin is also quantised, so that changes in spin must also occur in h steps. Spin can be either a symmetric or anti-symmetric component of the quantised particle's wavefunction. Fermion objects have half-integral spin (S =( n +1/2)h)w h i l s tb o s o n sh a v ei n t e g r a ls p i n (S = nh) or positive integer n. It appears that the fermion and boson QAM distributions, whilst being very similar, turn out to have remarkably different angle degeneracy properties. A particle with rotational and spin angular momentum has a total angular momentum j, subject to the same uncertainty and quantisation, with j = L + S.

The length of the quantised angular momentum vector is given by |j| = j(j +1)h.I ti st h ej(j + 1), rather than a j 2 term, that endows interesting properties to the QAM angle distribution. The direction of j is defined by the angle θ measured in cylindrical coordinates relative to the direction of the xy plane. The z-projection of j, j z is constrained to have |j z | = mh,w h e r em is integral for integer j and half-odd integral for half-integer j values.

QAM Vectors as Points in the Plane

If we denote the radius of the projection of the vector j on to the x-y plane as r m ,thenr m = j(j +1)-m 2 , see Figure 4(a). The angle θ of the vector j from the origin to (r m ,m) is given by the arctangent of the gradient, taken here as g(j, m)=m/r m .

The equations for lines of constant |m| form simple parabolas, for example, r m = √ j for |m| = j. The half-integer points fall exactly in between the integer points, as both sets of allowed (r m ,m) values lie on the same quadratic curves as seen in Figure 3(c).

The r m values will not, in general, be integers, but can be scaled and rounded to the nearest integer at any required precision. Here (r m ,m)ha v ethesamerole as the (a, b) integers used in Section 2 for the Farey and FRT examples. The z-projection m is an integer or half-integer but in either case changes in integer steps. As r m may be irrational, the exact process used to convert this value to an integer may be important in any application. For j ≫ 1, j(j + 1) asymptotes to j + 1 2 (and so can never be an integer, nor an exact 1 2 integer). Hence the direction of any QAM vector j cannot ever be aligned exactly along the z-axis, as j z = m must be an integer or halfinteger. The maximum angle of θ occurs when m = j,s ot h a tr m = √ j and g = tan(θ)=j/ √ j = √ j. The maximum angle is 88.7 • for j = 2000 (and reaches only 70.5 • for j =8,foraj typical of the atomic case).

The minimum angle for integer QAM is zero degrees, which occurs when m =0withr m = j(j + 1). The minimum half integer angle occurs at |m| = 1 2 . For large j, the minimum gradient g then approaches 1/(2j). This corresponds to a minimum angle of 1.8 • for j =1 5 /2, the typical maximum half integer j for the atomic case, falling to 0.014 • when j = 3999/2.

QAM Angle Degeneracy

The QAM angle set can adapt to the size and shape of a discrete array, as |j| can be matched to the array size and |m| can be used to accomodate asymmetry in the array shape. The QAM digital angles then meet the design properties 1,2 and 3 that are outlined in Section 1. Are there integer j vectors that link the origin to points (r m ,m)a n d( r m ′ ,m ′ ) that have the same angle? For this to be true, the values j and j ′ (and hence m ′ and m,o rr m ′ and r m )m u s ts c a l e according to

(m ′ /m) 2 = j ′ (j ′ +1) j(j +1) . (1)
The next section examines integer j values where ( 1) is satisfied.

6.1 Degeneracy of the QAM Angles for Integer j

For r m = nm, that is a gradient of g = m/(nm)=1 /n, (i.e. 1 : n), then j(j +1)-m 2 =(nm) 2 means that N = 1+4(n 2 +1)m 2 needs to have integer solutions, where N =2j+1. N also corresponds, incidentally, to the total number of allowed QAM vector projections for a given j.

Then j(j +1) = (n 2 +1)m 2 =( n 2 +1)k 2 l 2 ,w h e r em = kl,w i t hk,a n dl integers. We choose to factor m into a product of two integers, kl, because then we can identify j and j + 1 as separate squared quantities that enable m to be an integer. The values of k and l (and hence j and m)c a nb ed e t e r m i n e d recursively, starting with the values k = l =1.Then j =(n 2 +1)k 2 and j +1=l 2 .

(

The ratio of k/l approximates 1/ √ n 2 + 1 because

l 2 =(n 2 +1)k 2 +1. (3)
The values obtained for k and l turn out to be exceptionally good integer approximations for the irrational number √ n 2 + 1. The next integral solution at j ′ , m ′ turns out to given by

j ′ =(n 2 +1)k ′2 and j ′ +1=l ′2 with m ′ = k ′ l ′ , where (4)
k ′ = nk + l and l ′ =(n 2 +1)k + nl. (5)
The recursive relationship (5) determines all of the redundant solutions for gradients (1 : n). Table 1 shows example redundant (j, m) values for the gradients 1:1and1:3. Finding integers k,l that give integer values of r m has a parallel to the approximation of surds (such as a+ √ b) using continued fractions, where the sequence of continued fraction values is periodic. For example, the value of √ 2-1(whichcorresponds to (n 2 +1)= √ 2forn = 1) can be found as 1/(2 + 1/(2 + 1/(2 + ...), which can be written recursively as a r+1 =1 /(2 + a r )w h e r ea 0 =0 .T h ef r a ctions k/l used to find (j, m) values that each have exactly the same gradient 1:n can also be found using a r+1 =1/(2n + a r ) and adding n to each fraction.

The j values that replicate a given gradient have quadratic separation in j. The gap between the next (j ′ ,m ′ ) with the same slope as (j, m) grows very rapidly. The ratio or scale, s, between successive values of j is given by m ′ /m (or j ′ /j or N ′ /N ) and can be shown to be s =2n 2 +1+2n √ n 2 +1. The ratios for gradients n : 1 have the same values of r m and m in exchanged r o l e sa sf o r1:n, so these angles and their j,m values can be found from the 1:n result. The tabulated N and j values for n : 1 are the same as for 1 : n, but the value of m is just n times that for 1 : n.

For g = p : q (as well as g = √ p/q) there is also a similar sparse redundancy in the representation of angles. Results for g = 3 : 5 are given in Table 2. The total number of redundant angles for the integer QAM case appears to increase approximately linearly with increasing j.F o r0<j≤ 2000, we found just 242 redundant angles in the integer QAM angle distribution out of a total number of j(j +1)/2=2, 001, 000 angles.

Degeneracy of the j QAM Angles for Half-Integer j

For half integer QAM, j and m are both required to be odd, so we write:

(2n + 1)(2n +3)/4 -(2m +1) 2 /4=r 2 m , (6)
where n and m are any positive integers. Then r 2 m =(n 2 -m 2 )+(2n -m)+1/2, so that r m can never have integer values for the half-integer QAM case. For the gradient to be degenerate, i.e.

g = m/r m = m ′ /r m ′ ,then[j(j+2)/4-m 2 /4]m ′2 = [j ′ (j ′ +2)/4 -m ′2 /4]m 2 ,a n d (m ′ /m) 2 = j ′ (j ′ +2) j(j +2) . (7)
Condition ( 7) turns out to be much harder to satisfy than (1) for the integer QAM case. Scaled solutions for j and j ′ that give redundant angles are only possible if j ′ = α 2 j and j ′ +2= β 2 (j +2)with m ′ =( αβ)m,a n dw h e r eα, β are integral, which is similar to the constraint (4) for integer QAM. For 0 <j≤ 3999/2, we found only 16 redundant angles in the first 2 million possible angles (as compared to 242 redundant angles for j up to 2000 for the integer case). Note the repeated occurrence of the integers 845, 2023 and 3969 in Table 3. These integers have highly composite forms, for example 3969 = 3 4 .7 2 . If the half integer and integer angle distributions are pooled (this does not occur in real quantum systems), the number of redundant angles increases, at approximately double the integer rate, to reach a total of 510 redundancies for j up to 2000. The QAM angle set hence satisfies the fourth digital angle property as set out in Section 1. The uniformity of the QAM angle distribution has been examined using the same unevenness criterion (D) as used for the Farey angle set [START_REF] Svalbe | On correcting the unevenness of angle distributions arising from integer ratios lying in restricted portions of the farey plane[END_REF]. The QAM sets are much more uneven than the Farey sets because of the smooth decrease in angle density as the z-axis projected value of j increases. For the half integer case at j = 199/2, D =6 5 .09 over 2203 angles (selected for m<r m and r m <j/ √ 2=7 0 .71 here, from a total of 5050 angles). For the Farey sequence F 86 (1/86, 1/85,...85/86) which has 2273 angles, D =2.106.

Figure 4(b) shows the density of m, r m points from the QAM distribution in the 2D plane which satisfies the property 3 qiven in Section 1. Note that the points "missing" for the low angles are reflections of the positions of the real points at large angles.

Conclusions and Further Work

The QAM vectors produce an interesting set of digital angles with remarkably little degeneracy, especially for the half-integer angular momenta where the conditions required for integer-based solutions are harsher. In physical systems, the QAM directions are limited to either whole or half-integral values, but the above investigation can be extended to include angles for 1/3 (quark-like) or other fractional quantisation values. More work is needed on how to best round or interpolate the (often irrational) values of r m to integers when applying these digital angle sets.

Fig. 1 .

 1 Fig. 1. (a) real space CT projections of a continuous object at three analogue angles. (b) discrete projections of a simple digital object taken at four "grid-friendly" angles.

Fig. 2 .

 2 Fig. 2. (a) Farey angles (F4) made of co-prime integer ratios a/b (a<b ). (b) the set F20 depicted as co-ordinates (a, b) giving θ ab (grey points), and the reflected set (b, a) giving 90 • -θ ab (black points). (c) The first half (0 • <θ ab < 90 • )o ft h eF R Ta n g l e set is a subset of the extended Farey set, shown here for p = 379 and F20.

Fig. 3 .

 3 Fig. 3. (a) Depiction of the allowed quantum angular momentum (QAM) vector precession orientations for angular momentum j = 3. (b) The quantised z-projection of the QAM vectors for j = 5. (c) The full integer j QAM angle set for |j|≤19.

Fig. 4 .

 4 Fig. 4. (a) Definition of the QAM projections m and rm for a given j. (b) Plot of the density of the QAM angle distribution for j = 199/2. Each white point represent a selected angle.

Table 1 .

 1 Examples of degenerate QAM angles for integer j and gradients of 1 : n

		g =1:1	j(j +1)=(2)k 2 l 2
	N	j	m	k	l
	3	1	1	1	1
	17	8	6	2	3
	99	49	35	5	7
	577	288	204	12	17
	3363	1681	1189	29	41
	19601	9800	6930	70	99
		g =1:3	j(j + 1) = (10)k 2 l 2
	N	j	m	k	l
	19	9	3	1	3
	721	360	114	6	19
	27379	13689	4329	37	110

Table 2 .

 2 QAM angle redundancy for the gradient 3:5

		g =3:5	j(j +1)=(34/9)k 2 l 2
	N	j	m	k	l
	35	17	9	9	1
	2449	1224	630	18	35

Table 3 .

 3 Degenerate values of j, m for half-integer QAM angles

	G = m/rm 0.101015 0.127001 0.101499 0.267261 0.316228 0.408248 0.581238 0.588348 0.707107 0.988849 1.224745 1.581139 1.870829 2.12132	G 1/ √ 1/ √ 1/ √ 1/ √ 1/ √ 1/ √ 5/ √ 3/ √ 1/ √ 7/(5 2.7 2 2.31 2.17 2.7 2.5 2.3 2.37 2.13 2 √ 2) Ô 3/2 Ô 5/2 Ô 7/2 3/ √ 2	j, m 9/2,1/2 7/2, 1/2 5/2, 1/2 3/2, 1/2 9/2, 3/2 7/2, 3/2 9/2, 5/2 5/2, 3/2 1/2, 1/2 9/2, 7/2 3/2, 3/2 5/2, 5/2 7/2, 7/2 9/2,9/2	j ′ , m ′ 3969/2, 399/2 2023/2, 255/2 845/2, 143/2 243/2, 63/2 3969/2, 1197/2 2023/2, 765/2 3969/2, 1995/2 845/2, 429/2 25/2, 15/2 and 361/2, 209/2 3969/2, 2793/2 243/2, 189/2 845/2, 715/2 2023/2, 1785/2 3969/2, 3591/2
	7 Distribution of QAM Angles

Acknowledgements

SC holds a postgraduate scholarship from the Faculty of Science and School of Physics at Monash University. AK holds a postdoctoral position supported by a grant from Pays De Loire, France.


Sacha Duverger
email: sacha.duverger@inrae.fr

Vasileios Angelidakis
email: vasileios.angelidakis@fau.de

Sadegh Nadimi
email: sadegh.nadimi-shahraki@newcastle.ac.uk

Stefano Utili
email: stefano.utili@newcastle.ac.uk

Stéphane Bonelli
email: stephane.bonelli@inrae.fr

Pierre Philippe
email: pierre.philippe@inrae.fr

Jérôme Duriez
email: jerome.duriez@inrae.fr

Yukio Nakata

Shuji Moriguchi

Shin- 1359 Taro Kajiyama

Ryunosuke Kido

Hidetaka Kikkawa

Daiki Saomoto

Takano

Michael Rackl

Kevin J Hanley

Thomas Roessler

André Christian Richter

Investigation techniques and physical aspects of the angle of repose of granular matter

Keywords: Angle of repose, Discrete Element Method, Non-spherical particles

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

 [START_REF] Matuttis | Discrete element simulations of dense packings and heaps made of spherical and non-spherical particles[END_REF]and possible geometrical effects [START_REF] Zhou | Numerical investigation of the angle of repose of monosized spheres[END_REF][START_REF] Miki Y Matsuo | Geometric effect of angle of repose revisited[END_REF]. The latter comes in addition to the more natural influence of physical microscopic properties such as particle shape [START_REF] Pöschel | Static friction phenomena in granular materials: Coulomb law versus particle geometry[END_REF][START_REF] Hao Chen | Dem investigation of angle of repose for super-ellipsoidal particles[END_REF] and contact friction [START_REF] Pöschel | Static friction phenomena in granular materials: Coulomb law versus particle geometry[END_REF][START_REF] Zhou | Numerical investigation of the angle of repose of monosized spheres[END_REF][START_REF] Hao Chen | Dem investigation of angle of repose for super-ellipsoidal particles[END_REF], as well as mesoscopic ones, such as fabric [START_REF] Hao Chen | Dem investigation of angle of repose for super-ellipsoidal particles[END_REF]. Last, it should be noted that granular heaps may not systematically conform to a linear slope [START_REF] Topić | Nonuniformities in the angle of repose and packing fraction of large heaps of particles[END_REF][START_REF] Akbar | The slope of dry granular materials surface is generally curved[END_REF] which may prevent one to define a single-valued AOR.

Following up on these previous works, the aim of the present manuscript is twofold. First, rigorous simulation and measurement methods are proposed in order to ease evergoing AOR studies. Second, with the help of these methods, an in-depth study is conducted in order to gain further insights on the AOR variations with respect to physical parameters. The present analysis combines the use of two Discrete Element Method (DEM) approaches and existing experimental results recently proposed by the Japanese Geotechnical Society (JGS) as part of a round robin series of tests [4,[START_REF] Saomoto | Round robin test on angle of repose: Dem simulation results collected from 16 groups around the world[END_REF].

The remainder of the article is as follows.

Section 2 first recalls the JGS reference experiments [4] and the two DEM formulations which are both executed within the YADE code [START_REF] Smilauer | Yade Documentation 3rd ed. The Yade Project[END_REF] but differ in the way the exact shape of the grains is described. Section 3 then introduces new methods enhancing AOR studies, namely a systematic definition of the AOR value after detection of the external slope and versatile measurement methods of the packing compacity (void ratio) of the heap since the latter is a fundamental property of granular matter. It also provides a discussion on computational aspects of the two DEM approaches used to simulate the same JGS experiments. Section 4 finally provides new insights on the role of some physical parameters on the AOR value, after conducting a large number of DEM simulations interpreted in a statistics fashion for the sake of robust conclusions. Mechanics and Geotechnical Engineering (ISS-MGE) [4] and will serve here as reference. In a first step, data only included an experimental characterization of the granular material at hand, together with properties of the two experimental setups used for AOR measurement, before that blind DEM predictions of the AOR values could be proposed by international participants to 89 the round robin and compared with experimental 90 values [START_REF] Saomoto | Round robin test on angle of repose: Dem simulation results collected from 16 groups around the world[END_REF].

91

An artificial granular material was considered 92 with non-spherical particles made of 3D-printing 93 resin. Particles constituting the mono-dispersed 94 material resemble a tetrahedral arrangement of 95 four spheres clumped together (see Figure 1). Indi-96 vidual spheres have a radius of r s = 0.3101 cm, 97 while each global particle is inscribed in a radius 98 r clump = 0.5 cm. spheres is created to reflect the particle morphol-130 ogy as a so-called clump, e.g. as shown in [START_REF] Angelidakis | CLUMP: A Code Library to generate Universal Multi-sphere Particles[END_REF]. properties of a clump, many DEM codes still simply add the masses of the clump members and directly combine their inertia matrices, which leads to an overestimation in the case of clumps with overlapping members, like the one adopted to simulate the present 3D-printed particle. To mitigate this issue, methods to adjust the density of each sphere-member have been proposed in the literature, such as the one of Ferellec and McDowell [START_REF] Ferellec | A method to model realistic particle 1433 shape and inertia in DEM[END_REF] to correct mass and inertia at the cost of some pre-processing efforts. YADE, along with PFC, provide an alternative solution, where a three-dimensional grid of voxels is generated in the bounding box of the particle, and it is evaluated for each voxel whether it belongs to at least one sphere-member of the clump. For the particles at hand in this study, a grid size of 1000×1000×1000 voxels is used to estimate the volume (and thus the mass) and inertia tensor, with negligible discretisation error induced by the grid resolution, since finer grids led to the same inertial properties.

Potential particles

While the above clump approach is a straightforward DEM strategy for describing the physical particles at hand (Figure 1), a comparison is carried out with a second approach using the so-called "potential particles" introduced by Houlsby [START_REF] Tinmouth | Potential particles: 1436 a method for modelling non-circular parti-1437 cles in DEM[END_REF], and extended to three-dimensions by Boon et al. [START_REF] Weng Boon | A new contact detec-1441 tion algorithm for three-dimensional non-1442 spherical particles[END_REF]. The potential particles are generalised convex non-spherical particles, assembled as a combination of 2 nd degree polynomial functions and a 167 fraction of a sphere, while their edges are rounded 168 with a user-defined radius. In line with their inher-169 ent restriction to convexity, rather common in 170 DEM with complex shapes, e.g. as in [START_REF] Matuttis | Discrete element simulations of dense packings and heaps made of spherical and non-spherical particles[END_REF], the For the exact definition of a potential particle, as detailed in Boon et al. [START_REF] Weng Boon | A new contact detec-1441 tion algorithm for three-dimensional non-1442 spherical particles[END_REF], a set of N planes are assembled such that their normal vectors point outwards, with their interior forming a convex polytope. These planes are summed quadratically and expanded by a distance r, which is also related to the radius of the curvature at the corners.

V cyl = 6.67 • 10 -4 m/s H cyl = 9 • 10 -2 m R cyl = 8 • 10 -2 m Plane strain N part = 2, 150 V par = 4.3 • 10 -2 m/s H par = 1.9 • 10 -1 m L par = 1 • 10 -1 m
Furthermore, a 'shadow' spherical term is added,

183
where R is its radius and 0 < k ≤ 1 denotes the 184 fraction of sphericity of the particle. A value of 185 k ≈ 0 corresponds to a nearly sharp polyhedron,

186
while k = 1 corresponds to a perfectly spherical 187 particle.

188

A potential particle is eventually defined by a 189 potential function f as in Equation 1:

190 f (x, y, z) =(1 -k) N i=1 ⟨a i x + b i y + c i z -d i ⟩ 2 r 2 -1 + k x 2 + y 2 + z 2 R 2 -1 (1)
where (a i , b i , c i ) is the normal vector of the i th Here, the mathematical formulation of the potential particles enables one to approximate the given particle shape by a rounded tetrahedron. To decide which planes to use in order to assemble the potential particle of the 3D-printed material, two criteria were considered, a physical and a practical one, with the latter aiming to achieve post-processing convenience: (1) First, the potential particle should capture the morphology of the physical particle as faithfully as possible in terms of size, surface curvature, mass and inertia of the given physical particle, or other shape descriptors such as the sphericity; (2) To achieve comparable results with the clump models, for the evaluation of the AOR, it is convenient for each potential particle to be monitored via four points being located at the same positions than the centers of the four spheres making the tetrahedron. Thus, it is sought that the potential particle has a straightforward analogy to this format. To satisfy these criteria, dimensions, while also approximating its curvature. The parameters r, R and k were chosen via a trial-and-error procedure. Figure 4 demonstrates visually the geometrical faithfulness of the generated potential particle to the shape of the real, physical particle.

In addition to modelling the rounded, tetrahedral-like particles, the potential particles also serve to simulate cuboidal elements of various sizes, making the moving and still parts of the plane-strain and axisymmetric devices, enabling one to build YADE models using a single, unified approach and contact detection algorithm.

Particle shape characterisation

As demonstrated in Figure 4, the selected potential particle can approximate the morphology of the physical particle faithfully, as it qualitatively represents the main dimensions of the particle, determining particle form, along with the curvatures of its edges/corners, relating to particle roundness. However, the potential particles modelling approach cannot represent the concavity of the physical particle. A quantitative characterisation of particle form was also performed using SHAPE [START_REF] Angelidakis | SHape Analyser for Particle Engineering (SHAPE): Seamless characterisation and simplification of particle morphology from imaging data[END_REF], an open-source shape analysis software for three-dimensional particles, in order to quantify in Table 3 the similarity between the physical particle and its two numerical replicates. To this end, the surface mesh of the physical particle was first tessellated from its corresponding DEM Investigation techniques and physical aspects of the angle of repose of granular matter Bringing the error on mass down to zero through this scaling, the error in inertia values for potential particles drops from 27.17 % down to 7.96 %. It is interesting to note that using overlapping spheres with no correction for uniform

DEM contact formulation

At each contact, kinematics is defined with the normal and tangential relative displacements of the particles, u n and u t respectively. For the clump model, contacts are detected between spheres belonging to different clumps and u n is computed as the norm of the branch vector to the spheres' radii, while u t is computed incrementally, see e.g. [START_REF] Duriez | Pre-1457 cision and computational costs of Level Set-1458 Discrete Element Method (LS-DEM) with 1459 respect to DEM[END_REF]. For the potential particle model, u n is computed using a bracketed line-search algorithm as detailed in Boon et al. [START_REF] Weng Boon | A new contact detec-1441 tion algorithm for three-dimensional non-1442 spherical particles[END_REF], deployed along the contact normal direction and starting from the contact point, to detect two points on the surface of each particle, forming a branch vector, the norm of which is considered as the sought approaching distance. The shear increment of u t is calculated in a similar manner as for spheres, i.e. via time integration of the shear component of the relative velocity during contact.

The same contact model applies to these kine-346 matic quantities for both the clump and the 347 potential particle approaches, accounting for lin-348 ear visco-elasticity and friction (Figure 5). In the normal direction, a spring with a normal stiffness 350

K n is associated in parallel with a viscous damper 351 of coefficient c n , as formulated in Equation 2. In 352 the tangential direction, a spring with a tangential 353 stiffness K t is associated in series with a frictional 354 slider (contact friction angle φ), see Equation 3.

F n = max(K n u n + c n un , 0) (2)
|F t | = min(K t |u t |, F n tan(φ)) (3)
One should note that different YADE classes implement the above Eqs. 2-3 for clumps and potential particles with different methods of expressing the viscous damping coefficient c n . In all cases, a desired normal restitution coefficient e n serves a starting point before some differences appear in the YADE workflow, as detailed in Appendix A. Nevertheless, Figure 6 illustrates the common dissipative behavior of both models with the same F n (u n ) curves in the case of two colliding spheres (obtained after using k = 1 in Equation 1for the PP approach) with an initial relative normal velocity V , demonstrating the consistency of the two implementations of visco-elasticity. In the framework of the round robin test, the JGS measured the contact friction angle φ and the normal restitution coefficient e n for resin against acrylic contacts and for resin against resin contacts, as well as the normal stiffness K n for resin spheres. Experimental measurements exhibited a variability and are thus given as distributions (see Figure 7).

409

The default set of parameters for this numeri-410 cal setup is the experimental one previously given 411 in Table 1. spherical members, the increased effort in computational time when using potential particles is in a classical order of magnitude for DEM approaches for non-spherical particles [START_REF] Duriez | Pre-1457 cision and computational costs of Level Set-1458 Discrete Element Method (LS-DEM) with 1459 respect to DEM[END_REF].

A systematic determination of the angle of repose

This section proposes two rigorous methods to measure the AOR, first, by defining an outer surface of particles and second, by computing an angle from these particles positions.

Outer surface detection

In the axisymmetric case (respectively plane strain case), the 3D space is discretized in several subdo-

mains {r; θ ∈ [θ a , θ b ]; z ∈ [z a , z b ]} (respectively {x ∈ [x a , x b ]; y ∈ [y a , y b ], z})
, giving an intersection with the outer surface at max(r) (respectively max(z)) in each subdomain. The extent of each interval is selected such that only one particle 458 should be therein detected as belonging to the 459 outer surface. For such a purpose, length scales L η 460 are used for the coordinates θ, z in the axisymmet-461 ric case and x, y in the plane strain case. The index 462 η can represent each of these coordinates. The 463 number of intervals on each coordinate is then:

464 N η = L η d clump -1 (4) with L θ = 2πR cyl , L z = H cyl , L x = L par , and 465
L y = H par . 466 θ i a = i 2π N θ ; θ i b = (i + 1) 2π N θ (5) with i ∈ 0; N θ η i a = i L η N η ; η i b = (i + 1) L η N η (6) with i ∈ 0; N η , η ∈ {x, y, z}

Angle of repose measurement 471

From this point the method is the same in both Assuming a z-invariance of the heaps, we project the spheres on the ( x, y) planes (see Figure 10) and perform a linear regression on the resulting points to determine the AOR α. Letting the linear regression be y 1 = a 1 x + b 1 , one has:

α = arctan(a 1 ) (7)
Consistently with [START_REF] Topić | Nonuniformities in the angle of repose and packing fraction of large heaps of particles[END_REF][START_REF] Akbar | The slope of dry granular materials surface is generally curved[END_REF], one can notice that the surface isn't exactly flat but slightly curved (especially in the axisymmetric configuration). It can thus be useful to compute a second degree regression as well in order to fit the outer surface in the best possible way. Letting the second degree regression be y 2 = a 2 x 2 +b 2 x+c 2 , one can compute a local angle:

α( x) = arctan(2a 2 x + b 2 ) (8)
To make the measurement more meaningful, one may naturally restrict the procedure to a smaller zone of the heap: the particles considered in the regressions would only be the ones inside an interval [ x min , x max ]. Indeed, the lower particles may be abruptly blocked by the bottom ridge of the container devices and should be excluded from the measurement. Also, particles with a high

x, away from the opened boundary, could be unaf- [START_REF] Miki Y Matsuo | Geometric effect of angle of repose revisited[END_REF] which, considering Equation 7, gives the standard deviation on the angle, StD(α):

fected
StD(a 1 ) = 1 N -2 N i=1 (a 1 x i + b 1 -y i ) 2 N i=1 ( x i -x) 2 (9)
StD(α) = StD(a 1 ) 1 + a 1 2 (11)
Figure 10 shows the regressions made on the projection of the outer surface in both configurations and the resulting angle for x max / x box = 0.4, The error on the measurement is very low but increases with x max . The measurement is more stable for low x max , specially in the axisymmetric case. From now on, the measurements will be performed on most of the outer surface, using

with x box ∈ {R cyl , L par }.
x min = 0.32d clump and x max = x box .

Error due to repeatability

The simulations performed with the clump model include two sources of randomness. The first one is the initial configuration of the sample, with random positions for the particles in the initial cloud. In order to quantify the repeatability error, a series of simulations was performed with the clump model using 30 different values for the seed parameter, the particles in the initial samples of each simulation thus have different positions and contact properties. This series will be called Investigation techniques and physical aspects of the angle of repose of granular matter CLP1 and uses the default parameters of Table 568 4. Among those parameters, the time step is 569 computed from contact stiffnesses and particle 570 masses following [START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soilinclusion problems[END_REF]. Note that all samples have 571 approximately the same initial densities. Figure • the tetrahedra method applies for any geometry of sample with no requirements on the geometry. It is based on a triangulation of the sample.

• the sub-volume method requires to define an homothetic sub-volume inside the sample, 659

V tet part = Nmc i=1 H {p} cut χ p (x i ) N mc × V tet (12)
A local void ratio can then be computed for 660 the tetrahedron:

661 e tet = V tet -V tet part V tet part ( 13
)
And globally:

V part = {tet} V tet part (14)
V tot = {tet} V tet (15) e = V tot -V part V part (16)
Taking advantage of the independence between In the case of the plane strain configuration, the homothetic sub-volume can be determined as follows (Figure 13): In the case of the axisymmetric configuration, χ sub (x) function, the number of vertices inside the sub-volume for a particle p reads:

s min = (1 -C)(S max -S min ) + S min (17)
s max = C(S max -S min ) + S min (18
N p in = 8 i=1 χ sub (x p i ) (23)
If N p in = 8, the particle p is completely inside the sub-volume while if N p in = 0 the particle p is completely outside the sub-volume.

Denoting V p the volume of the particle p, the total volume of particles completely inside the sub-volume is: determine if they are simultaneously inside the volume then reads:

V in part = {p | N p in =8} V p (24
736 V part = V in part + V cut part (27)
Finally, the void ratio is determined by:

737 e = V sub -V part V part (28)
For the simplest sub-volume geometries the expression of V sub is trivial. In more complex situations it can be determined using once again the Monte Carlo method inside the sub-volume bounding box:

V bb sub = (x max -x min )(y max -y min )(z max -z min ) (29)
V sub = Nmc i=1 χ sub (x i ) N mc × V bb sub (30)
Both methods can be optimized when used with simple shapes (e.g. spheres): one could detect more precisely which particle may be cut. Also, one may be able to draw uniformly points directly inside the particle instead of the bounding box, making it possible to set aside the function χ p and giving a more accurate Monte Carlo method.

Parallel implementation of the tetrahedra method

The independence of the processing of each tetrahedron makes it possible to parallelize this method. A series of measurements was performed on 30 clump samples at their initial states (showing different individual locations of particles) T Ncores ± ∆T Ncores as follows:

S = T 1 T Ncores (31) ∆S = S ∆T 1 T 1 + ∆T Ncores T Ncores (32)
Since the CPU cache was not precisely con-784 trolled, the total CPU load had an influence on the 785 computation speed, which might lead to a speed-786 up seemingly above perfection in the eventuality of 787 on 1 core and not for more cores.

V tot = (Xmax -X min )(Ymax -Y min )(Zmax -Z min ) (33) V part = {p} V p ,
with {p} the set of all particles (34) 5).

e REF = V tot -V part V part (35)
Results are given in Figure 18 for what concerns the initial and final states of the samples.

Most importantly, the two different values of tangential stiffness are shown to result in virtually the same AOR distribution. The K s = 240N.m -1 value will thus be kept in the remainder of the sequel for it results in a higher critical time step.

One may furthermore note that the initial coordination number is slightly lower with a higher K s , which is expected since stiff particles tend to be further away from each other, even when constrained. However, at the final state, the average coordination number is unaffected by K s , certainly because they are not constrained enough for their relative distance to depend on K s .

(Non-)Sensitivity to the particle mass density

While the AOR α refers to a static condition, the mass density of particles ρ physically affects the prior dynamic evolutions of the system. On the other hand, from a computational standpoint, the density also controls the critical time step of the

894

Corresponding parameters are all given in Table 6. 5). 100 kg/m 3 0.652 ± 0.011 10 500 kg/m 3 0.633 ± 0.010 1, 111 kg/m 3 0.623 ± 0.009 2, 000 kg/m 3 0.618 ± 0.010 4, 000 kg/m 3 0.607 ± 0.011 8, 000 kg/m 3 0.603 ± 0.012 10, 000 kg/m 3 0.603 ± 0.010 Axisymmetric 2, 468 1, 111 kg/m 3 0.743 ± 0.022 4, 000 kg/m 3 0.723 ± 0.015 8, 000 kg/m 3 0.709 ± 0.019 10, 000 kg/m 3 0.694 ± 0.010

Numerical angle of repose vs experimental one

The numerical simulations are now compared with the experimental results provided at the end of the JGS round-robin and in [4]. In this framework, a simpler method was adopted to compute the AOR, considering only the highest particle instead of the whole external surface as in previous Section 3.2, for sake of simplicity during the experiments. In the axisymmetric configuration, slopes are actually determined in 360 directions being not exactly radial and their average is used to compute the AOR, while in the plane strain configuration the AOR is computed using only one slope direction in the plane.

In this subsection, the exact same method is adopted to interpret our numerical results for a consistent comparison. The set of parameters used for the clump model is the same as for CLP1 (see Table 4), and the set of parameters used for the potential particle model is given in Table 7.

Table 8 compares the obtained experimental 938 and numerical results. In the plane strain config-939 uration, the experimental AOR is approximately 940 8% higher than the one obtained for the clump 941 model and 16% higher than the one obtained for 942 the potential particle model. In the axisymmetric 943 configuration, the experimental AOR is approx-944 imately 4% higher than the one obtained with 945 the clump model and 16% higher than the one 946 obtained for the potential particle model. Also, 947 one should notice that in the plane strain con-948 figuration the JGS method measures an AOR 949 higher than the method presented in this paper, 950 and lower in the axisymmetric configuration (see 951 Figure 12). This changes the conclusion on the 952 influence of the configuration: with our measure-953 ment method both configurations gives the same 954 AOR (difference of approximately 1% with the 955 clump model), while the JGS method gives a 956 difference of approximately 11%. The differences in AOR observed in Table 8 959 between the clump and potential particle (PP) 9 and Table 10.

Figure 20 characterizes the initial and final states of these PP-CLP simulations. First and foremost, it is to notice that the AOR is approximately 14% lower with the potential particle model. This difference can be considered as significant and is even greater with respect to experiments even though the physical particles show a fairly high convexity of 0.954. In line with additional possibilities of interlocking for non-convex particles, while convex particles fall more easily from the heap, the number of lost particles is approximately 23% lower with the clump model.

One can also note that the final void ratio is approximately the same with both models and 11.

1035

The initial void ratio e 0 is controlled by momen- shows α against e 0 and reveals some decreasing tendency of α for e 0 < 0.55 in both configurations that would be more consistent with an interpretation of the AOR in terms of a porosity-dependent maximum friction angle. However, the significant dispersion of the results prompts the need for further investigations in the following.

In addition, it is to note that a critical state interpretation of the repose would impose a correlation between critical state values of porosity (or void ratio) and mean pressure p in the form Fig. 22: Heaps' states after the collapse compared to the critical state line as determined in [START_REF] Duverger | Rattlers' involvement for possibly looser critical states under higher mean stress[END_REF] from a large set of triaxial tests with different initial states in terms of void ratio and/or confining pressure P c heaps is lower than the mean pressure in plane strain heaps, and the dispersion in final void ratio is smaller. Most importantly, the heap states are clearly not consistent with the blue-colored (rattlers excluded) CSL serving as reference, which already suffices to exclude the assumption that a heap of particles under gravity is at critical state.

For completeness, the angle of repose of these CLPX series is still furthermore directly compared in Figure 23 with the critical state friction angle ϕ crit and the (porosity;mean stress)-dependent peak friction angle ϕ peak of the material, determined on the triaxial simulations from [START_REF] Duverger | Rattlers' involvement for possibly looser critical states under higher mean stress[END_REF]. The AOR is therein shown to be significantly different (higher from approx. 10 degrees) than ϕ crit . It actually lies in the observed interval for ϕ peak , even though both are observed to be essentially

ϕ mob = arctan   σ glob 1 -σ glob 3 2 σ glob 3 σ glob 1   (36)
, no obvious correlation is actually found in 1113 1148 • S is the set containing all particles,

N ref part = N part V box 1 3 L η (38) with V box ∈ {V cyl , V par } and L η ∈ {R cyl , L par }.
• C ext is the set containing all the contacts between particles and boundaries,

• the upper-script • p specifies that the quantity is taken for a particle p,

• the upper-script • c specifies that the quantity is taken for a contact c,

• the sub-script • ,xi denotes the derivative with respect to x i ,

• the total volume of the heap is noted V and can be determined from the volumes of all tetrahedrons V tet given by the Delaunay triangulation on particles centers (see previous section 3.3.1): 1254 V = p∈S V tet , 1255

• the number of underline denotes the order of a 1256 tensor (• for vectors and • for matrices),

1257

• the Kronecker symbol δ ij and Einstein's nota-1258 tion will be used,

1259

• classical sign convention for stress is adopted, 1260

where the traction vector t = σ.n applies onto 1261

the system for an outwards normal n.

1262

The global stress tensor σ glob can be expressed 1263 according to the local stress tensor σ:

1264 σ glob = 1 V V σdV (B4)
One can compute σ using the divergence of the 1265 third order tensor σ⊗x (with x the position of any 1266 point in V with respect to a given, even though 1267 arbitrary, origin):

1268 (σ ik x j ) ,k = σ ik,k x j + σ ik x j,k (B5)
Since the measurement is made when the heap 1269 is under equilibrium, the following equation holds, 1270 denoting g the gravitational acceleration and ρ the 1271 particle density:

σ ij,j = -ρg i (B6)
Moreover, x j,k = δ jk , thus:

σ ik x j,k = σ ij (B7)

  experiments Measurements of AOR data have been recently proposed by the JGS as part of a round robin test organised within the activities of Technical Committee 105 (TC105: Geo-Mechanics from Micro to Macro) of the International Society for Soil

Fig. 1 :Figure 2 )Figure 1 )

 121 Fig. 1: Physical particle made of 4 clumped spheres.

  129

131

  This technique leads to an increased total number132 of discrete elements in a simulation, compared to 133 the number of physical particles, however it ben-134 efits from the low computational cost of collision 135 detection among spheres. To define the inertial 136

Fig. 2 :

 2 Fig. 2: Initial (left) and final (right) states of the heap in the axisymmetric configuration.

Fig. 3 :

 3 Fig. 3: Initial (left) and final (right) state of the heap in the plane strain configuration.

171

  additional consideration of using potential parti-172 cles will illustrate the mechanical implications of 173 neglecting the concavity of the physical particles to the AOR.

191

  plane in local particle coordinates, d i is the dis-192 tance of the plane to the local origin and ⟨ ⟩ are 193 Macaulay brackets, i.e., ⟨x⟩ = x for x > 0; ⟨x⟩ = 0 194 for x ≤ 0. 195 This potential function takes zero values (f = 196 0) on the particle surface, negative values (f < 0) 197 inside the particle and positive values (f > 0) out-198 side. In this sense, some similarity can be found 199 with the Level-Set Discrete Element Method (LS-200 DEM) [21, 22, 23] where the potential is the actual 201 distance function, unlike here. The contact point 202 between two potential particles is found as the 203 optimal point of a Second Order Conic optimisation Problem (SOCP) describing the contact detection problem, representing a point nearest to both the particles, based on their potential functions.

236

  To match the local surface curvature of the 237 physical particle, a radius r = r s was chosen in 238 Equation 1 to control the roundness of the edges 239 and corners of the potential particle consistently 240 with the r s radius of each individual sphere in the 241 physical particle. The radius of the shadow par-242 ticle was assigned to R = √ 2r s , to capture the 243 curvature of faces of the given particle shape. The 244 remaining parameter needed to be calibrated in 245 order to match the given particle shape was the 246 parameter k, which controls the curvature of the 247 faces. A value of k = 0.65 led to a good match 248 with the target geometry, i.e. it achieves an ade-249 quate representation of both the overall form of 250 the real particle and features such as its main 251

Fig. 4 :

 4 Fig. 4: Clumped tetrahedral particle (left); fitted potential particle (middle); overlap of the two (right).

Fig. 5 :

 5 Fig. 5: Contact model with visco-elasticity and friction.

Fig. 6 :

 6 Fig. 6: Contact behaviour for different impact velocities in the two DEM approaches.

398

  The actual AOR simulation starts from this 399 initial state by displacing the moving parts of the 400 container in a manner equivalent to the exper-401 iments. Particles leaving the container from its 402 periphery are counted as so-called lost particles 403 and erased for computational efficiency. The sim-404 ulation continues until the sample finds a new 405 equilibrium in the form of a static heap. It is then 406 possible to measure the angle between its exte-407 rior surface and the horizontal plane following the 408 procedures discussed below.

Fig. 8 :

 8 Fig.8: Computation speed statistics during 30 simulations with each DEM approach (either potential particles or clump, see Section 4.3 for details). Dots represent the mean speed value with the surrounding filled area corresponding to its standard deviation.

Figure 9 Fig. 9 :

 99 Figure9shows a typical result after detect-

472 heap configurations except for the orientation 473 of

 473 the horizontal axis. The coordinates ( x, y, z) 474 will thus denote respectively (-r, z, θ) in the 475 axisymmetric case or (x, z, y) in the plane strain 476 case. The width of the container x box for instance 477 stands for R cyl in the axisymmetric case and L par 478 in the plane strain case.

  Figure 11 shows measurements performed for several x max in both configurations. The error bars represent the error computed with Equation 11. One can see that the AOR increases with x max , except for very high values of x max where the part of the outer surface considered is very small compared to its size. This may be caused by the ridge on the bottom of the open container that maintains some particles, affecting the geometry of the outer surface.

Fig. 10 :

 10 Fig. 10: Outer surface regressions for an intermediate x min and the associated measurement ((a), (c)) in both configurations.

Fig. 11 :

 11 Fig. 11: Average slope as measured for different values of x max .

  572 12 shows the AOR measured using CLP1 heaps 573 and one can see that the variation in the mea-574 surement is lower than 3%. Even though such a 575 repeatability error is low, it will be systematically 576 given for all series of simulations in this paper as 577 error bars on the AOR charts.

Fig. 12 :

 12 Fig. 12: Distributions of measured angles of repose when investigating repeatability in the CLP1 configurations of Table4.

3. 3

 3 Measuring the void ratio for any geometry of assemblyWith respect to the objective of discussing the possible constitutive nature of the AOR determined as per Section 3.2, it is interesting to characterize the state of the heap in terms of density or void ratio e, as a fundamental parameter of granular materials. This density characterization is not straightforward because of the irregular geometry of the heap along its free surface, and possible bias caused by an excess of void near the walls As such, two methods are proposed below to compute the void ratio inside a granular assembly with a complex geometry, while avoiding the boundary effects: a so-called "tetrahedra method" and a "sub-volume method". Both methods provide local values for e and rely on a Monte Carlo procedure to compute volume proportions, combined with (straightforward, here) tests to determine whether a random point in space is inside a physical particle. The following differences still exist, though:

611 3 . 3 . 1

 331 The tetrahedra method 612 In order to compute a void ratio on a heap with 613 a random geometry, the tetrahedra method starts 614 by a triangulation of the heap. The Monte Carlo 615 method is then used to determine the proportion 616 of particles inside each tetrahedron resulting from 617 the triangulation, leading to an expression for the 618 is done using Delaunay's trian-623 gulation on the centers of all particles, although it 624 could be done using another set of relevant points 625 (e.g. the center of all spheres for clump simula-626 tions). Also, one should keep in mind that when 627 triangulating using the center of the particles a 628 small part of the sample is ignored: all particles on 629 the outer surfaces are cut by the boundary tetra-630 hedra. This should effectively remove the excess of 631 void near the walls of any sample. The set of all tetrahedra will be denoted {tet}.3.3.1.2 Detecting which particle may be partially inside each tetrahedronAll particles bounding boxes are tested to determine if they overlap a tetrahedron bounding box. If so, the particle is further checked for intersected volume with the Monte Carlo method, forming a set of particles that is denoted {p} cut . This step is not mandatory but it drastically reduces the computation time.3.3.1.3 Computing the total volume of particle inside each tetrahedronIn this final step, N mc points {x i , i ∈ [1, N mc ]∩ N} are uniformly drawn inside the tetrahedron, following[START_REF] Rocchini | Generating random points in a tetrahedron[END_REF]. Each point is tested to determine if it is located inside any of the particles potentially cut {p} cut . Denoting V tet the volume of a tetrahedron tet computed using its vertices' coordinates; χ p (x) the Boolean test function equal to Investigation techniques and physical aspects of the angle of repose of granular matter 1 if the point x is inside the particle p, 0 otherwise; 656 and H(n) the Heaviside function, the Monte-Carlo 657 method gives the total volume of particle inside 658 the tetrahedron as follows:

662 3 . 3 . 2 3 . 3 . 2 . 1

 3323321 figurations. The sub-volume and the total volume of the sample will be denoted V sub and V , respectively. At the final state, the geometry of the sample is assumed to be a half parallelepiped (respectively a cone) for the plane strain (respectively axisymmetric) configuration. The sub-volume is defined using a parameter C that pilots the homothetic transformation. The coordinates of the subvolume axis aligned bounding box are denoted x min , y min , z min and x max , y max , z max and depend on the coordinates of the sample axis aligned bounding box: X min , Y min , Z min and X max , Y max , Z max .

  ) with ∀(s, S) ∈ {(x, X), (y, Y ), (z, Z)} and ∀C ∈ ]0.5, 1].

9 Fig. 13 :

 913 Fig. 13: Illustration of sub-volume for several C values in the plane strain configuration.

) 3 . 3 . 2 . 3

 3323 Counting the volume of particles partially inside the sub-volume If 0 < N p in < 8, the particle may be cut by the faces of the sub-volume. The proportion of the particle volume inside the sub-volume is again determined using the Monte Carlo method: N mc points, {x i , i ∈ 1; N mc }, are uniformly drawn inside the particle bounding box and tested to

Fig. 14 :

 14 Fig. 14: Illustration of sub-volume for several C values in the axisymmetric case.

3. 3 . 3

 33 Examples of void ratio measurements 3.3.3.1 Local void ratio The tetrahedra method makes it possible to establish directly a local representation of the void ratio, as illustrated in Figure 15 for one plane strain final heap. One can notice that the geometry of the final heap is accurately captured by the triangulation, giving a rounded half parallelepipedic boundary surface. The density range is quite wide: some tetrahedra located on the outer surface, where the particles moved, contain approximately 1000 times more voids than other tetrahedra located where the particles almost didn't move. Note that this figure represents the void ratio directly interpolated from the centroids of each tetrahedron and thus should be interpreted carefully.

Fig. 15 :

 15 Fig. 15: Local void ratio in a plane strain final heap as measured with the tetrahedra method.

Fig. 16 :Figure 16

 1616 Fig. 16: Parallelization speed-up for the tetrahedra method.

  were performed for 802 the clump model on all 30 initial samples of 803 CLP1 plane-strain series of simulations, discussed 804 in more details in Section 4.3. Because of the 805 simple parallelepipedic geometry of these granular assemblies, a reference void ratio can be easily computed using the sample bounding box:

Figure 17 (

 17 Figure 17 (a) illustrates the comparison between the mean values and standard deviation over the 30 samples of e REF together with e T ET for the tetrahedra method and e SU B for the subvolume method. The latter has been computed for 3 values of N mc and 40 values of C. For the lowest values of C, the measured e SU B void ratio varies a lot among the 30 simulations and in function of N mc . Between C ≈ 0.7 and C ≈ 0.9, e SU B is constant and its standard deviation gets lower, being furthermore little dependent on N mc . For C > 0.9, its mean value and standard deviation finally start to increase as expected due to the rigid boundaries constraining the granular assembly and favoring voids to form near the outer surfaces. Finally, for C = 1, the sub-volume method gives by definition the exact same values for void ratio than when using the global bounding box: e SU B = e REF .

8494. 1

 1 Parametric study 4.1.1 (Non-)Sensitivity to the tangential stiffness A first series of simulations investigates the role of tangential stiffness when using the potential particles model and two different values of K s : 240 N/m and 444 N/m (see sets B and C of Table

(Fig. 17 :

 17 Fig.17: Void ratio measurements with both methods for 30 parallelepipedic initial samples differing in individual locations of particles. On figure (b), the y axis is broken at two places: first between 40 ns and 50 ms, and second between 250 s and 5200 s. The three parts of the y axis do not have the same scale.

879cost.

  Other series of simulations with different ρ 880 are thus proposed to check whether a variation 881 from the experimental reference ρ = 1, 111 kg/m 3 882 would affect the AOR results. 883 Using potential particles, two values for the 884 particle density are considered in the framework 885 of the PP1 series (sets A and B of Table 5): 886 the experimental one, ρ = 1, 111 kg/m 3 , and 887 ρ = 943 kg/m 3 that would confer the potential 888 particle the same mass as the physical parti-889 cle in spite of the volume differences discussed 890 in the above Section 2.2.3. Using clumps in a 891 CLP2 series, four to seven different values for 892 ρ ∈ [100 kg/m 3 ; 10, 000 kg/m 3 ] are considered, 893 with 10 different initial samples in each case.

Figure 19 shows(Fig. 18 :

 1918 Figure19shows the resulting angles of repose in the CLP2 and PP1 series, together with time costs of CLP1 series measured for a sequential execution on the same machine presented in previous Section 3.1. The dots correspond to the mean measurement over all samples for a given particle density and the error bars represent the standard deviation. One can see that all error bars share a common zone for a given shape description. As such, it is herein concluded, consistent to[START_REF] Zhou | Numerical investigation of the angle of repose of monosized spheres[END_REF], that particle density does not impact the AOR. During DEM simulations, one can thus adopt, when necessary, an artificial ρ = 10, 000 kg/m 3 , multiplying the critical time step by a factor of 10000 1111 ≈ 3

Fig. 19 :

 19 Fig. 19: Particle density influence on the angle of repose and on the time cost -CLP2 and PP1 series

  960 models, with a higher discrepancy for PP towards 961 experiments, certainly arise from the convex sim-962 plification of potential particles, with respect to 963 the concavities of the physical particles which 964 allow them to interlock better. In order to gain 965 more insights into the influence of particle con-966 cavity, a rigorous comparison between the two 967 numerical models is led in this subsection, adopt-968 ing the same parameters for both models (except 969 for the time step for computational efficiency) 970 and determining the AOR using the more reli-971 able method presented in Section 3.2. This series 972 is called PP-CLP, with all parameters being listed in Table

Fig. 20 :

 20 Fig. 20: Characterization of the PP-CLP series on particle concavity in terms of initial and final states

1036

  tarily altering the inter-particle friction angle dur-1037 ing the generation of the sample, φ gen , whereby 1038 lower φ gen -values lead to denser packings, as plot-1039 ted in Figure 21 (a) where the error bars represent 1040 the standard deviation of e 0 on all 10 simulations performed at the same φ gen . Figure 21 (b)

Fig. 21 :

 21 Fig. 21: Influence of the initial void ratio e 0 -CLP4 series

  ) P c = 0.25 kPa P c = 2 kPa P c = 10 kPa P c = 40 kPa P c = 100 kPa P c = 200 kPa P c = 400 kPa Angle of repose (Axisymmetric) Angle of repose (Plane strain)

σ glob 1 ≥ σ glob 3 in

 13 different.From the observations that the repose states are not consistent with the shear strength properties of the granular material, neither the critical one nor the maximum one, the AOR is concluded to bear no constitutive nature. Interpreting the repose stress state σ glob ij with its extreme principal stresses terms of a mobilized friction angle ϕ mob , given in Equation 36 using the soil mechanics sign convention:

Figure 24 ,respect to e 0 Fig. 23 :Fig. 24 :

 2402324 Figure 24, no matter the shape model, between 1114 the mechanics of the heap, ϕ mob , and its geom-1115 etry, α, which would have been mandatory for a 1116 constitutive interpretation.

Figure 25 Fig. 25 :

 2525 Figure 25 shows the AOR values obtained in this CLP5 series, with error bars from the standard deviation computed on the ten simulations performed for each value of N part . An exponential model is proposed to fit the data and provide an extrapolated value of α for an infinite number of particles:

1149

  For bigger systems with a higher N part a 1150 clear difference appears, with a negligible stan-1151 dard deviation. According to the exponential 1152 decay model, an asymptotic difference difference 1153 between the two configurations would be 1.19 • 1154 which corresponds to approximately 3.59% of the 1155 measurement in the plane strain configuration. 1156 The present dependency to N part constitutes 1157 a last argument against the constitutive nature 1158 of the AOR, making fortuitous the near-similarity 1159 of α initially observed between the two configura-1160 tions, in connection with the particular values of 1161 N part used in the reference experiments.

-β ω π -arctan 2βω ω 2 -β 2 for β < ω0 √ 2 exp -β ω arctan 2βω ω 2 -β 2 for β ∈ ω0 √ 2 ;= ω 0 2 -β 2

 22222222 In the framework of a round-robin activity providing an experimental reference, the angle of repose of a granular material has been studied with DEM, adopting two distinct characterizations for particle shape: concave clumps of spheres and convex-simplified potential particles, with a quantification of the morphological differences between the two (and the experimental reference), in terms e.g. of convexity and sphericity.A methodological discussion has then been first proposed for generic angle of repose studies, designing systematic measurement procedures of the slope angle and of the void ratio of the heap, the latter being a possible factor of influence onto the former, as a fundamental property of granular matter.Physically, a thorough analysis provided a number of consistent observations that the AOR measured on the considered devices does not bear a constitutive nature but is instead processdependent. In the comparison with experiments, while adopting a simpler measurement method of the AOR due to experimental limitations, the clump approach successfully predicted the AOR within a 8% tolerance. On the other hand, the potential particles underestimated to a greater extent the AOR, as expected due to their artificial convexity. Even though the material particles had a fairly high convexity value (C = 0.954), neglecting their local concavities brought down the AOR from 35.95 ± 0.88 • to 31.26 ± 0.95 • . It is interesting to note that both the clump and the potential particle shape descriptions share the same dimensions in terms of a minimal bounding box and thus the same flatness and elongation values, prompting the need for a systematic investigation of other particle-scale shape indices that would possibly affect the AOR. Appendix A YADE implementations of visco-elasticity With reference to the normal contact law in Equation 2 and for the clump approach (which use here the Ip2 ViscElMat ViscElMat ViscElPhys and Law2 ScGeom ViscElPhys Basic YADE classes), c n is computed from a given normal restitution coefficient e n according to the mass-1210 dependent expression of [30], recalled in Equation 1211 1213 and Ω = β 2ω 0 2 , with m being the harmonic-1214 average particle mass and K n the normal contact 1215 stiffness. Equation A1 is solved inside YADE using 1216 a small number of Newton-Raphson iterations to 1217 make the inverse calculation of c n based on the 1218 desired value of e n . A straightforward calculation 1219 of c n for the linear contact model has been pro-1220vided in Equation (B4) of[START_REF] Thornton | An investigation of the 1495 comparative behaviour of alternative con-1496 tact force models during inelastic collisions[END_REF] via curve-fitting of 1221 the exact solution of[START_REF] Schwager | Coefficient of restitution and linear-dashpot model revisited[END_REF], which is not however 1222 employed here.1223For the potential particles approach (through, 1224 e.g., Ip2 FrictMat FrictMat KnKsPhys and 1225 Law2 SCG KnKsPhys KnKsLaw classes), a viscous 1226 damping parameter β n serves as input for deriv-1227 ing c n , consistently with a desired e n and[START_REF] Antypov | On an analytical 1499 solution for the damped hertzian spring[END_REF] via 1228 viscous damping coefficient is calculated as in Equation A3: c n = 2β m • K n evaluate the stress state of the heap, one has to compute the stress tensor from the contact forces of all the contacts. Moreover, gravity being present in the simulations and at the origin of the movement, gravitational forces should thus be accounted for. In this subsection:

Fig. C1 :

 C1 Fig. C1: Wall velocity influence on the angle of repose and on the time cost -CLP3

Table 1 :

 1 Default configuration of AOR simulations

	Configuration Initial number of particles	Side wall velocity	Container height	Container width
	Axisymmetric	N part = 2, 468

Table 2 :

 2 Coefficients defining the planes making the faces of the tetrahedral potential particle as described in Equation1.

		Plane coefficient	Plane 1	Plane 2	Plane 3	Plane 4
		a b	0 0	√	2/3 2/3	2	0 √ 2/3	-2/3 √ 2/3
		c	-1	1/3		1/3	1/3
		d (cm)	0.063299	0.063299	0.063299	0.063299
	229	the planes used to assemble the potential parti-
		cle were chosen as the faces of the tetrahedron

230

connecting the centers of the spheres making the 231 physical particle (see Table

2

). This approach can 232 be generalised to approximate any convex shape, 233 given a tessellation of its surface, or a multi-sphere 234 representation of a particle made of spheres with 235 equal radii.

Table 3 :

 3 Shape parameters of the physical particle in comparison with various DEM approaches

	Shape characteristics	(1) Physical particle or present clump approach		(2) Potential Particle		(2)-(1) (1)		(3) Clump approach with non-uniform density	(3)-(1) (1)
	Volume (m 3 ) Surface area (m 2 )	3.3304 × 10 -7 2.491 × 10 -4 2.584 0  0 		3.9248 10 -7 2.632 × 10 -4 3.286 0 0 	17.85% 5.66%		4.9965 × 10 -7 2.491 × 10 -4 3.123 0 0 	50.03% 0
	Inertia tensor/ρ (m 5 )		0 2.584 0 0 0 2.584  × 10 -12		0 3.286 0 0 0 3.286	 × 10 -12	27.17%		0 3.123 0 0 0 3.123  × 10 -12	20.86%
	Convexity		0.954		1		4.82%		0.954	0
	True sphericity		0.9328		0.9849		5.59%		0.9328	0
	density i.e. inner overlaps would lead to an error
	of 50.03 % for the volume and 20.86 % for the
	eigenvalues of the principal inertia tensor.

  Unless specified otherwise, the DEM clump simulations are defined accordingly, assign-378 ing to all contacts random values of φ and c n 379 which respect the same distributions.

	385	Section 2.1, so as to form a extremely loose assem-
	386	bly of non-overlapping particles. The assembly
	387	is then deposited under its own weight until it							380
	389	becomes stable, and is saved to be subsequently used under different conditions. Different samples can be obtained starting from different initial par-	0.025 0.050 0.075 0.100 0.125 Relative frequency					resin against resin resin against acrylic
			0.000	20	25	30	( )	35	40	45
			(a) Contact friction angles for particle/particle
			(φ p/p ) or particle/wall (φ p/w ) contacts
			0.1 0.2 0.3 0.4 Relative frequency 0.5 0.6		resin against resin resin against acrylic
			0.0	en 0.70 0.72 0.74 0.76 0.78 0.80 0.82 0.84
						(b) Restitution coefficient
			0.05 0.10 0.15 0.20 Relative frequency
			0.00	Kn (N/m) 35000 40000 45000 50000 55000 60000 65000 70000 75000
						(c) Normal stiffness
			Fig. 7: Distributions of contact properties as
			experimentally measured by the JGS (adapted
			from [4])
			2.4 DEM simulation workflow	381
			Building DEM samples starts with generating 382
			randomly in space particles inside the cuboidal 383
			or cylindrical containers mentioned in the above 384

390 ticle arrangements. For this first step that has no 391 experimental counterpart, an extra, non-physical, 392 damping source is added in the local, non-viscous, 393 form (Cundall's damping) to speed up the genera-394 tion. For the rest of the simulations, contact-scale 395 viscous damper and friction solely ensure the sta-396 bilisation of the simulations and no other source 397 of damping is used.

  by the discharge and still form a horizontal Investigation techniques and physical aspects of the angle of repose of granular matter surface, especially in the axisymmetric configura-

	504
	505	tion. Indeed, finite particle-size effects necessarily
	506	exist and affect the transition from one side of the
	507	slope to another.
	508	Excluding from the bottom of the heap the
	509	few particles that are stuck by the ridge, and only
		those, is obtained choosing: x min = 0.32d clump .
	527

510

An appropriate value for x max is sought by mea-511 suring α for several x max . The best x max is the 512 smallest for which the measurement does not 513 change. The error on the measurement is also a 514 criterion to choose the best x max . This method 515 should be specially relevant in the axisymmetric 516 case since the outer surface is curved, but it should 517 work on the plane strain heap as well.

518 3.2.3 Error on the measurement 519 For a given heap, the dispersion of positions data 520 induces some error on the linear regression and the 521 measurement of α. As an alternative to the corre-522 lation coefficient R 2 , this error can be quantified 523 from a standard deviation on the slope a 1 of the 524 fitting line, StD(a 1 ). If N is the number of points 525 and ( x i , y i ) are the coordinates of the i th point, 526 one has:

Table 4 :

 4 Parameters for heap simulations investigating repeatability (CLP1 series, 60 simulations in total).

		Configuration N part	K n (N.m -1 )	K s /K n	ρ (kg.m -3 )	∆t (s)	φ	e 0	Number of samples
		Plane strain Axisymmetric 2, 468 2, 150	58, 250	0.37	1, 111	7.86 × 10 -5	see Figure 7 0.744 ± 0.028 0.622 ± 0.012	30
	609	which can be difficult if the latter adopts a pecu-
	610	liar geometry. However, it is substantially faster
		than the tetrahedra method.

Table 5 :

 5 Used parameters for the parametric study with potential particles -series PP1.Set id Configuration K n (N.m -1 ) K s /K n ρ (kg.m -3 ) ∆t (s) φ p/p ( • ) φ p/w ( • )

	β n	Number of samples

present explicit DEM scheme and the total time

Table 6 :

 6 Parameters of heap simulations investigating ρ influence (CLP2 series, 110 simulations in total).

	Configuration N part	ρ	e 0	Number of samples
	Plane-strain	2, 150

Table 7 :

 7 Material properties used in the potential particles models.K n (N.m -1 ) K s /K n ρ (kg.m -3 ) ∆t (s) φ p/p ( • ) φ p/w ( • )

	β n

Table 8 :

 8 Angle of repose as per the JGS measurement method.

			Number of samples	Average ( • )	Standard Deviation	( • ) Minimum ( • ) Maximum ( • )
		Plane strain configuration
		Experiments [4]	400	41.4	1.28	38.3	46.3
		Clumps	100	38.1	1.14	35.0	41.3
		Potential particles	18	34.8	1.61	32.5	38.0
		Axisymmetric configuration
		Experiments [4]	50	35.3	0.9	33.3	37.3
		Clumps	100	33.9	0.8	32.0	36.1
		Potential particles	19	29.7	0.78	28.5	31.2
	958	4.3 Role of particle concavity

Table 9 :

 9 Contact parameters of the PP-CLP series focusing on particle concavity Clump 1.2 kN.m -1 0.24 kN.m -1 0.8 35.5 • 27.2 • Potential particle

	Model	K n	K s	e n	φ p/p φ p/w

Table 10 :

 10 Other simulation parameters of the PP-CLP series on particle concavity

	Model	Configuration	ρ	∆t	Number of samples
	Clump Potential particle	Plane strain	1111 kg.m -3 ≈ 78.5 µs ≈ 85.2 µs	30
	that the final average number of contact points	they are free to move, they do not fill the voids 1010
	per particleis approximately 25% lower with the	better than convex particles. On the other hand, 1011
	potential particle model since two convex parti-	the difference on the average number of contact 1012
	cles can form only one contact point, unlike the	points is approximately the same at the initial and 1013
	concave clump.			final states.		1014
	It is worth noticing that the experimental mea-
	surement method gives a gap between the axisym-
	metric and plane strain configurations' angle of
				4.4 (Non-)Constitutive nature of	1015
	repose. The curvature of the axisymmetric heap
	explains the lower measurement obtained with the	the angle of repose	1016
	JGS method.			The compatibility of AOR values measured (with 1017
	Looking at initial stages, one observes that the	the present method, see Figures 12 and 19) in both 1018
	initial void ratio is approximately 5% lower with	axisymmetric and plane strain configurations of 1019
	the clump model. The difference between the void	the reference device could confirm a constitutive 1020
	ratio at the initial state and final state is interest-	nature of the AOR inferred in, e.g., [5, 6, 7], sug-1021
	ing: at the initial state the sample is constrained	gesting to compare the latter with shear strength 1022
	by four side walls and one bottom wall, while at	properties of the granular material. Generally 1023
	the final state one of the four side of the sample is	speaking in solid-like granular mechanics, these 1024
	free. This suggests that concave particles are more	shear strength properties may refer either to a crit-1025
	likely to fill the voids when there are surrounded	ical state or a state of maximum stress ratio, the 1026
	by walls (parallel to the gravity axis), but when	two being possibly different depending on initial 1027
				porosity.

Table 11 :

 11 Parameters used when investigating a possible influence of e 0 (CLP4 series, 200 simulations in total)

	Configuration N part	ρ			e 0	Number of samples for each e 0
	Plane strain	2, 150	10, 000 kg/m 3	0.414 ± 0.010 0.459 ± 0.009 0.480 ± 0.008 0.495 ± 0.006 0.504 ± 0.008 0.528 ± 0.006 0.536 ± 0.008 0.553 ± 0.010 0.567 ± 0.012 0.574 ± 0.007	10
	Axisymmetric 2, 468			0.419 ± 0.006 0.486 ± 0.006 0.525 ± 0.009 0.553 ± 0.006 0.577 ± 0.011 0.607 ± 0.018 0.626 ± 0.014 0.632 ± 0.010 0.655 ± 0.012 0.671 ± 0.016
	Configuration	N part	ρ		e 0	Number of samples
	Plane strain	2, 150 2, 468 5, 000 10, 000	10, 000 kg/m 3	0.603 ± 0.010 0.604 ± 0.009 0.596 ± 0.005 0.596 ± 0.003 0.594 ± 0.002 0.591 ± 0.003 0.588 ± 0.001	10
	Axisymmetric	20, 000 30, 000 50, 000			0.705 ± 0.011 0.693 ± 0.010 0.729 ± 0.017 0.741 ± 0.006 0.758 ± 0.008 0.746 ± 0.016 0.736 ± 0.004

Table 12 :

 12 Parameters of heap simulations investigating the influence of N part (CLP5 series, 140 simulations in total)

Acknowledgements

We thank the Japanese Geotechnical Society for organizing the DEM round robin test that motivated this study. The efforts of Dr Shuji Moriguchi (Tohoku University) are in particular gratefully acknowledged. We also acknowledge the support from the French "Sud" region to the recent LS-ENROC project with an acquisition of a server machine that could accomodate the 980 simulations presented herein.

one has to compute the x and y coordinates of the 698 center, x Ω and y Ω respectively, and the maximum 699 radius r c of the cone (Figure 14). The homothetic 700 sub-volume is then:

s Ω = S min + S max 2 (21)

with ∀(s, S) ∈ {(x, X), (y, Y )} and ∀C ∈]0.

The total volume of particles partially inside 733 the sub-volume is:

The total volume of particle inside the sub-Equation B5 one gets:

Equation B4 then gives:

Particles having an homogeneous density, one 1276 furthermore has, with m p and x p the mass and 1277 center of p :

The part of σ due to gravity (σ G ij ) can thus be 1279 written:

As for the part due to contacts (σ C ij ), Green-

Ostrogradski theorem gives:

Considering the traction vector t = σ.n, one has:

The traction vector is not nil only on contact points. Since the system is closed, contact forces between particles cancel each other leaving only forces comming from outside of V . As a consequence, one can only consider the contact forces between particles and walls. For these contacts f c denotes the contact force exerted by the wall on the particle and x c the contact point. One has:

Finally, the global stress tensor for a stable heap of particles made of homogeneous particles and subjected to gravity is:

The mean stress can then be computed as p = Tr(σ glob ) 3 .

wall velocity

In both configurations a wall holding the particles moves in order to let them fall.