{'fraction_non_alphanumeric': 0.024718921936049756, 'fraction_numerical': 0.008532014990830077, 'mean_word_length': 5.091306459446333, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discrete simulation of railway ballast shear test: spherical and polyhedral grain shapes

Dynamical behaviour of railway ballast is important in design and maintenance of railroads. Various modelling approaches are available. The least phenomenological but also the most computational demanding approach is direct representation of every ballast grain and its interaction with other grains in the model. The computational complexity is partially reduced by assuming non-deformable bodies with simplified shapes. Two different shapes of the bodies are considered: spherical and polyhedral. Spherical shapes are advantageous because of their great computational simplicity; however, missing shape information must be compensated by adding phenomenological rolling resistance. These two model variants are used to simulate different shear test. Results of the models are compared to each other and also to the experimental data from the literature.


Introduction
An assumption of material homogeneity is generally used to investigate many complex problems in the railway engineering, for example, railway sleepers [1], wheel/rail sliding contact [2] or fatigue of railway axles [3]. Homogeneous continuous formulation in connection with complex phenomenological constitutive relations can be suitable tool for mathematical modelling of the observed phenomena at the macroscopic scale. But with such an approach, it is difficult to take into account phenomena such as dilatancy, pressure sensitivity, compaction, crushing etc. The discontinuous heterogeneous nature plays an important role in the behaviour of granular materials.
The Discrete Element Method (DEM) is an approach that takes an advantage of the granular nature of the material. Nowadays, DEM is one of the most useful and efficient tools to describe the behaviour of discontinuous materials. In the Discrete Element Method, originally proposed by Cundall and Strack [4], every grain is represented by an ideally rigid and independent body which can interact with the other bodies through forces and moments at their contacts. The rigidity of the bodies is motivated by assumption that macroscopic material behaviour is mostly dictated by rearranging of the grains, not by the actual grain deformations. The DEM typically uses explicit integration of the motion equations. The time is discretized time into small time steps ∆t of size limited by stiffness and size of the grains. In every time step, it is necessary to detect all contacts between bodies and evaluate contact forces. The computational cost of this operation significantly limits the applicability of the DEM to some practical, large scale problems, where millions of particles are typically involved. The difficulty of contact detection and force evaluation depends on shape of the bodies that represent the real grains of the railway ballast. We consider two common possibilities of shape simplification. one is to use convex polyhedral shape which closely corresponds to real shape of ballast grains [5,6,7]. The simple one is to use spherical bodies [8,9]. Other options, not considered here, involve e.g. clumps of spheres [10], ellipsoids [11] or superquadrics [12]. The simple spherical shape requires enhancement or constitutive law to provide sufficient mutual locking of bodies. The enhancement consists in additional phenomenological rolling resistance. The rolling resistance mimics resistance of two grains of complex shape in contact rolling relatively to each other [13 14, 15].
Motivated by difficulties in adjusting the rolling model parameters [16,17], we are presenting this study. It aims at comparison of performance of polyhedral and spherical shapes in simulation of shear test. The shearing is performed after compaction that is realized via application of large pre-stressing.

Derivation of equivalent constitutive relation for spheres and polyhedras
The polyhedral and spherical shape use completely different formulation of constitutive relations. For comparison purposes, material behavior for both models must be approximately equal. Simple relation between material parameters can be established under assumption of ideal polyhedral particle of spherical shape. In this section, we analyze in detail contact of spherical particles using two formulations of the constitutive law, the one for spheres and the one used by polyhedrons. Relations between material parameters for polyhedras and spheres are then established by comparison of these two cases. Superscripts (s) and (p) are used to distinguish between the constitutive law formulations.
Two spherical particles in contact can mutually interact in different modes by forces and moments that are sketched in Figure 1. The constitutive formulation for spheres is based on Hertz-Mindlin [18,19] with additional rolling resistance. The polyhedral formulation uses intersecting volume to determine repulsive force [20]. The twist mode is neglected in both cases. All the presented equations are taken from open source platform YADE [21].
For the Hertz-Mindlin contact type, one can express the normal force as R is equivalent radius of spheres is Young's modulus, ν (s) is ration between normal and shear stiffness, R=R 1 +R 2 is radius of both spheres, G is mutual overlapping of particles and n is a normal in direction determined by particle centroids. The equations (1)-(3) combined together gives For polyhedral shapes, we can express the normal force as with effective Young's modulus for polyhedrons The cap of sphere with depth 2 and the volume of the intersection of spheres is then Then the equations (5), (6) and (8) provides Requiring the same force produced by both constitutive formulation gives us After the modification of equation (10), we can compute required ratio of stiffness parameters, k ( ) The same approach can be repeated for tangential direction. The expressions for shear directions [21] shall produce the same force.
Both for spheres and polyhedras, the maximum shear force is limited by Coulomb friction model [21] with additional parameter, φ, called angle of internal friction same for both constitutive formulations.
The both derived relations are dependent on the radius R and mutual overlapping δ. The dependence on overlapping δ is disadvantageous because it varies during the simulations and has to be approximately chosen. Moreover, the polyhedral particle shapes are far from the spherical shapes. Thus, the derived relations cannot ensure equal results of the simulation. However, they still provide reasonable guess of the values for both models.
The spherical constitutive formulation involves also the rolling resistance. The rolling resistance [13,14,15] is an artificial moment arising on the contact of two discrete elements that mimics resistance of two grains of complex shape in contact rolling relatively to each other. It is necessary to use rolling resistance when spherical elements are employed to reproduce behaviour of non-spherical particles such as grains in railway ballast. The model of rolling resistance assumes linear dependence of the rolling moment increment, ∆M R , on the rolling angle increment, ∆ω R , via the stiffness in rolling, k R [Nm] [21].
The stiffness in rolling is calculated based on parameter α R [-] Similarly to the Coulomb shear, the upper limit is imposed on the rolling moment magnitude.
Parameter μ R [-] controls the maximum rolling moment similarly to the tangent of internal friction angle in shear.

Time step
The estimation of time step is based on the contact stiffness, the material density and particle size. For spherical particles the critical time step is typically considered as where R eq is the equivalent radius of sphere, ( ) eff s E is effective Young's modulus and ρ is the density. For polyhedrons the equivalent volume radius R V is set and the critical time step is consequently estimated.
In case of spheres, the equation (18) is rigorously derived. It ignores the shear and rolling stiffness, that are typically much lower compared to the normal stiffness. The estimation of critical time step for polyhedrons is rather empirical and needs further improvements.

Influence of particle shape in the shear test
The comparison between spherical and polyhedral discrete elements is shown on simulation of the shear test. The experiments published in [22] are used for comparison. Two separate steel boxes are filled by a granular material (the spherical/ polyhedral discrete elements in simulations). A steel loading plate that provides vertical surcharge during shearing is placed on the top. The bottom part of the box is horizontally displaced up to 37 mm, while the upper part remains at the original position. The dimensions of a large scale shear box ( Figure 3) were a = 300 mm, b = 300 mm, and c = d = 100 mm. The vertical pressure p, that is kept constant during the loading, is 15 kPa. Before the shear loading starts, the sample is compacted by an incresed vertical pressure. The influence of this compacting pressure is studied, the results are compared with the experiments [22]. Three variants differing in shape of the grains are numerically simulated: (i) spherical shapes of constant radius R, (ii) convex polyhedrons of random geometry of size comparable to spheres with radius R (iii) spheres with radii randomly sampled from a uniform distribution over interval from 0.5R mm to 1.5R mm. The radius R corresponds to the mean radius R 50 from the experiment [22], but the experimental sieve curve was not represented. The maximal radius of laboratory tested ballast was 53 mm.
The material a geometrical properties are shown in Table 1. The elastic parameters for the polyhedrons are estimated based the equations (11) and (13) assuming δ=1.42×10 -4 m. The remaining properties of both particles were identified in previous research or assumed from the literature.  Figure 4; all the curves are averages over three samples computed. The first row of graphs shows dilatancy -evolution of vertical displacement of loading plate v on horizontal displacement of the bottom box h. Both spherical model variants show similar results; there is an initial sharp peak dependent on compaction level followed by monotonic slow increase of volume. Initially, the sample has to rearrange the grains and the loading plate goes up. The later stable growth is not sufficiently explained yet, but it does correspond to the experimental measurement. The polyhedrons are much more sensitive to compaction -omitting compaction leads to large contraction. For compacted samples, one can see the same compaction dependent peak as in the case of spheres, but followed by almost constant level of vertical shift v. Theoretically, after jumping behind the first peak, the compaction effect should diminish. This can be seen in graphs in the third row, where the curves of the bulk density (total mass of the sample divided by total volumes including pores) are getting close to each other.
The results for shear stress show disorder and large variability among the various samples. This is attributed to relatively low number of particles in the boxes. In case of polyhedrons, the measurement of shear stress exhibits lot of jumps and sharp peaks. The reason is probably inadequate length of the time step ∆t. Two polyhedral bodies may come into contact by protruding parallel facets, giving rise to sudden increase in intersecting volume and large repulsive force. The smaller time step may reduce the sudden increase of volume and reduce the jumps in the shear stress. The authors intent to improve the estimation of time step necessary to obtain stable simulation in their future work.

Conclusion
The normal and shear stiffness of polyhedral grains were estimated based the comparison of their volumetric constitutive law and Herz-Mindlin constitutive laws for spheres. The calculation of corresponding parameters depends on mutual overlapping of particles δ, which has to be estimated.
Simulations of shear test using polyhedral and spherical particles were performed. The shearing was realized after application of various vertical compaction pressures: no compaction (15 kPa), compaction 200 kPa and compaction 2000 kPa. The polyhedral model yields the highest sensitivity to compaction, however also spherical grains were significantly influenced -polydisperse assemblies more than monodisperse. The simulations did not consider the breakage of ballast grains. Including this phenomena may lead to improvement of the results.