Discrete element modelling of rubber-protected ballast performance subjected to direct shear test and cyclic loading

: The rubber-protected ballast (RPB) is made from natural ballast particles and crumb rubber particles. The crumb rubber is shredded waste tires. RPB was chosen to replace the ballast as it has higher resistance to breakage and abrasion. However, the static and dynamic performance of the RPB has not been conﬁrmed yet. Towards this end, experimental tests and numerical simulations were utilized to study the feasibility of RPB application. Direct shear tests (DSTs) were performed and a DST model and three-sleeper track model with the discrete element method (DEM) were built. The shear strength, settlement, displacement, and acceleration of the RPB were studied. The results show that the RPB has the advantage of increasing the force (stress) distribution and that the smaller crumb rubber size was more suitable for replacing the ballast particles.


Introduction
Railways provide the biggest universal network for rapid, economic, and safe passenger and freight transportation [1]. Currently, ballasted tracks are a widely-used infrastructure worldwide, as this type of track has a number of advantages compared with the slab track, including fast construction, low construction cost, high drainage capacity, and good noise and vibration absorption [2]. The ballast layer is one of the most essential components in a ballasted track, and it is built by placing ballast particles underneath and between the sleepers as the load-transmitting platform and for conveniently restoring track geometry [3,4]. Other ballast layer functions are also significant, such as providing adequate lateral and longitudinal resistance and sufficient drainage [5,6].
For the better performance of these functions, the ballast is carefully selected, considering such qualities as material, strength, and erosion resistance, however, after undergoing cyclic loadings, ballast particles become severely deteriorated (e.g., breakage and abrasion) [7,8]. This deterioration is exacerbated by increasing axle load (freight line) and train speed (passenger transport) [9]. Furthermore, the deteriorated ballast particles lead to shear-strength reduction and drainage failure [10,11]. Ballast fouling (powder and small sized particles) is traditionally considered a hazardous contamination to the track structure because it can increase permanent deformation and induce differential track settlement [12,13]. To combat ballast fouling, frequent maintenance is performed to restore the track geometry, and ballast replacement is needed as soon as the ballast layer fails [14].
Reducing ballast degradation is necessary for solving safety and economic problems [15,16]. More importantly, a new problem has occurred in some areas where the lack of high-quality parent rock has The main methodology applied in this study is numerical simulation with DEM models, since the DEM has been a viable tool for railway ballast simulation and successfully applied in many studies, e.g., [1,15,[28][29][30][31][32][33][34][35][36][37]. Because railway ballast is one kind of granular material, continuum models such as the finite element method or finite difference method are not able to present realistic ballast characteristics (e.g., movements, morphology and degradation). Dissimilarly, the DEM models can present not only the ballast characteristics but also ballast performances from the mesoscopic level (e.g., contact force chains, accelerations and displacements).
The commercial DEM software, Particle Flow Code in two dimensions (PFC2D) is utilized in this study. The calculation cycle performed in PFC2D is via a time-stepping algorithm that repeatedly applies (1) the Newton's second law to every particle, (2) a force-displacement relationship to every contact (3) and constant wall position updates. More specific explanations about the time-stepping algorithm can be found in [33].

Materials
The RPB applied in the laboratory DSTs was made from ballast particles glued with the CR particles using polyurethane. The ballast material is crushed volcanic basalt, provided by Tangshan Quarry, Hebei Province. Ballast material properties were examined based on the British standard, including the durability, mineralogy and morphology (size and shape) [38]. The properties of the ballast material, polyurethane and the CR particles are given in Table 1. Ballast and CR particles were washed and dried at room temperature before bonding them together. Two types of RPB were made with two different size ranges of CR particles, i.e., 0.0-0.25 mm and 2.5-5.0 mm, as shown in Figure  2a. The main methodology applied in this study is numerical simulation with DEM models, since the DEM has been a viable tool for railway ballast simulation and successfully applied in many studies, e.g., [1,15,[28][29][30][31][32][33][34][35][36][37]. Because railway ballast is one kind of granular material, continuum models such as the finite element method or finite difference method are not able to present realistic ballast characteristics (e.g., movements, morphology and degradation). Dissimilarly, the DEM models can present not only the ballast characteristics but also ballast performances from the mesoscopic level (e.g., contact force chains, accelerations and displacements).
The commercial DEM software, Particle Flow Code in two dimensions (PFC2D) is utilized in this study. The calculation cycle performed in PFC2D is via a time-stepping algorithm that repeatedly applies (1) the Newton's second law to every particle, (2) a force-displacement relationship to every contact (3) and constant wall position updates. More specific explanations about the time-stepping algorithm can be found in [33].

Materials
The RPB applied in the laboratory DSTs was made from ballast particles glued with the CR particles using polyurethane. The ballast material is crushed volcanic basalt, provided by Tangshan Quarry, Hebei Province. Ballast material properties were examined based on the British standard, including the durability, mineralogy and morphology (size and shape) [38]. The properties of the ballast material, polyurethane and the CR particles are given in Table 1. Ballast and CR particles were washed and dried at room temperature before bonding them together. Two types of RPB were made with two different size ranges of CR particles, i.e., 0.0-0.25 mm and 2.5-5.0 mm, as shown in Figure 2a. After the materials are prepared, RPB is produced with three components: crumb rubber, ballast particles and polyurethane. The manufacturing process includes three steps. Firstly, the polyurethane is sprayed on the ballast particles, and the polyurethane is like glue, which needs to mix two kinds of liquid. Secondly, the ballast particles (with polyurethane) are dropped into the crumb rubber chip, as shown in Figure 2b. Finally, the finished particles are left on the table to be dried up. The RPB mass ratio of the ballast particles to the CR particles is 5:0.16 (0-0.25 mm CR) and 5:0.67 (2.5-5 mm CR), respectively. To obtain the same particle size distribution (PSD) as the samples of ballast particles, RPB particles at different size fractions were weighed and mixed according to the PSD of the ballast sample, as shown in Figure 3.

Test Setup
A set of direct shear tests were performed with the large direct shear test rig as shown in Figure 2.
The tests were performed on ballast and RPB, respectively.
The DST rig is larger than the common ones, which can minimize the boundary effects sufficiently. As reported in [33,39], on the condition that specimen dimension is over eight times larger than particle size of the majority, the boundary effects can be ignored.
The DST rig consists of three main parts: two steel square boxes (shear boxes), two hydraulic jacks and two dial indicators (Figure 2d). The shear boxes consist of an upper steel square box with the dimension size (length × width × height) at 600 × 600 × 250 mm 3 , a lower steel square box (dimension size: 700 × 600 × 300 mm 3 ) and a steel loading plate (size: 600 × 600 × 20 mm 3 ). The steel wall thickness of the test rig is 20 mm. The maximum relative horizontal displacement of the two shear boxes is 100 mm, which is enough to reach the maximum shear stress.
Vertical and lateral hydraulic jacks can provide the maximum loading of 30 and 10 ton, respectively ( Figure 2d). The vertical jack actuator was used to apply a constant normal stress to the ballast assemblies, and the lateral jack actuator was applied to slowly move the lower shear box. A pressure sensor was attached next to the lateral jack actuator, which was used to measure the shear stress. The measuring range of the pressure sensor is 0-50 ton, and the resolution is ± 0.1%.  Vertical and lateral hydraulic jacks can provide the maximum loading of 30 and 10 ton, respectively ( Figure 2d). The vertical jack actuator was used to apply a constant normal stress to the ballast assemblies, and the lateral jack actuator was applied to slowly move the lower shear box. A pressure sensor was attached next to the lateral jack actuator, which was used to measure the shear stress. The measuring range of the pressure sensor is 0-50 ton, and the resolution is ± 0.1%.
The dial indicators were utilized to measure vertical and lateral displacements. The measuring range of the indicators is 0-30 mm, and the resolution is 0.001 mm. The measured vertical displacements were used to calculate the dilation of the samples, and the measured lateral displacement was used for the displacement-stress curve.

Test Procedure
The test procedure includes three steps. Firstly, two types of RPB particles were made with different CR sizes. Afterwards, one type of RPB particle weas placed in the DST rig by three layers. After filling in each layer, the assemblies were compacted with a compactor (Figure 2c). The compaction procedure was performed with a heavy steel weight, and the steel weight was dropped on RPB 50 times for each layer. The bulk density of the final samples was 1.59 g/cm3 (CR size 0-0.25 mm) and 1.42 g/cm3 (CR size 2.5-5 mm), respectively. The bulk density was calculated by the volume of the direct shear box and the weight of RPB particles. Finally, after the specimen was compacted The dial indicators were utilized to measure vertical and lateral displacements. The measuring range of the indicators is 0-30 mm, and the resolution is 0.001 mm. The measured vertical displacements were used to calculate the dilation of the samples, and the measured lateral displacement was used for the displacement-stress curve.

Test Procedure
The test procedure includes three steps. Firstly, two types of RPB particles were made with different CR sizes. Afterwards, one type of RPB particle weas placed in the DST rig by three layers. After filling in each layer, the assemblies were compacted with a compactor (Figure 2c). The compaction procedure was performed with a heavy steel weight, and the steel weight was dropped on RPB 50 times for each layer. The bulk density of the final samples was 1.59 g/cm3 (CR size 0-0.25 mm) and 1.42 g/cm3 (CR size 2.5-5 mm), respectively. The bulk density was calculated by the volume of the direct shear box and the weight of RPB particles. Finally, after the specimen was compacted with a flat ballast surface (for uniform vertical loading), the steel plate was placed on the top. Afterwards, the lower box was pushed until reaching 60 mm (10% of shear strain), and tests were under the normal stresses of 50, 100 and 200 kPa, respectively. The lower box was pushed with the speed rate of 1 cm/min, and a Sustainability 2020, 12, 2836 6 of 31 servo-controlled confining pressure was applied on the top steel plate. The same procedure at the second and third steps was performed on the other type of RPB particles.

RPB Particle Model
The basic elements to simulate ballast particles in PFC2D are discs. Due to the insufficient interlocking and unavoidable excessive rolling, only using discs is, in most cases, not accurate enough to present the natural characteristics of railway ballast, e.g., irregularity and angularity [33]. A solution has been proposed to model irregular particle shapes by the Clump or Cluster [31,[40][41][42], as shown in Figure 3. The Clump or Cluster is created by using two or more discs to present one particle. The difference between the Clump and the Cluster is whether the particle can break. The Clump is a rigid particle that cannot break no matter how large a force is applied on it. The Cluster is able to break because the component discs are bonded together by the parallel bonds. The Cluster breaks when the force applied on them is larger than the defined value [35,43]. The Clump or Clusters are better than using the discs to present ballast particles. Nevertheless, due to the computational costs, the Clump or Cluster normally cannot be applied in the large model or on conditions of numerous cyclic loading cycles. with a flat ballast surface (for uniform vertical loading), the steel plate was placed on the top. Afterwards, the lower box was pushed until reaching 60 mm (10% of shear strain), and tests were under the normal stresses of 50, 100 and 200 kPa, respectively. The lower box was pushed with the speed rate of 1 cm/min, and a servo-controlled confining pressure was applied on the top steel plate.
The same procedure at the second and third steps was performed on the other type of RPB particles.

RPB Particle Model
The basic elements to simulate ballast particles in PFC2D are discs. Due to the insufficient interlocking and unavoidable excessive rolling, only using discs is, in most cases, not accurate enough to present the natural characteristics of railway ballast, e.g., irregularity and angularity [33]. A solution has been proposed to model irregular particle shapes by the Clump or Cluster [31,[40][41][42], as shown in Figure 3. The Clump or Cluster is created by using two or more discs to present one particle. The difference between the Clump and the Cluster is whether the particle can break. The Clump is a rigid particle that cannot break no matter how large a force is applied on it. The Cluster is able to break because the component discs are bonded together by the parallel bonds. The Cluster breaks when the force applied on them is larger than the defined value [35,43]. The Clump or Clusters are better than using the discs to present ballast particles. Nevertheless, due to the computational costs, the Clump or Cluster normally cannot be applied in the large model or on conditions of numerous cyclic loading cycles.  To balance the simulation accuracy and computational costs, a simplified particle shape with a modified contact model (introduced in Section 2.3) was applied in our models. This method has been validated and effectively applied in many studies, e.g., [45,46]. As shown in Figure 4, the ballast particle was modelled with two discs, and RPB was modelled by bonding small discs to the ballast particle with the parallel bonds. Particle size distributions (PSDs) are shown in the figure as well, and the PSDs in the tests and models were accordant and meet the British standard [38]. The PSD in the model was obtained by generating RPB particles in different size fractions to the required mass percentages, and the required mass percentages were the same as the particle size distribution of the experimental DST.
In order to obtain the required PSD, two RPB templates were created for each size fraction (eight templates in total). The creation process includes three steps. Firstly, the ballast particle (made by the large disc and medium disc) was designed to make its size within the size fraction (e.g. 31.5-40 mm). Specifically, the large disc and medium disc sizes of the eight templates particles are given in Table  2. Afterwards, based on the perimeter of the large and medium discs, the small discs that were used  To balance the simulation accuracy and computational costs, a simplified particle shape with a modified contact model (introduced in Section 2.3) was applied in our models. This method has been validated and effectively applied in many studies, e.g., [45,46]. As shown in Figure 4, the ballast particle was modelled with two discs, and RPB was modelled by bonding small discs to the ballast particle with the parallel bonds. Particle size distributions (PSDs) are shown in the figure as well, and the PSDs in the tests and models were accordant and meet the British standard [38]. The PSD in the model was obtained by generating RPB particles in different size fractions to the required mass percentages, and the required mass percentages were the same as the particle size distribution of the experimental DST.
In order to obtain the required PSD, two RPB templates were created for each size fraction (eight templates in total). The creation process includes three steps. Firstly, the ballast particle (made by the large disc and medium disc) was designed to make its size within the size fraction (e.g. 31.5-40 mm). Specifically, the large disc and medium disc sizes of the eight templates particles are given in Table 2. Afterwards, based on the perimeter of the large and medium discs, the small discs that were used to simulate CR chips were determined at the disc numbers and positions, as shown in Table 2. The positions were confirmed through the relative coordinates to the ballast particle (large and medium discs). Finally, according to the diameters and relative coordinates, RPB templates were created. . Direct shear test model, particle size distribution and one RPB particle.
The CR particles were selected as 4 mm, which is in the range of 2.5-5 mm. Afterwards, the parameters for the model were confirmed by comparing the numerical simulation results with the experimental ones. Finally, another DST model was developed with RPB made by 2 mm CR particles. The DST results of different CR size RPB (2 or 4 mm) were compared. It needs to note that RPBs made by the 0-0.25 mm CR particles are not modelled in this study because it is nearly impossible to model that due to the huge computational costs. When the CR size is quite small (almost powder), then many discs are needed to cover the modelled ballast particle (two-disc). This kind of RPB particle is not applicable for the large-scale DST model (four times larger by volume than the normal DST model), not to mention the three-sleeper full track model ( Figure 6).

Direct Shear Test Model
The rolling resistance contact model (with simple particles) was applied in the DST and threesleeper track models. The rolling resistance is applied by adding rolling friction at contacts between modelled ballast particles, as shown in Figure 5. Compared with the widely-used linear contact model (in almost all earlier studies), it is better at providing a realistic performance of ballast assemblies by restricting relative particle rotation as proved in [45].
In the PFC, the rolling resistance contact model was developed by improving the linear contact model. In other words, it was created by adding a new algorithm to the linear contact model. This applies a turning moment to the contact area to resist relative rotation. It has one more parameter (i.e., rolling friction) compared with the linear contact model. Specifically, the rolling friction is used to resist the particle rotation. The maximum rotation restriction is equal to the product of the rolling friction with the corresponding normal force. The restriction effect can be regarded as the rolling stiffness, which is similar to the clockwork spring ( Figure 4). The parameters of the rolling resistance contact model are given in Table 3.  . Direct shear test model, particle size distribution and one RPB particle. It is worth noting that in the PFC2D, the "Generate" command can generate the Clump using the RPB templates according to the user-designated size fraction to produce the required PSD. However, RPB is the Cluster, which cannot be generated using the "Generate" command according to the templates. In addition, the density of the component discs in one Clump is the same value, and the other material characteristics in one Clump (e.g., shear modulus) are also the same value. Because the ballast and crumb rubber in RPB have different material characteristics, the Clump is not suitable to simulate RPB particles. Therefore, the model is built with Clumps initially, and afterwards the Clumps are replaced by the Clusters. The Cluster can present RPB, because the discs in one Cluster have various parameters to present different characteristics for different materials. The replacement process of Clumps to Clusters is named as "Particle-replacing".
The specific process of the Particle-replacing is: after the model was built with the Clumps (the model creation process is explained in Section 2.3), the coordinates and diameters of discs in every Clump were obtained. Afterwards, the Clump was deleted, and according to the coordinates and diameters of discs in the clump, a new particle was created at the same position with same discs, but the new particle was a Cluster. The Cluster was made by two parts, the ballast particle (two overlapped discs, uncrushable) and the CR chips (bonded to the ballast particle with parallel bonds).
Particularly, the parallel bonds gives the physical performance of a cement-like substance sticking together the two contacting particles [47]. When a force is acted on a parallel bonded particle, the particle develops a force and moment within the bond due to a relative motion between the corresponding two spheres. When the force applied on the particle exceeds the bond strength, the parallel bonds are removed, together with the corresponding force and moment [47]. The effectiveness of the parallel bond is quite similar to the polyurethane. The parameters of the modelled particle (ballast and CR) are given in Table 3. The parameters of the parallel bond (to simulate the RPB binders) were decided according to the studies on the polyurethane [48][49][50]. In these studies, the glued ballast particles were built with parallel bonds to simulate the glued ballast particles with the polyurethane, and the parallel bond parameters of the polyurethane were used in building RPB particles, as shown in Table 3.
The CR particles were selected as 4 mm, which is in the range of 2.5-5 mm. Afterwards, the parameters for the model were confirmed by comparing the numerical simulation results with the experimental ones. Finally, another DST model was developed with RPB made by 2 mm CR particles. The DST results of different CR size RPB (2 or 4 mm) were compared. It needs to note that RPBs made by the 0-0.25 mm CR particles are not modelled in this study because it is nearly impossible to model that due to the huge computational costs. When the CR size is quite small (almost powder), then many discs are needed to cover the modelled ballast particle (two-disc). This kind of RPB particle is not applicable for the large-scale DST model (four times larger by volume than the normal DST model), not to mention the three-sleeper full track model.

Direct Shear Test Model
The rolling resistance contact model (with simple particles) was applied in the DST and three-sleeper track models. The rolling resistance is applied by adding rolling friction at contacts between modelled ballast particles, as shown in Figure 5. Compared with the widely-used linear contact model (in almost all earlier studies), it is better at providing a realistic performance of ballast assemblies by restricting relative particle rotation as proved in [45].  . Direct shear test model, particle size distribution and one RPB particle.
The CR particles were selected as 4 mm, which is in the range of 2.5-5 mm. Afterwards, the parameters for the model were confirmed by comparing the numerical simulation results with the experimental ones. Finally, another DST model was developed with RPB made by 2 mm CR particles. The DST results of different CR size RPB (2 or 4 mm) were compared. It needs to note that RPBs made by the 0-0.25 mm CR particles are not modelled in this study because it is nearly impossible to model that due to the huge computational costs. When the CR size is quite small (almost powder), then many discs are needed to cover the modelled ballast particle (two-disc). This kind of RPB particle is not applicable for the large-scale DST model (four times larger by volume than the normal DST model), not to mention the three-sleeper full track model ( Figure 6).

Direct Shear Test Model
The rolling resistance contact model (with simple particles) was applied in the DST and threesleeper track models. The rolling resistance is applied by adding rolling friction at contacts between modelled ballast particles, as shown in Figure 5. Compared with the widely-used linear contact model (in almost all earlier studies), it is better at providing a realistic performance of ballast assemblies by restricting relative particle rotation as proved in [45].
In the PFC, the rolling resistance contact model was developed by improving the linear contact model. In other words, it was created by adding a new algorithm to the linear contact model. This applies a turning moment to the contact area to resist relative rotation. It has one more parameter (i.e., rolling friction) compared with the linear contact model. Specifically, the rolling friction is used to resist the particle rotation. The maximum rotation restriction is equal to the product of the rolling friction with the corresponding normal force. The restriction effect can be regarded as the rolling stiffness, which is similar to the clockwork spring ( Figure 4). The parameters of the rolling resistance contact model are given in Table 3.   In the PFC, the rolling resistance contact model was developed by improving the linear contact model. In other words, it was created by adding a new algorithm to the linear contact model. This applies a turning moment to the contact area to resist relative rotation. It has one more parameter (i.e., rolling friction) compared with the linear contact model. Specifically, the rolling friction is used to resist the particle rotation. The maximum rotation restriction is equal to the product of the rolling friction with the corresponding normal force. The restriction effect can be regarded as the rolling stiffness, which is similar to the clockwork spring ( Figure 4). The parameters of the rolling resistance contact model are given in Table 3.
The DST model is shown in Figure 4, and the dimension size is based on the original test rig, which is 700 × 250 mm2 (lower box) and 600 × 300 mm2 (upper box). The procedure of building the DST model with the specimen has three steps.
1. A taller container (than the DST test rig) was built for containing the Clumps, which were generated according to the PSD of the experimental DST. While the Clumps were generating, the RPB templates were used (eight templates in total, Section 2.2) and the total area of the Clumps was calculated by the desired porosity and the modelled DST rig. The "Generate" command generates Clumps that do not have any overlaps between each other; because of this, a bigger container is needed; 2. As the Clumps were randomly generated in the bigger container, the Clumps needed to drop to the bottom and settle. For this, the frictions (including rolling friction and translation friction) were set to 0.0 and the gravity was set at 9.81 m/s2 to make the Clumps fall down. Afterwards, all the Clumps dropped to the container bottom, and hundreds of cycles were performed until the specimen was eventually settled. Finally, the bigger container was deleted and the DST rig was generated with the rigid walls, which were also a basic element in PFC2D. The specimen porosity was 0.1, which was lower than the experimental one (0.42). This is due to the fact that the 2D model has much smaller voids than models in 3D, and the 2D samples are easier to compact. This is a normal phenomenon in most of the 2D DEM models [15,31,40]. The Clumps that were out of the DST rig were deleted. Moreover, the frictions were set to the normal values, as shown in Table 3; 3. Hundreds of cycles were performed until the specimen is settled. After settling, RPB particles out of the test rig are deleted. Afterwards, the "Particle-replacing" process was performed for the settled specimen (introduced in Section 2.2). After the DST model was built, a further shearing process was performed. The normal stress was applied on the specimen with the servo control mechanism. After the specimen was settled, one of the normal stresses (50 kPa) was applied upon the sample. Afterwards, the lower box slowly moved at 5 mm/s. The other two tests under normal stresses (100,200 kPa) were performed, subsequently. The DST model with 2 mm CR RPB was performed using the same model procedure as above. During the simulated tests, lower box displacements were recorded, as well as the shear stress, contact force and particle displacement. Figure 6a shows the three-sleeper track model with dimensions at 2100 × 500 mm 2 , and the dimension of sleepers is at 250 × 200 mm 2 . The three-sleeper track model was built using the ballast particles (two-disc Clumps) at the beginning, afterwards, according to different simulation conditions (i.e., RPB thickness) the ballast particles under the sleeper were replaced by RPB particles (Clusters, introduced in Section 2.2). Finally, the cyclic loadings were applied to ballast-RPB layer by the three sleepers. Figure 6a shows the three-sleeper track model with dimensions at 2100 × 500 mm 2 , and the dimension of sleepers is at 250 × 200 mm 2 . The three-sleeper track model was built using the ballast particles (two-disc Clumps) at the beginning, afterwards, according to different simulation conditions (i.e., RPB thickness) the ballast particles under the sleeper were replaced by RPB particles (Clusters, introduced in Section 2.2). Finally, the cyclic loadings were applied to ballast-RPB layer by the three sleepers.  The three-sleeper track model was firstly built with the ballast particles (two-disc Clumps). The ballast particles were generated in a taller container and then dropped to the bottom as the gravity, The three-sleeper track model was firstly built with the ballast particles (two-disc Clumps). The ballast particles were generated in a taller container and then dropped to the bottom as the gravity, which is the same as first procedure, as introduced in Section 2.3. Afterwards, the ballast particles that were at the positions of sleepers were deleted, and the sleepers were generated at those positions. After that, hundreds of cycles were performed to settle the ballast bed.

Three-sleeper Track Model
Afterwards, three steps of changing were performed on the model as follows. Firstly, part of the ballast particles under the sleeper was replaced by RPB particles, which made a thickness of two layers together at 30 cm. The two layers are the RPB layer (under sleeper) and ballast layer (under RPB). The same RPB particles (Figure 4) whose parameters were confirmed (by DSTs) were applied in the three-sleeper track models using the same PSD. In order to obtain the optimal RPB thickness, six models were built including different size CR RPBs (2 or 4 mm) with corresponding different RPB thicknesses (20, 25 or 30 cm).
Secondly, the cyclic loadings were the 90 • out-of-phase loading, which is the same as the study in [28], as shown in Figure 6b. Specifically, the forces applied to the three sleepers were based on the equation below. In the equation, A is the amplitude, 30 kN; f is the frequency, 5 Hz; T is the period, 0.2 s. The final sinusoidal loading is in the range of 5 kN to 65 kN. A total of 100 cycles are simultaneously applied to each sleeper.
Finally, the necessary results were obtained, including the settlement, contact force, particle displacement and acceleration. The settlement was obtained by recording the sleeper positions, and the contact force and displacement were shown by the software PFC. It should be noted that the particle accelerations were recorded at certain positions (little differences between the six models) during the cyclic loadings, as shown in Figure 6a. The recorded particles (in red color) can be regarded as the acceleration sensors, measuring the accelerations of the whole loading process.

Model Parameters
After calibration, the parameters of the ballast, CR and polyurethane applied in the following models are given in Table 3. In the table, the PB is short for the parallel bond, which was used to simulate the polyurethane. The rubber parameters are from the reference [51], and the polyurethane parameters are from the reference [49]. The test rig parameters and the ballast parameters were confirmed by comparing the DST simulation results with the experimental results, which are explained in detail in the next section.
The parameters in Table 3 were confirmed according to the earlier studies. For example, in [1], the normal and shear stiffnesses were set as the 2e6 and 1e6, respectively. In [31], the normal and shear stiffnesses were 2.5e8 and 2e8, respectively. In [52], the normal and shear stiffnesses were 4.82e8 and 2.41e8, respectively. In [53], the normal and shear stiffnesses were 5.2e7 and 5.2e7, respectively. From the above values, it can be seen that the parameters are within a range. According to the direct shear test results, the proper values that can match the test results and the simulation results were chosen.

Displacement-Stress and Dilation
As shown in Figure 7, the parameters used for the ballast and test rig are given in Table 3 based on the illustrated results. It can be seen that the simulation and test results can acceptably be matched. Specifically, for the tests and simulations of ballast particles, the result differences are within 10% (displacement-shear stress) and 7.5% (dilation). This is acceptable for the following simulations.     Particularly, in the experimental tests, the CR size significantly influences RPB resilience and also the interaction between RPB particles, consequently, the dilation of RPB has a large range of variation (Figure 7d). From Figure 7a,c,e, it can be observed that RPB particles provide 1/3~1/6 shear stress as the ballast particles. In addition, smaller CR RPB particles can provide a higher shear stress, but still lower than the shear stress of ballast particles.
The experimental tests can provide a macro performance of RPB, and the simulation is able to show the mesoscopic performance, which is given in the following sections. Figure 8 presents the contact force results of the condition that is 60 mm shear displacement under the 100 kPa normal stress. The other results of contact force (under 50 or 200 kPa normal stress) can be found in Figure A1. Particularly, in the experimental tests, the CR size significantly influences RPB resilience and also the interaction between RPB particles, consequently, the dilation of RPB has a large range of variation (Figure 7d). From Figure 7a/c/e, it can be observed that RPB particles provide 1/3~1/6 shear stress as the ballast particles. In addition, smaller CR RPB particles can provide a higher shear stress, but still lower than the shear stress of ballast particles.

Contact Force
The experimental tests can provide a macro performance of RPB, and the simulation is able to show the mesoscopic performance, which is given in the following sections. From the figure, it can be observed that the biggest contact force reduces from 27.1 kN (ballast) to 18.6 kN (2mm CR RPB) or 13.9 kN (4mm CR RPB). In addition, for the 4mm CR RPB, the force distribution is more homogeneous and the force chain is not obvious during RPB shear test. This is due to the fact that the CR can soften the contacts by increasing the contact numbers and areas. Figure 9 presents the particle displacements of the ballast or RPB with the shear displacement at 60 mm and under the normal stress 100 kPa. The other simulation results of the particle displacements are given in Table A2. From the figure, it can be seen that the CR size can influence the particle translation direction. To be more specific, the 4 mm CR RPB particles have lower value of particle displacements (in upper box) and the direction is approximately from left to right. This is quite different from the ballast and 2 mm CR RPB particle translation direction, which is in the up direction. Moreover, it can be observed that the 2 mm CR RPB has the more large-displacement particles at the shearing surface, which means it can transmit the forces not only to the up direction, but also to the down direction. Moreover, from Figure 9c, it can be seen that the particle displacements are From the figure, it can be observed that the biggest contact force reduces from 27.1 kN (ballast) to 18.6 kN (2 mm CR RPB) or 13.9 kN (4 mm CR RPB). In addition, for the 4mm CR RPB, the force distribution is more homogeneous and the force chain is not obvious during RPB shear test. This is due to the fact that the CR can soften the contacts by increasing the contact numbers and areas. Figure 9 presents the particle displacements of the ballast or RPB with the shear displacement at 60 mm and under the normal stress 100 kPa. The other simulation results of the particle displacements are given in Figure A2. From the figure, it can be seen that the CR size can influence the particle translation direction. To be more specific, the 4 mm CR RPB particles have lower value of particle displacements (in upper box) and the direction is approximately from left to right. This is quite different from the ballast and 2 mm CR RPB particle translation direction, which is in the up direction. Moreover, it can be observed that the 2 mm CR RPB has the more large-displacement particles at the shearing surface, which means it can transmit the forces not only to the up direction, but also to the down direction. Moreover, from Figure 9c, it can be seen that the particle displacements are horizontal instead of going upwards, which is a reason that the dilation of 4 mm CR RPB has a negative value (Figure 7d).

Three-sleeper track model results
The three-sleeper track model results include the settlement, the contact force, particle displacement and particle acceleration. The settlement is obtained from the sleeper positions after 100 cyclic loading cycles. The contact forces and the displacements of all the particles are shown at the 50th cyclic loading cycle, while the particle accelerations are recorded during the 100 cycles.

Settlement
The settlements of the three-sleeper track models with ballast or ballast and RPB are shown in Figure 10. The figure presents only parts of the results; more results can be found in Table Figure A.3 (Appendix). Figure 10a presents the applied force-settlement curves of the middle sleepers, and it shows the settlements (100 cycles) of three conditions: 1) only ballast; 2) only 2 mm CR RPB under the sleeper; and 3) only 4 mm CR RPB under the sleeper. From Figure 10a, it can be observed that the settlement increases the CR size, which means the 4 mm CR RPB has the biggest settlement (around 70 mm) after cyclic loadings. Moreover, the 2 mm CR RPB has a smaller settlement (around 27 mm), which is still larger than the only ballast (around 4 mm).
From Figure 10b, it can be seen that the track stiffness reduces a lot when replacing the ballast under the sleeper with RPB, and the 4 mm CR RPB with thickness at 30 cm has the lowest stiffness. However, it has a lower settlement than 4 mm CR RPB with thickness at 25 cm (Figure 10c). Figure  10c also presents that 2 mm CR RPB (27, 40, 53 mm) has less settlement than the 4 mm CR RPB (48, 70, 85 mm), and the 2 mm CR RPB with a thickness at 30 cm can be the optimal choice. Surprisingly, the settlements for the 2 mm CR RPB are not as imagined. The minimum settlement is RPB thickness at 30 cm, and the maximum settlement is the 20 cm thickness. This suggest that the settlement value is reduced when increasing RPB layer thickness. However, this is not observed in the 4 mm CR RPB. Two reasons can be considered for this phenomenon. On one hand, this may result from the random RPB movements, because the interlockings between 4 mm CR RPB particles are not strong. This will lead to contacts between sleeper and RPB particles being insufficient. For this, higher

Three-sleeper Track Model Results
The three-sleeper track model results include the settlement, the contact force, particle displacement and particle acceleration. The settlement is obtained from the sleeper positions after 100 cyclic loading cycles. The contact forces and the displacements of all the particles are shown at the 50th cyclic loading cycle, while the particle accelerations are recorded during the 100 cycles.

Settlement
The settlements of the three-sleeper track models with ballast or ballast and RPB are shown in Figure 10. The figure presents only parts of the results; more results can be found in Figure A3. Figure 10a presents the applied force-settlement curves of the middle sleepers, and it shows the settlements (100 cycles) of three conditions: (1) only ballast; (2) only 2 mm CR RPB under the sleeper; and (3) only 4 mm CR RPB under the sleeper. From Figure 10a, it can be observed that the settlement increases the CR size, which means the 4 mm CR RPB has the biggest settlement (around 70 mm) after cyclic loadings. Moreover, the 2 mm CR RPB has a smaller settlement (around 27 mm), which is still larger than the only ballast (around 4 mm).
From Figure 10b, it can be seen that the track stiffness reduces a lot when replacing the ballast under the sleeper with RPB, and the 4 mm CR RPB with thickness at 30 cm has the lowest stiffness. However, it has a lower settlement than 4 mm CR RPB with thickness at 25 cm (Figure 10c). Figure 10c also presents that 2 mm CR RPB (27, 40, 53 mm) has less settlement than the 4 mm CR RPB (48, 70, 85 mm), and the 2 mm CR RPB with a thickness at 30 cm can be the optimal choice. Surprisingly, the settlements for the 2 mm CR RPB are not as imagined. The minimum settlement is RPB thickness at 30 cm, and the maximum settlement is the 20 cm thickness. This suggest that the settlement value is reduced when increasing RPB layer thickness. However, this is not observed in the Sustainability 2020, 12, 2836 15 of 31 4 mm CR RPB. Two reasons can be considered for this phenomenon. On one hand, this may result from the random RPB movements, because the interlockings between 4 mm CR RPB particles are not strong. This will lead to contacts between sleeper and RPB particles being insufficient. For this, higher impact loadings are randomly applied on the RPB layer, and eventually this causes the uncertain relationship between settlement and RPB thickness. On the other hand, the initial compaction of the RPB bed is not easy to keep the same, and RPB beds have differences. In addition, the interactions between RPB particles are not strong, which will result in some RPB particles with randomly higher accelerations and displacements. This can also be reflected in Figure 10b, which presents that the 2 mm CR RPB displacement of 25 cm thickness is smaller than that of 20 cm thickness. Afterwards, after 100 cycles, the settlement of 25 cm thickness (2 mm CR RPB) is smaller than that of 20 cm thickness. From this, it can be also be concluded that the effect factors for ballast/RPB studies are a lot and it is nearly impossible to consider all of them, e.g., initial condition, bulk density, particle rearrangement and particle movements. How they influence the macro performance of the ballast/RPB bed needs deeper understanding, especially it is necessary to understand how they influence each other and their evolution processes.
The above-mentioned two reasons can influence the results, but the comparison among the ballast, the 2 mm CR RPB and 4 mm CR RPB are not influenced. For example, even though RPB influences the interlockings between RPB particles, the 2 mm CR RPB particles have stronger interlockings than the 4 mm CR RPB. Therefore, it can be seen that the ballast has the least settlement, and the 4 mm CR RPB has the largest settlement. The smaller CR RPB has better performance at the settlement. Due to several effect factors influencing the settlement results, the following mesoscopic performance analysis is presented to confirm the optimal RPB thickness and CR size. On the other hand, the initial compaction of the RPB bed is not easy to keep the same, and RPB beds have differences. In addition, the interactions between RPB particles are not strong, which will result in some RPB particles with randomly higher accelerations and displacements. This can also be reflected in Figure 10b, which presents that the 2 mm CR RPB displacement of 25 cm thickness is smaller than that of 20 cm thickness. Afterwards, after 100 cycles, the settlement of 25 cm thickness (2 mm CR RPB) is smaller than that of 20 cm thickness. From this, it can be also be concluded that the effect factors for ballast/RPB studies are a lot and it is nearly impossible to consider all of them, e.g., initial condition, bulk density, particle rearrangement and particle movements. How they influence the macro performance of the ballast/RPB bed needs deeper understanding, especially it is necessary to understand how they influence each other and their evolution processes.

Contact Force
The above-mentioned two reasons can influence the results, but the comparison among the ballast, the 2 mm CR RPB and 4 mm CR RPB are not influenced. For example, even though RPB influences the interlockings between RPB particles, the 2 mm CR RPB particles have stronger interlockings than the 4 mm CR RPB. Therefore, it can be seen that the ballast has the least settlement, and the 4 mm CR RPB has the largest settlement. The smaller CR RPB has better performance at the settlement. Due to several effect factors influencing the settlement results, the following mesoscopic performance analysis is presented to confirm the optimal RPB thickness and CR size. Figure 11 presents the contact force after 50 loading cycles, and three types of tracks are shown, including 1) only ballast (30 cm thickness), 2) 2 mm CR RPB (30 cm thickness) and 3) 4 mm CR RPB (30 cm thickness). Other results can be found in Figure A4. settlement. Due to several effect factors influencing the settlement results, the following mesoscopic performance analysis is presented to confirm the optimal RPB thickness and CR size. Figure 11 presents the contact force after 50 loading cycles, and three types of tracks are shown, including 1) only ballast (30 cm thickness), 2) 2 mm CR RPB (30 cm thickness) and 3) 4 mm CR RPB (30 cm thickness). Other results can be found in Table Figure  From Figure 11, it can be seen that the contact force distribution of RPB is bigger than that of the only ballast. This can be reflected by the distribution angles, which are 58° (ballast), 42° (2 mm CR RPB) and 45° (4 mm CR RPB). The angle is based on the direction of the large contact forces. Particularly, for the 2 mm CR RPB, the contact forces are more homogeneous. This is due to its contact areas and numbers hibeinggher, which can also be observed in the DST contact force results ( Figure  8). Thus, a smaller CR size is recommended for RPB application.

Contact Force
Moreover, the maximum contact forces of RPB (28.0, 16.1 kN) are higher than that of the ballast (15.4 kN). This is due to the fact that the CR can induce insufficiently soft contacts between the sleeper and RPB particles, which has been proved in [54], where the soft contacts led to higher sleeper accelerations. However, in some earlier studies on the under sleeper pads (USPs) [21,22,55], they argued that the soft contacts can provide a better ballast bed performance by reducing the ballast degradation. Particularly, the difference between the USPs and RPB is that the USPs attach to the sleeper bottom without any movement and the ballast particles' rearrangement is slow, whereas RPB particles can move randomly after applied loadings and the interaction between particles is not strong enough to restrict RPB particles. Because of this, after cyclic loadings the contacts remain From Figure 11, it can be seen that the contact force distribution of RPB is bigger than that of the only ballast. This can be reflected by the distribution angles, which are 58 • (ballast), 42 • (2 mm CR RPB) and 45 • (4 mm CR RPB). The angle is based on the direction of the large contact forces. Particularly, for the 2 mm CR RPB, the contact forces are more homogeneous. This is due to its contact areas and numbers hibeinggher, which can also be observed in the DST contact force results (Figure 8). Thus, a smaller CR size is recommended for RPB application.
Moreover, the maximum contact forces of RPB (28.0, 16.1 kN) are higher than that of the ballast (15.4 kN). This is due to the fact that the CR can induce insufficiently soft contacts between the sleeper and RPB particles, which has been proved in [54], where the soft contacts led to higher sleeper accelerations. However, in some earlier studies on the under sleeper pads (USPs) [21,22,55], they argued that the soft contacts can provide a better ballast bed performance by reducing the ballast degradation. Particularly, the difference between the USPs and RPB is that the USPs attach to the sleeper bottom without any movement and the ballast particles' rearrangement is slow, whereas RPB particles can move randomly after applied loadings and the interaction between particles is not strong enough to restrict RPB particles. Because of this, after cyclic loadings the contacts remain almost the same for the USPs, however, for RPB, the contacts are random and stress concentration will possibly happen. Therefore, it is not recommended that RPB particles are directly placed under the sleeper. Figure 12 presents the particle displacement after 50 cyclic loadings of the ballast and RPB (2 mm or 4 mm CR). The other displacement results of different RPB thicknesses can be found in Figure A5. From the figure, it can be observed that the maximum displacement of RPB (10.3, 11.3 cm) is around six times larger than the ballast (1.7 cm). The maximum displacements of RPB happen under the sleeper, and the displacement of the ballast under the sleeper is round 1.1 cm. This is due to the RPB being soft and the contacts being weakened, causing the large displacements. Therefore, replacing the whole ballast layer into the RPB layer is not recommended. almost the same for the USPs, however, for RPB, the contacts are random and stress concentration will possibly happen. Therefore, it is not recommended that RPB particles are directly placed under the sleeper. Figure 12 presents the particle displacement after 50 cyclic loadings of the ballast and RPB (2 mm or 4 mm CR). The other displacement results of different RPB thicknesses can be found in Table  Figure A.5 (Appendix). From the figure, it can be observed that the maximum displacement of RPB (10.3, 11.3 cm) is around six times larger than the ballast (1.7 cm). The maximum displacements of RPB happen under the sleeper, and the displacement of the ballast under the sleeper is round 1.1 cm. This is due to the RPB being soft and the contacts being weakened, causing the large displacements. Therefore, replacing the whole ballast layer into the RPB layer is not recommended.

Acceleration
The particle acceleration at the position beneath the sleeper (15 cm; Position 5) is shown in Figure  13 as an example. The positions of the acceleration-measured particles can be found in Figure 6a. All the acceleration results can be found in Table Figure A.6 and Table A7 (Appendix). The ballast acceleration results are compared with the experimental results in [56] and they match well, which demonstrates that the numerical acceleration results are reliable.

Acceleration
The particle acceleration at the position beneath the sleeper (15 cm; Position 5) is shown in Figure 13 as an example. The positions of the acceleration-measured particles can be found in Figure 6a. All the acceleration results can be found in Figure A6 and Supplementary Materials. The ballast acceleration results are compared with the experimental results in [56] and they match well, which demonstrates that the numerical acceleration results are reliable.

Acceleration
The particle acceleration at the position beneath the sleeper (15 cm; Position 5) is shown in Figure  13 as an example. The positions of the acceleration-measured particles can be found in Figure 6a. All the acceleration results can be found in Table Figure A.6 and Table A7 (Appendix). The ballast acceleration results are compared with the experimental results in [56] and they match well, which demonstrates that the numerical acceleration results are reliable. From the figure, it can be seen that the acceleration of RPB is much higher than the ballast, and large acceleration values can be found frequently for RPB. In addition, at the beginning of the cyclic loadings, a peak acceleration value can be observed for RPB. This demonstrates that replacing the ballast layer with an RPB layer should be carefully tested before being applied in the field because of the high acceleration. The acceleration is due to the fact that the high RPB particle resilience results in the hanging sleeper and insufficient contacts (low interlocking) between RPB particles. Therefore, the RPB and ballast mixture can be a solution for both ballast protection and sufficient performance.
Moreover, the 2 mm CR RPB with 30 cm thickness has the lowest and most stable acceleration, as shown in Figure 13a. However, for the 4 mm CR RPB, the layer thickness seems to have few influences on the acceleration. This demonstrates that a large CR size leads to the RPB layer having a more complex situation that is not easy to control. Moreover, the smaller CR size RPB is recommended to be used. From the figure, it can be seen that the acceleration of RPB is much higher than the ballast, and large acceleration values can be found frequently for RPB. In addition, at the beginning of the cyclic loadings, a peak acceleration value can be observed for RPB. This demonstrates that replacing the ballast layer with an RPB layer should be carefully tested before being applied in the field because of the high acceleration. The acceleration is due to the fact that the high RPB particle resilience results in the hanging sleeper and insufficient contacts (low interlocking) between RPB particles. Therefore, the RPB and ballast mixture can be a solution for both ballast protection and sufficient performance.

Conclusions and perspectives
Moreover, the 2 mm CR RPB with 30 cm thickness has the lowest and most stable acceleration, as shown in Figure 13a. However, for the 4 mm CR RPB, the layer thickness seems to have few influences on the acceleration. This demonstrates that a large CR size leads to the RPB layer having a more complex situation that is not easy to control. Moreover, the smaller CR size RPB is recommended to be used.

Conclusions and Perspectives
In this paper, the performance (static and dynamic) of RPB is compared with that of the ballast with the DEM models, i.e., direct shear test model and three-sleeper track model. The direct shear test model is validated by the experimental tests, afterwards, the parameters are confirmed and applied in the three-sleeper track model. In the models, the contact force, displacement and acceleration results are measured and studied. From the results, the following conclusions can be made.

1.
RPB with a smaller crumb rubber size can provide higher shear strength, but still lower than the ballast. However, RPB can provide more homogeneous contact forces; 2.
Specifically, during the direct shear tests, the contact forces between RPB are more homogeneous than normal ballast; 3.
The RPB with smaller crumb rubber size has less settlement and for the 2 mm crumb rubber size RPB, the 30 cm thickness is recommended because it has the least settlement; 4.
An RPB of bigger crumb rubber size can soften the layer, and the smaller crumb rubber size RPB can provide better particle interlockings; 5.
The force distribution of RPB is better than the ballast, due to the soft contacts leading to more contact areas and numbers. However, the soft contacts also induce large particle displacements and accelerations. Therefore, completely replacing the ballast layer with RPB is not recommended, and an RPB and ballast mixture can be a better solution.

Perspective
Replacing the whole ballast layer with an RPB is studied in this paper, however, mixing RPB with ballast has not been studied to date. The experimental tests of a cyclic test on RPB should be performed to validate the numerical results. Moreover, the cyclic loading cycles are 100 in this study, which can only present short-term settlement performance. The long-term settlement performance should also be checked. It needs to be noted that the 2D model is still not adequate, and further studies will focus on building 3D models. Consequently, further studies will be performed at these directions.

Statement
The aim of this paper is only for the demonstration of the DEM method, and RPB is not for any commercial purposes, only for the research aim. The technique or process of producing RPB is similar with the Neoballast, however, different materials (ballast, binder and rubber) are used. The study is only an example of showing how to use the DEM for rubber-ballast related studies.    Figure A3. Applied force-settlement curve. Figure A2. Displacements of ballast or RPB at 60 mm shear displacement.       Figure A6. Particle acceleration during the cyclic loadings.