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

Authors: Hongyi Zhao; Jing Chen
Journal: Advances in Materials Science and Engineering
 (2020)

---

## Abstract

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 different normal pressures. The results show that the contact between ballast particles intensifies 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 different 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 identified 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.


---

## Body

## 1. 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. The 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. The 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. This 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–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. The 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–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). The 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.


## 2. Numerical Simulations

### 2.1. Discrete Element Modeling

The 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–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–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. The 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. This continues until the specified requirements are met.Figure 1
A schematic illustration of PFC3D.


### 2.2. The Linear Contact Model

The 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. The material parameters used to define the model are the normal and tangential spring constantskn and ks, the normal and tangential critical damping ratios βn and βs, and the coefficient of friction μ. They appear in the force-displacement relations for normal and tangential force calculations as(1)Fn=minFnl0+knΔδn+2βnmcknδ˙n,0,Fs=minFsl0−ksΔδs+2βsmcksδ˙s,−μFn,withmc=m1m2m1+m2ball−ball,m1ball−facet,where Fln0 and Fls0 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).


### 2.3. Generation and Calibration of Ballast Particles

In this research, the six ballast particles shown in Figure2 were created using the method proposed by Chen et al. [18]. Each ballast clump consists of around 100 spheres. The 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). The 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). These tests were repeated three times to reduce any experimental error. The angle of repose of the ballast in these experiments was 32.5°. The 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. The 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, kn, ks, knw, and ksw, have negligible influence on the formation of ballast piles. In this regard, the values of kn, ks, knw, and ksw are directly inherited from the existing research. The 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.Figure 2
The six ballast clumps used in DEM models.Figure 3
The cylinder lifting test: (a) loaded ballast in the cylinder; (b) ballast pile under gravity; (c) DEM model; (d) measurement cone in PFC3D.
(a)(b)(c)(d)Table 1
The parameters of the linear contact model.
ParametersValueParticle density (kg/m3)2700Contact normal and shear stiffness,kn and ks (N/m)1.0 × 107Contact normal and shear stiffness of wall-particleknw and ksw (N/m)5.2 × 106Wall-particle coefficient of frictionµw0.3Particle-particle coefficient of frictionµp0.7Dampc0.2


### 2.4. DEM Models of Direct Shearing

The 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 Figure2 were randomly replicated in a shear box. The 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 Chinese Railway Ballast Standard (TB/T 2140-2008); the size of the ballast aggregates and some critical parameters are shown in 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. The porosity of the aggregates after compaction was measured with a measuring ball [26] and reached 0.43. The 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%. The 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 4
The size distribution of ballast aggregates in DEM.


## 2.1. Discrete Element Modeling

The 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–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–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. The 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. This continues until the specified requirements are met.Figure 1
A schematic illustration of PFC3D.


## 2.2. The Linear Contact Model

The 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. The material parameters used to define the model are the normal and tangential spring constantskn and ks, the normal and tangential critical damping ratios βn and βs, and the coefficient of friction μ. They appear in the force-displacement relations for normal and tangential force calculations as(1)Fn=minFnl0+knΔδn+2βnmcknδ˙n,0,Fs=minFsl0−ksΔδs+2βsmcksδ˙s,−μFn,withmc=m1m2m1+m2ball−ball,m1ball−facet,where Fln0 and Fls0 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).


## 2.3. Generation and Calibration of Ballast Particles

In this research, the six ballast particles shown in Figure2 were created using the method proposed by Chen et al. [18]. Each ballast clump consists of around 100 spheres. The 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). The 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). These tests were repeated three times to reduce any experimental error. The angle of repose of the ballast in these experiments was 32.5°. The 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. The 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, kn, ks, knw, and ksw, have negligible influence on the formation of ballast piles. In this regard, the values of kn, ks, knw, and ksw are directly inherited from the existing research. The 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.Figure 2
The six ballast clumps used in DEM models.Figure 3
The cylinder lifting test: (a) loaded ballast in the cylinder; (b) ballast pile under gravity; (c) DEM model; (d) measurement cone in PFC3D.
(a)(b)(c)(d)Table 1
The parameters of the linear contact model.
ParametersValueParticle density (kg/m3)2700Contact normal and shear stiffness,kn and ks (N/m)1.0 × 107Contact normal and shear stiffness of wall-particleknw and ksw (N/m)5.2 × 106Wall-particle coefficient of frictionµw0.3Particle-particle coefficient of frictionµp0.7Dampc0.2


## 2.4. DEM Models of Direct Shearing

The 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 Figure2 were randomly replicated in a shear box. The 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 Chinese Railway Ballast Standard (TB/T 2140-2008); the size of the ballast aggregates and some critical parameters are shown in 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. The porosity of the aggregates after compaction was measured with a measuring ball [26] and reached 0.43. The 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%. The 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 4
The size distribution of ballast aggregates in DEM.


## 3. Results and Discussion

### 3.1. Macroscopic Characteristics

Figure5 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. The peak shear strength is reached when εs reaches 8%. There 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. The 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.Figure 5
The shear stress and normal strain of ballast aggregates under various normal pressures.Figure 6
The contacts between ballast particles atεs of 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)


### 3.2. Microscopic Characteristics

#### 3.2.1. Particle Contacts

The 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. The total number of contacts and some statistical values of four samples are listed in Table 2. The 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. This 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.Table 2
The statistical values of contacts between ballast particles atεs of 10%.
σnNcFmaxNFaveNFmedN15 kPa7679631.9892.9828.1335 kPa7703630.82112.9643.5755 kPa7854634.53135.6164.0375 kPa7959634.10139.0270.28To 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 Figures 7–9. 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. The contact information within each bin limit was collected. The 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. The 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:(2)Eθ=12π1+accos2θ−θc,f¯nθ=f¯01+ancos2θ−θn,f¯tθ=f¯0−atsin2θ−θt,f¯0=∫02πf¯θndθ,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 ac, an, at 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.Figure 7
Rose diagrams of contact normal distribution atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figure 8
Rose diagrams of contact normal force distribution atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figure 9
Rose diagrams of contact tangential force distribution atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figures7 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. The values of ac 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. The orientation of the contact normal forces is similar to the contact normal, whereas θn is slightly lower than θc with values ranging from 32.3° to 37.3°. The 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 at, an and at 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. These intensified contact forces provided sufficient shear resistance as the shearing propagated.


#### 3.2.2. Particle Rotation

In PFC3D, the particles could rotate about theX-, 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 Y-axis. 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. The 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. The 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.Figure 10
The rotation of particles atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.Figures11–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. The 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%. This was probably due to the positive rotation of particles in local areas being thwarted by their surroundings, which reduced their maximum positive rotations.Figure 11
The total number of particle rotations atεs of 2%, 6%, and 10%: (a) about the negative y-axis; (b) about the positive y-axis.
(a)(b)Figure 12
The maximum, medium, and average values of particle rotation about the negativey-axis: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figure 13
The maximum, medium, and average values of particle rotation about the positivey-axis: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)


#### 3.2.3. Particle Displacement

Figure14 shows diagrams of the cumulative displacement of particles from the beginning till the end of shearing. The 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. The 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. This 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]. This 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 dilation of the system is comprehensively influenced by particles moving upwards and downwards. The greater downward movement of particles in ‘Region A’ balances out the volume dilation induced by particles moving up in ‘Region B’; therefore the change in volume was smaller under higher normal loading pressures.Figure 14
The displacement of particles: (a)σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)


## 3.1. Macroscopic Characteristics

Figure5 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. The peak shear strength is reached when εs reaches 8%. There 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. The 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.Figure 5
The shear stress and normal strain of ballast aggregates under various normal pressures.Figure 6
The contacts between ballast particles atεs of 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)


## 3.2. Microscopic Characteristics

### 3.2.1. Particle Contacts

The 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. The total number of contacts and some statistical values of four samples are listed in Table 2. The 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. This 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.Table 2
The statistical values of contacts between ballast particles atεs of 10%.
σnNcFmaxNFaveNFmedN15 kPa7679631.9892.9828.1335 kPa7703630.82112.9643.5755 kPa7854634.53135.6164.0375 kPa7959634.10139.0270.28To 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 Figures 7–9. 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. The contact information within each bin limit was collected. The 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. The 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:(2)Eθ=12π1+accos2θ−θc,f¯nθ=f¯01+ancos2θ−θn,f¯tθ=f¯0−atsin2θ−θt,f¯0=∫02πf¯θndθ,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 ac, an, at 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.Figure 7
Rose diagrams of contact normal distribution atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figure 8
Rose diagrams of contact normal force distribution atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figure 9
Rose diagrams of contact tangential force distribution atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figures7 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. The values of ac 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. The orientation of the contact normal forces is similar to the contact normal, whereas θn is slightly lower than θc with values ranging from 32.3° to 37.3°. The 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 at, an and at 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. These intensified contact forces provided sufficient shear resistance as the shearing propagated.


### 3.2.2. Particle Rotation

In PFC3D, the particles could rotate about theX-, 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 Y-axis. 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. The 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. The 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.Figure 10
The rotation of particles atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.Figures11–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. The 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%. This was probably due to the positive rotation of particles in local areas being thwarted by their surroundings, which reduced their maximum positive rotations.Figure 11
The total number of particle rotations atεs of 2%, 6%, and 10%: (a) about the negative y-axis; (b) about the positive y-axis.
(a)(b)Figure 12
The maximum, medium, and average values of particle rotation about the negativey-axis: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figure 13
The maximum, medium, and average values of particle rotation about the positivey-axis: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)


### 3.2.3. Particle Displacement

Figure14 shows diagrams of the cumulative displacement of particles from the beginning till the end of shearing. The 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. The 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. This 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]. This 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 dilation of the system is comprehensively influenced by particles moving upwards and downwards. The greater downward movement of particles in ‘Region A’ balances out the volume dilation induced by particles moving up in ‘Region B’; therefore the change in volume was smaller under higher normal loading pressures.Figure 14
The displacement of particles: (a)σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)


## 3.2.1. Particle Contacts

The 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. The total number of contacts and some statistical values of four samples are listed in Table 2. The 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. This 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.Table 2
The statistical values of contacts between ballast particles atεs of 10%.
σnNcFmaxNFaveNFmedN15 kPa7679631.9892.9828.1335 kPa7703630.82112.9643.5755 kPa7854634.53135.6164.0375 kPa7959634.10139.0270.28To 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 Figures 7–9. 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. The contact information within each bin limit was collected. The 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. The 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:(2)Eθ=12π1+accos2θ−θc,f¯nθ=f¯01+ancos2θ−θn,f¯tθ=f¯0−atsin2θ−θt,f¯0=∫02πf¯θndθ,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 ac, an, at 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.Figure 7
Rose diagrams of contact normal distribution atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figure 8
Rose diagrams of contact normal force distribution atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figure 9
Rose diagrams of contact tangential force distribution atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figures7 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. The values of ac 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. The orientation of the contact normal forces is similar to the contact normal, whereas θn is slightly lower than θc with values ranging from 32.3° to 37.3°. The 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 at, an and at 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. These intensified contact forces provided sufficient shear resistance as the shearing propagated.


## 3.2.2. Particle Rotation

In PFC3D, the particles could rotate about theX-, 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 Y-axis. 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. The 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. The 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.Figure 10
The rotation of particles atεs of 0% and 10%: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.Figures11–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. The 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%. This was probably due to the positive rotation of particles in local areas being thwarted by their surroundings, which reduced their maximum positive rotations.Figure 11
The total number of particle rotations atεs of 2%, 6%, and 10%: (a) about the negative y-axis; (b) about the positive y-axis.
(a)(b)Figure 12
The maximum, medium, and average values of particle rotation about the negativey-axis: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)Figure 13
The maximum, medium, and average values of particle rotation about the positivey-axis: (a) σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)


## 3.2.3. Particle Displacement

Figure14 shows diagrams of the cumulative displacement of particles from the beginning till the end of shearing. The 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. The 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. This 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]. This 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 dilation of the system is comprehensively influenced by particles moving upwards and downwards. The greater downward movement of particles in ‘Region A’ balances out the volume dilation induced by particles moving up in ‘Region B’; therefore the change in volume was smaller under higher normal loading pressures.Figure 14
The displacement of particles: (a)σn = 15 kPa; (b) σn = 35 kPa; (c) σn = 55 kPa; (d) σn = 75 kPa.
(a)(b)(c)(d)


## 4. Limitations

(1)
The ballast clumps used in this research are unbreakable. The 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. The 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.


## 5. 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. The parameters of the contact model were calibrated by matching the predicted angle of repose for ballast in the cylinder lifting test. These 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. The key findings are summarized as follows:(1)
The 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)
The 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 negativeY-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. The overall volume dilation of the aggregate is caused by the combination of particle rising and falling in the upper shear box.



---

*Source: 3404208-2020-11-23.xml*