A Numerical Study of Railway Ballast Subjected to Direct Shearing Using the Discrete Element Method

Railway ballast is a coarse granular material used to carry train loads and provide drainage for the rail tracks. This study presents numerical explorations of the mechanical performance of ballast aggregates subjected to direct shear tests. The discrete element method (DEM) was used to investigate the microscopic characteristics of ballast aggregates during shearing while considering contact distribution, particle rotation, and particle displacement. By testing the angle of repose of ballast aggregates, the parameters for the DEM contact model could be calibrated. Four specimens were prepared and then subjected to diﬀerent normal pressures. The results show that the contact between ballast particles intensiﬁes in terms of the amount and magnitude as the normal pressure increases. A Fourier analysis was applied to investigate the anisotropy of contact normal and the contact forces for ballast aggregates at diﬀerent shearing phases. The rotational and translational movements of ballast particles were investigated, and this investigation revealed that particle rotation gradually increased as the shearing propagated. Four regions in the aggregates were identiﬁed according to the translational pattern of ballast particles. The results of this research provide an in-depth analysis of microscopic characteristics from a particulate scale.


Introduction
Railways are one of the major means of transportation worldwide because they play a crucial role in influencing the economic development of a country. Ballasted tracks make up a large proportion of existing railway lines due to their reliable operation, low construction cost, and simple maintenance requirements. Ballast usually consists of medium-to-coarse crushed rock fill that provides a platform for free drainage in railways. e upper loads generated by passing trains are also dissipated to an acceptable level and are transmitted to the subballast and subgrade via the ballast layer. It is very important to understand the mechanical properties of ballast aggregates theoretically and practically.
Direct shear tests are a traditional and widely used method of studying the mechanical properties of geomaterials [1,2]. While a lot of existing research has used direct shear tests on ballast aggregates [3,4], and the mechanical properties of ballast have been fully examined in the laboratory or in the field assessment [5], the limitations of experimental techniques mean that only the macroscopic characteristics of ballast aggregates can be investigated, but limiting research to the macrolevel is not good enough, especially for coarse granular materials such as ballast. e problem is individual ballast particles change continually during their service life due to abrasion and inevitable breakdown. It has been suggested, and investigated, that the degradation of ballast is associated with insufficient drainage and differential settlement, issues that lead to high maintenance costs but safety risks. is is why ballast aggregates need to be examined from a microscopic perspective.
To explore the microscopic performance of ballast, many studies used the discrete element method (DEM) proposed by Cundall and Strack in 1979 [6], and it has been successfully applied by many researchers in exploring the behavior of ballast aggregates [7][8][9][10]. In DEM, particles are treated as an individual body, which suits the granular characteristics of ballast. Unlike sand or coal fines, the morphologies of ballast particles are extremely irregular, as well as the friction induced by the surface texture of granular material. e interlocking effect due to the angularity of grains is much more pronounced between ballast particles. It is acknowledged that particle shape is a key factor affecting the macromechanical behavior of coarse aggregates [11][12][13][14][15][16], which is why existing research has emphasized particle shapes in DEM modeling. Lu and McDowell [17] carried out large-scale box test simulations with different particles and found that models with clumps exhibited more realistic results than with simple spheres. By carrying out a series of simulated calibrations in DEM, Coetzee found that clumps with complex shapes performed better in validation tests [14], but even so, particles with simple shapes are still preferred in prevailing research because of computational efficiency. Despite the satisfying results that have been obtained with simplified particle shapes, large discrepancies would emerge in numerical simulations when irregular shaped particles are involved [18]. It is therefore imperative to explore the microscopic mechanical characteristics of the aggregates using particles that resemble real ballast shapes.
In this research, the behavior of ballast aggregates subjected to direct shearing was explored numerically using DEM. Particles that resemble real ballast were created in the commercial DEM software known as Particle Flow Code 3D (PFC3D). e parameters of the DEM model were calibrated through a series of angle of repose tests, and then the model was sheared under various normal pressures so that the microscopic performance of the aggregates could be investigated.

Discrete Element Modeling.
e unique feature of DEM is that it can capture the behavior of individual particles and their interaction from a particulate scale [19], and unlike traditional numerical methods such as the finite element method (FEM) [20][21][22], it can analyze large discontinuous deformation. Boosted by advanced commercial software such as PFC3D, which was developed by Itasca, DEM has been successfully applied in various research disciplines [23][24][25]. Despite the direct shear response can be regarded as planar, there is no coupling between the translational and rotational degrees of freedom using 2D DEM simulations as stated by O'Sullivan [19], which may lead to the incorrect analysis on particle rotation. In this regard, the 3D DEM simulation was adopted in the current research to explore the mechanical behavior of ballast aggregates. As Figure 1 shows, the basic element in PFC3D is a sphere, from which complex particles can be created. Clumped spheres act as a rigid body while their in-between interactions are ignored. Certain contact models must be specified to calculate the contact forces between two bodies. Based on Newton's second law, the acceleration of individual bodies is calculated, and the particle velocity is integrated. e incremental displacement of individual bodies over a small time step can be determined, and then the position of each particle at the end of the time increments is updated; the whole system moves forward to the next calculation sequence. is continues until the specified requirements are met.

e Linear Contact Model.
e force transmitted by ball-to-ball contact can be described using combinations of linear springs and dashpot elements in parallel, where the linear springs provide elastic (no-tension) and frictional behavior and the dashpot element provides viscous behavior. e material parameters used to define the model are the normal and tangential spring constants k n and k s , the normal and tangential critical damping ratios β n and β s , and the coefficient of friction μ. ey appear in the forcedisplacement relations for normal and tangential force calculations as where F l n0 and F l s0 are the linear normal and shear force at the beginning of the time step, Δδ n is the relative normal-displacement increment, Δδ s is the relative shear-displacement increment, δ n and δ s are the relative translational velocities normal and tangential to the contact plane, and m (b) is the mass of the body (b).

Generation and Calibration of Ballast Particles.
In this research, the six ballast particles shown in Figure 2 were created using the method proposed by Chen et al. [18]. Each ballast clump consists of around 100 spheres. e linear contact model described in Section 2.2 [26] was utilised for contact resolution. To calibrate the parameters needed for the contact model, a cylinder lifting test (CLT) was carried out to measure the angle of repose for ballast aggregates. In these laboratory tests, 100 kg of ballast was loaded into a 400 mm diameter cylinder, as shown in Figure 3(a). e cylinder was raised at a constant velocity of 5 mm/s, and then the ballast particles fell down and formed a pile (Figure 3(b)).
To determine the angle of repose for this ballast, the angles of the piles were measured using the image-based method shown in Figure 3(b). ese tests were repeated three times to reduce any experimental error. e angle of repose of the ballast in these experiments was 32.5°. e parameters of the linear contact model were calibrated by establishing a corresponding DEM model (Figure 3(c)). After being generated inside the cylinder, the ballast clumps settled down under gravity, as per the experiment.
e whole system was then cycled to reach the state of equilibrium where the unbalanced force ratio reached 1e −5 , and then the cylinder was raised up at the same velocity to form piles of ballast. Using the subroutines developed by the authors in PFC3D, the smallest cone that could contain all the particles was generated as shown in Figure 3(d). As investigated by Nakashima et al. [27], the contact normal and tangential stiffness of ballast particles and walls, k n , k s , k nw , and k sw , have negligible influence on the formation of ballast piles. In this regard, the values of k n , k s , k nw , and k sw are directly inherited from the existing research. e friction coefficient was calibrated by generating the ballast repose angle similar to that measured in the CLTs with DEM models. It should be noted that the viscous parameter, c, was included in the current model to reflect the loss of energy due to collisions between the ballast particles. Table 1 summarizes the parameters of the current model.

DEM Models of Direct Shearing.
e calibrated DEM model was used to investigate the microscopic behavior of the ballast aggregates involved in a series of direct shearing tests. In this model, the six clumps of ballast shown in Figure 2 were randomly replicated in a shear box. e shear box is 600 mm wide, 600 mm long, and 500 mm high. Each ballast particle was scaled so that the size of the aggregates would comply with the  Figure 4. Note here that the ratio of D50 to the width of the shearing box is more than 6 so that the boundary effect during shearing would be eliminated. After being generated in the box, the aggregates were compacted under a normal pressure of 100 kPa to reach a dense state. e porosity of the aggregates after compaction was measured with a measuring ball [26] and reached 0.43. e box was divided into two equal halves, and the bottom shearing box moved at a constant velocity of 3 mm/min. A baffle wall was installed at the shearing plane to stop the ballast from falling from the upper half of the box, as shown in Figure 4. Four normal pressures, 15, 35, 55, and 75 kPa were applied onto the specimens. During shearing, the top wall can move in a vertical direction, and its movement is controlled by a servocontrol system to maintain a constant normal pressure. All the specimens were sheared at a shear strain of 10%. e total force acting on the bottom box along the shearing direction was recorded as the shearing force. In the simulation, the movement of ballast particles was also recorded for later analysis. Figure 5 shows the relationship between shear stress τ s and shear strain ε s . As ε s developed, the shear stress of all samples gradually increased in a declining gradient. e peak shear strength is reached  2700 Contact normal and shear stiffness, k n and k s (N/m) 1.0 × 10 7 Contact normal and shear stiffness of wall-particle k nw and k sw (N/m) 5.2 × 10 6

Macroscopic Characteristics.
Wall-particle coefficient of friction µ w 0.3 Particle-particle coefficient of friction µ p 0.7 Damp c 0.2 when ε s reaches 8%. ere were some fluctuations in the shear stress after it reached its peak value. As expected, the shear stress τ s was larger at a higher normal stress σ n . With the onset of shearing, the interaction between ballast was initiated by the angularity of the particles, which increased in the shear strength; this interaction became stronger at higher external loading pressures. At the end of shearing, the ballast particles tended to move up and over their neighbours, which attenuated the shearing resistance provided by particle interaction and slowed down the shear stress. However, as shearing continued, some new interactions developed and the shear stress increased slightly. For example, the shear stress continued to grow at higher normal pressures where σ n was 55 kPa and 75 kPa. e evolutions of normal strain ε n over shear strain ε s for the four samples are shown in Figure 6. Every sample experienced a short volumetric compaction stage at the beginning of shearing, but then the normal strain became positive and gradually increased until the shearing process ended. With this increase of normal pressure, the samples exhibited smaller volumetric dilations because the ballast particles move up and down under higher normal pressure.

Particle Contacts.
e contacts between particles at the end of shearing (i.e., 10% of ε s ) for ballast aggregates under various normal pressures are shown in Figure 6. In this figure, the contacts are shown as cylinders that connect the centroids of two contacting pebbles where the size of each cylinder denotes the magnitude of the contact force. Figure 6 also shows that the amounts and the magnitudes of the contacts are enhanced by the increase of normal pressure. Meanwhile, there is an obvious band in the samples where most strong particle contacts are located; this band developed from the right side of the top box to the left side of the bottom box and inclined at an angle to the horizontal plane. e total number of contacts and some statistical values of four samples are listed in Table 2. e normal pressures increased and so did the number of particle contacts. However, the maximum contact force remained almost the same in all four samples, but there was a larger average and medium value at higher normal pressures. is increasing external normal pressure enhanced the interaction between ballast particles, leading to an increase in the total number of contacts. As more particle contacts formed in the aggregates to bear the external forces, the local concentrations of force dispersed and the maximum contact forces in the four samples remained at the same level. However, this increasing average and medium value indicated that more particle contacts have larger magnitudes and therefore greater shearing resistance.
To explore the distribution of normal contact and contact forces inside the samples quantitatively, rose diagrams of these distributions at 0% and 10% of ε s are shown in . Note that all the contacts and corresponding contact forces have been projected onto the XZ plane for a better illustration. A bin angle of 10°was predefined, and 36 bins were specified, as shown in Figures 7-9. e contact information within each bin limit was collected. e sizes of the bins in the rose histograms are proportional to the     probable density of contacts falling into each bin angle. Notably, the contact normal and the contact forces were normalised with respect to the total number of contacts and the average normal contact force.
e Fourier series approach was used for the approximation presented as red lines in Figures 7-9. Based on the work by Rothenburg and Bathurst [28], the distribution of contact normal and the contact force can be fitted by the Fourier series as follows: where E(θ), f n (θ), and f t (θ) are the probable density of the contact normal where the normal and tangential contact force fall within the specified bin angles and a c , a n , a t and θ c , θ n , θ t are the anisotropy and the principal orientation angles of the contact normal, the normal and tangential contact forces, and f 0 is the average normal contact force. Figures 7 and 9 show that all four samples have similar contact distributions irrespective of whether there is 0% or 10% of ε s . At the beginning of shearing, most particle contacts are distributed vertically, and during this stage, the contacts between ballast particles need to bear the external normal pressures and their self-weight. As shearing propagates, the major principal direction of particle contacts rotates clockwise and inclines at angles ranging from 42.9°to 47.6°to the horizontal plane. e values of a c decreased as shearing developed in all the samples because that part of the vertically distributed contacts rotates towards the horizontal plane, so there is less anisotropy of particle contacts in a vertical direction. e orientation of the contact normal forces is similar to the contact normal, whereas θ n is slightly lower than θ c a n = 0.169 θ n = 73.1°τ e contact tangential force has four symmetrical peaks, as shown in Figure 9, as observed in many existing research [13,28,29]. Unlike the various trends of a t , a n and a t increased from the beginning to the end of shearing of the four samples, but during this process, the amount and magnitude of the contact forces along certain angles intensified, leading to an increased level of anisotropy of the contact normal and tangential forces. ese intensified contact forces provided sufficient shear resistance as the shearing propagated.

Particle Rotation.
In PFC3D, the particles could rotate about the X-, Y-and Z-axis; this rotation is usually quantified as the Euler angle. However, only rotation in the Z plane, i.e., about Y-axis, is analyzed in this research. Note that the sign of particle rotation obeys the right-hand rule, so the negative value denotes rotation about the negative Yaxis. Figure 10 shows the rotation of samples at 2%, 6%, and 10% of ε s , where the rotation of each ballast particle is represented as a circle located in the same position as the particle. e circles are scaled in size according to the Euler angle of their corresponding particle rotations; negative rotation is shown as solid circles, and positive rotation as hollow circles. Note that particle rotation intensifies as shearing propagates. Rotation initially emerges at the back part of the upper shearing box and the front of the bottom box. e particles near the shearing plane gradually rotate to form an obvious band that connects the upper back and the bottom front. As Figure 10 shows, particles within the band have mainly negative rotation, and they are much pronounced than the positive ones. As the bottom box moves towards the positive X-axis to exert a shearing force onto the whole system, shearing resistance would be introduced by the negatively rotating particles; however, regardless of the variations in normal pressure, there was no significant change in the pattern of particle rotation. Figures 11-13 show some statistics of particle rotation. For instance, Figure 11 shows that the total number of particle rotations about the negative Y-axis gradually increased when ε s increased from 2% to 10%. However, the number of positive particle rotations decreased as shearing developed, which indicates that the rotation of some particles reversed during the shearing process. e maximum, medium, and average values of particle rotation (i.e., θ max , θ med , and θ ave ) about the negative and positive Y-axis of different samples are shown in Figure 12 and Figure 13. For particle rotation about the negative Y-axis, θ max , θ med , and θ ave increased as ε s increased from 2% to 10% in all four samples. With regards to positive particle rotation, θ med and θ ave were similar to negative particle rotation; however, θ max increased when ε s increased from 2% to 6% and then decreased when ε s was 10%.
is was probably due to the positive rotation of particles in local areas being thwarted by their surroundings, which reduced their maximum positive rotations. Figure 14 shows diagrams of the cumulative displacement of particles from the beginning till the end of shearing. e starting and finishing points of the vectors in Figure 14 are the positions of particles at 0% and 10% of ε s , respectively. As Figure 14 shows, all four samples exhibited a similar translational pattern whereby four distinct regions can be identified where particles have different movements. In the back part of the upper shearing box which is denoted as 'Region A,' the particles move anticlockwise. In 'Region A', particles from the upper back gradually lose support from the bottom during the shearing and fall; this fall is blocked by the horizontal baffle, so they begin to move upwards. e other particles in the upper shearing box (denoted as 'Region B') move obliquely upwards, and as particles in the rear left part of the bottom box move forward, they will extrude particles in 'Region B' and force them to move upwards. is extrusion effect is mainly due to the irregular shapes of ballast particles, as observed by Cui and O'Sullivan [30] and Liu et al. [10]. is moving up of particles in 'Region B' accounts for a large proportion of volume dilation during shearing, whereas in the bottom shearing box, most particles move horizontally and parallel to the direction of shearing, as shown in 'Region C'. Without interference by particles from the upper box, those located in 'Region C' are mainly pushed forward by the left wall. Besides 'Region C', a small number of particles in the upper front part of the bottom shearing box moved obliquely downwards into what is called 'Region D'. As the front wall moves, it would generate a local void space which allows the extruded particles at the front part of upper box to be squeezed downwards. As Figure 14 shows, particle displacement increases under higher normal pressures. However, it should be noted that the overall volume

Limitations
(1) e ballast clumps used in this research are unbreakable. e breakage of ballast aggregates under various loading conditions cannot be captured by the established DEM models; therefore, breakable ballast particles must be created, and their breakage behavior must be explored from a particulate scale. (2) Although water heavily influences the mechanical performance of ballast aggregates, it is not considered in this research. e interaction between ballast particles will be affected as water is induced into the aggregates. Innovative DEM models must be established to investigate the behavior of ballast aggregates with the inclusion of water.

Conclusion
A series of DEM models were created with a commercial code PFC3D to explore the microscopic behavior of ballast aggregates subjected to direct shear testing. Irregular ballast particles were created by clumping around 100 spheres together. e parameters of the contact model were calibrated by matching the predicted angle of repose for ballast in the cylinder lifting test. ese calibrated clumps of ballast were used to establish large-scale direct shear models. By shearing the ballast aggregates under four different loads, their macro-and microscopic performance was explored. e key findings are summarized as follows: (1) e macroscopic behavior of ballast aggregates under four different normal pressures was explored using the calibrated DEM models. As the normal pressure increased, the shear stress of the aggregates increased and the normal strain decreased. (2) e major particle contacts locate within a shear band. With the increase of normal pressure, the maximum contact forces and the anisotropy of particle contact normal and contact forces stay identical, while the number of particle contacts and the average and medium values of the contact force increase. (3) As shearing propagates, the contact normal rotates from vertical direction towards horizontal plane with a decrease in the anisotropic level, and the distributions of contact normal and tangential forces intensified along certain angles in terms of the amount and the magnitude. (4) As with the distribution of major contact forces, the particles mainly rotated within a band with a much more pronounced particle rotation about the negative Y-axis. As shearing develops, the amount, the average, and the medium of particle rotation gradually increase, and the maximum rotation increases while the positive rotation decreases after 6% of ε s . (5) Four regions can be identified according to the different translational patterns of particles. e overall volume dilation of the aggregate is caused by the combination of particle rising and falling in the upper shear box.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors Hongyi Zhao and Jing Chen declare no conflicts of interest.