Discrete element modelling of scaled railway ballast under triaxial conditions

The aim of this study is to demonstrate the use of tetrahedral clumps to model scaled railway ballast using the discrete element method (DEM). In experimental triaxial tests, the peak friction angles for scaled ballast are less sensitive to the confining pressure when compared to full-sized ballast. This is presumed to be due to the size effect on particle strength, whereby smaller particles are statistically stronger and exhibit less abrasion. To investigate this in DEM, the ballast is modelled using clumps with breakable asperities to produce the correct volumetric deformation. The effects of the quantity and properties of these asperities are investigated, and it is shown that the strength affects the macroscopic shear strength at both high and low confining pressures, while the effects of the number of asperities diminishes with increasing confining pressure due to asperity breakage. It is also shown that changing the number of asperities only affects the peak friction angle but not the ultimate friction angle by comparing the angles of repose of samples with different numbers of asperities.


Introduction
The majority of railway tracks in the world are still using ballast in their design because of their relatively low cost compared to concrete slab tracks [1]. In order to design ballasted track and maintain it properly, it is important to fully understand the mechanical behaviour of the ballast material.
B H. Li lihuiqireal@outlook.com 1 Nottingham Centre for Geomechanics, University of Nottingham, Nottingham, UK With regards to most geomaterials, the triaxial test is one of the most useful laboratory tests for investigating the deformation and strength of railway ballast [2][3][4][5][6][7][8][9]. Indraratna et al. [2] performed large scale triaxial tests on latite basalt under various confining pressures and found the angle of internal friction was a function of the confining pressure and particle breakage was more pronounced at higher confining pressure. Aursudkij et al. [6] also performed large scale triaxial tests on limestone particles under confining stresses of 10, 30 and 60 kPa and the sample with lowest confining stress was found to give the largest amount of dilation but the smallest amount of breakage. Additionally, as a powerful numerical tool to understand the micromechanical behaviour of granular material, the discrete element method (DEM) [10] has shown some success in simulating triaxial tests on granular materials [11][12][13][14][15][16][17][18][19]. However, it is a challenge to consider both irregular particle shape and ballast degradation in DEM modelling of railway ballast; e.g. Lobo-Guerrero and Vallejo [20,21] published work on DEM of ballast degradation but only circular particles were used, and the studies where the complex ballast particle shapes were modelled, particle breakage was ignored [22][23][24][25]. Considering that most ballast degradation is not attributable to particle splitting but instead primarily particle abrasion [4,5], Lu and McDowell [14,15] introduced a tetrahedral shaped clump with small breakable asperities ( Fig. 1) to represent a ballast particle, which as far as the authors are aware, was the DEM study of abrasion with an irregular shaped particle. Although this asperity breakage model cannot represent bulk fracture, it conserves mass and presents the fragment movements as a result of the crushing mechanism. Ballast particles are usually angular,uniformly graded [26], tend to be blocky and equi-dimensional-so the tetrahedral shape was chosen to fit these criteria. By using this model they successfully explained firstly the difference in shear strength at low and high confining pressure, and the associated reduction in dilation at high confining pressure observed on full-sized ballast particles.
However, scaled ballast has also been used in some cases to eliminate the requirement for large triaxial apparatus. The triaxial tests on scaled ballast [27] show generally similar stress strain behaviour to the full size ballast experimental results [2,3,6,9], except that the Mohr-Coulomb friction angles corresponding to the peak strengths are less sensitive to the effective confining pressure, as Fig. 2 shows. This feature is not consistent with the simulation results obtained by Lu and McDowell [15]. The purpose of this study is to modify Lu and McDowell [15]'s breakable particle model (Fig. 1) to make it applicable to smaller (scaled) ballast and quantitatively compare the results, in terms of mobilized angle of friction and volumetric strain as a function of axial strain, with experimental tests [27]. The impact of the number of asperities and their strength are investigated and optimal parameters are deduced to reproduce the behaviour of ballast within the range of pressures typical of railway track conditions. The first section of this paper describes the discrete element model with breakable asperities while the second one summarises the triaxial test conditions and the parameters of the models. The third part compares the model and experimental results  [15] of the triaxial tests at confining pressures of 15 and 200 kPa, to represent the extreme cases for in-situ conditions. A real challenge is developing a model that is capable of capturing the peak strength of the ballast over a range of confining pressures, but which does not significantly affect the ultimate angle of shearing resistance as a function of confining pressure -this is addressed in this paper and is considered a novel contribution. It is important to note that this work was carried out alongside the University of Southampton as part of the 'Track to the Future' project. Each research group independently devised a method for modelling the behaviour of scaled ballast with reference to the same laboratory data [27]. The work at Southampton, in a paper submitted for review [28] simultaneously with this one, proposed a new contact damage model as an alternative method for modelling the observed behaviour.

DEM model
The discrete element method (DEM) considers a granular material like ballast as an assembly of objects interacting through a contact law. The contact interaction is generally modelled using a linear stiffness governing the overlap between the contacting objects. This normal and shear contact forces are characterised by normal and shear stiffness coefficients k n and k s respectively. A friction coefficient f based on the Mohr-Coulomb criterion limits the ratio between the shear and normal components of the contact force. Although more complex contact models have been developed, this basic linear elastic model is embedded in Fig. 2 The Mohr-Coulomb peak friction angles versus different confining pressure the code PFC3D [29] and used in the present study. The acceleration of the object is first deduced from the contact forces acting on it using Newton's second law and then its velocity and displacement are deduced by integration over a time increment. The displacements of the objects lead to new overlaps between them from which new contact forces arise forming a full DEM cycle calculation. The motion of the objects is determined using a succession of time increments. In order to avoid non-physical oscillation within the assembly of objects, non-viscous damping reduces their accelerations by 70 % (by default) to allow further dissipation of energy in addition to friction at contacts [29].
The most common object used in DEM is the sphere as it makes contact detection and overlap determination between two objects simple. However, because of their simplicity, spheres are only able to provide a qualitative description of granular material behaviour as they cannot provide enough particle interlocking and associated rotational resistance. Different approaches have been used to obtain quantitatively reliable results [30][31][32] The approach adopted in this study consists of using an unbreakable clump of ten spheres arranged in a tetrahedral shape ( Fig. 1) providing enough interlocking between the particles and hence shear strength under triaxial tests. Clumps of single size were used to model the uniformly distributed 1/3 scaled ballast (d 50 = 13.33 mm, c u = 1.5) used by Aingaran [27]. The diameter of each large sphere within the body of the clump is 5.4 mm so that the length of the longest axis of main clump body is the same as the d 50 size for the experimental particle. The asperities consisting of smaller spheres (2 mm diameter) are bonded to the main tetrahedral clump body using breakable bonds. These bonds are a feature of the software PFC3D and act as a cementing cylinder characterised by normal and shear stiffness coefficients k pbn and k pbs , normal and shear strengths σ n and σ s and bond diameter equal to the average diameter of the two contacting spheres multiplied by a factor r pb (equal to 1 in all simulations here) [29]. A similar model using eight asperities showed that the breaking of asperities can lead to a reduction of dilation at high confining pressure in triaxial tests [15]. In the present study, the effect of the number of asperities has been investigated using 1, 3, 6 and 8 asperities as well as different asperity bond strengths, in order to model the behaviour of the scaled ballast over a range of confining pressures. The asperities are randomly created at the vertices of the tetrahedral clump for different numbers of asperities (Fig. 3).

Triaxial test conditions
The laboratory triaxial tests used as a reference here [27] were performed on 1/3 scaled ballast obtained by crushing and sieving of standard Network Rail ballast in the United-Kingdom. This was done as part of the EPSRC collaboration project 'Railway Track for the twenty-first century'. The cylindrical sample has a 15 cm diameter and is 30 cm high. To reduce calculation time in the DEM simulations by decreasing the number of particles, the simulation sample has a 9 cm diameter and is 18 cm high. A cylindrical vertical wall is used to simulate a flexible membrane [15] and two horizontal walls were placed at the top and bottom of the sample respectively. The walls were set to be frictionless with zero shear stiffness during each simulation including both sample preparation procedure and loading. The two horizontal walls had the same normal stiffness as the particles (10 9 N/m); the normal stiffness of the cylindrical vertical wall (10 8 N/m) was set to be one tenth of the normal stiffness of the particles in order to simulate the effect of the flexible membrane. As the influence of the sample to particle size ratio on results becomes negligible when the size ratio is approximately more than 6 [33], this smaller sample size should not significantly affect the results. The experimental samples had initial voids ratios between 0.68 and 0.70. The samples in DEM were generated using the radius expansion method: an assembly of spheres was first generated within the cell randomly without overlapping; these spheres were initially generated with an artificially small size and then expanded progressively to create the desired voids ratio of the sample. The system was then cycled to equilibrium with all walls fixed and the spheres were then replaced by the clumps (Fig. 2) with random orientation. Finally the system was cycled to equilibrium again to release the high contact forces produced by particle replacement (Fig. 4). The samples were then sheared at low and high confining pressures 15 and 200 kPa respectively at a speed of 0.1 m/s. Preliminary simulations showed that this shearing speed, chosen as high as possible to reduce DEM calculation time, did not have a significant impact on the results (i.e. reducing the speed by a factor of two, say, would have negligible effect on the results. The confining pressure was applied using a cylindrical wall and a servo-function to Fig. 4 Assembly of ten-ball clumps with asperities bonded in the triaxial cel Wall friction 0 vary the radial velocity and two moving platens at the top and bottom of the sample in order to maintain a constant radial stress. Table 1 lists the basic sample parameters which are used in this study as well as in previous studies of DEM modelling of triaxial tests on ballast particles [15,34]. The large spheres forming the clump and asperities have the same stiffness, given in Table 1. However these are bonded to each other and in compression the contact stiffness acts in parallel with the parallel bond stiffness. In tension there is no "contact bond", so the stiffness in tension is due only to the parallel bond.

Simulations without crushing
The purpose of the results presented here is to show a concise way to obtain a mobilised friction angle φ mob and volumetric deformation ε v which are quantitatively consistent with experimental results and thereby providing a micro mechanical explanation for the observed behaviour. Ten simulations in total are presented in this study, and the sample parameters are listed in Table 2. As real and model ballast particles have different shapes, the initial voids ratio of the model sample is different to the experimental one. Therefore the authors firstly tried alternative initial voids ratios to choose that which provides a similar peak shear strength, using clumps without asperities; these are the first three simulations in Table 2. Figure 5 shows the mobilised friction angle φ mob and volumetric deformation ε v from triaxial tests on tetrahedral clumps without asperities, at low and high confining pressures of 15 and 200 kPa respectively, for a range of alternative initial voids ratios. As Fig. 5 shows, the denser sample results in a higher peak mobilised friction angle at low confining pressure; for the volumetric deformation ε v , a higher voids ratio leads to less dilation at both stress levels. The sample with an initial voids ratio of 1.02 provides a much lower peak shear strength and lower maximum dilation rate compared to the   Although the denser sample with an initial void ratio of 0.88 offered a very close volumetric curve, its shear strength is not as close as the densest sample to the experimental result. Therefore, the initial voids ratio of 0.72 was chosen as it provides the closest peak strength at both stress levels; the aim here was simply to choose an initial voids ratio as a benchmark for future simulations with different numbers of asperities and asperity strengths over the entire range of confining pressures used in the experimental tests.

Simulations with breakage
Further simulations were conducted with 8 breakable asperities (simulations 4-7). The bond parameters are given in Table 2. The parallel bond stiffnesses k pbn and k pbs were chosen as the same values as in Lu and McDowell [15] and the bond strength was given a range of values to investigate its effect. Figure 6 shows the mobilised angles of friction and volumetric behaviour for simulation 4 compared to simulation 1 with no asperities. Introducing breakable asperities reduces the dilation and leads to a better modelling of the volumetric deformation at both stress levels because of the introduction of particle abrasion. In addition, the volumetric deformation at high confining pressure is affected more by the presence of the crushable asperities than that at low confining pressure, due to more broken asperities at the higher confining pressure. However, there is roughly 10 • difference in the peak mobilised friction angles at low and high confining pressure according to Fig. 6, which is only observed in the experimental results on full size ballast [2] but not on scaled size ballast [27]; this must because at the stress levels being considered, the smaller scaled ballast particles are statistically stronger [35,36] and influenced less by the range of applied confining pressures. In order to model better the mobilised strength of the scaled ballast under the applied confining pressures, the effect of changing the bond strengths for the 8-asperity clump (simulations 5-7) and the number of asperities (simulations 8-10, with the same bond strength as simulation 4) will now be investigated.  Figure 7 shows the effect of changing the bond strength for a clump with 8 asperities. It can be seen that the peak mobilised friction angle increases and the sample becomes more dilative with increasing asperity bond strength at both stress levels. Figure 7c, f shows the percentage of broken asperity bonds during monotonic loading. As expected, the percentage of broken asperities increases with decreasing bond strength. It is clearly observed that the number of broken asperities affects both the mobilised friction and the volumetric deformation of the sample, which is consistent with Lu and McDowell [15]: the higher the number of broken asperities during shearing, the lower is the mobilised friction angle and the less dilative the sample. However, it appears that no single bond strength for the 8-asperity clump is capable of capturing the mobilized strength and dilation at both 15 and 200 kPa. Thus the effect of the number of the crushable asperities on the behaviour was investigated in simulations 8-10.
If the peak strength and dilation of the sample is governed by the behaviour of the asperities, then it would be anticipated that the addition of more small surface asperities would increase the peak strength and dilation at low stress, but that effect would reduce at high stress levels due to the breakage of these asperities. This is clearly shown in Fig. 8. It is evident that the use of three asperities (simulation 9) is able to capture the essential behaviour of the scaled ballast over the range of applied confining pressures. One problem with the simulations, however, is establishing the effect of the asperities on the ultimate or critical state angle of friction because in each of the simulations the sample has not reached a critical state; indeed they are still dilating and in each case at an approximately constant rate with axial strain, giving a wide peak strength as a function of axial strain. Thus, in order to simulate conditions at large strains to establish the critical state angles, all that is required is the angle of repose in each case-in this test particles are indeed rolling down the soil slope at large shear strains. If the size of the asperities is sufficiently small, this should have a negligible effect on the ultimate friction angle φ crit , but still be sufficient to provide additional interlock and higher peak strength at low stress levels and low shear strains. Figure 9 shows a simple heap test for each clump shape in order to determine the effect of the number of asperities on the angle of repose. Figure 9 shows an approximately equal angle of repose in each case. This is to say, by changing the number of surface asperities, the tetrahedral clump samples exhibit different peak strengths at low confining pressure, but still provide an approximately constant ultimate friction angle φ crit .Therefore it is possible to achieve a quantitative modelling of the peak mobilised friction angle and volumetric strain at each stress level for scaled ballast by considering the number of surface asperities, without adversely affecting the modelling of ultimate conditions.

Conclusion
This paper demonstrates that it is possible to reproduce the behaviour of ballast particles at different scales in DEM using tetrahedral clumps with breakable asperities. It has been shown that introducing asperities increases the peak shear strength and dilation of sample but the effect reduces with increasing confining pressure due to particle abrasion (breaking of asperities). By examining the effect of asperity strength and the number of asperities on the mobilised peak strength and dilation, it has been possible to model the behaviour of scaled ballast over a range of confining pressures under triaxial conditions. However, under such conditions the axial strains are insufficient for the samples to reach ultimate conditions and so the effect of the number of asperities on the angle of repose was determined to ensure that not only could the peak strength and volumetric behaviour be correctly captured, but that the critical state angle of friction was not significantly affected. The results have shown the importance of the modelling of asperity breakage if the correct peak strength and associated dilation together with ultimate conditions are to be correctly captured under different confining pressures. Future work will examine the effect of a flexible membrane on the triaxial test simulation results, although this is computationally much more intensive.