Discrete element modelling of railway ballast performance considering particle shape and rolling resistance

To simulate ballast performance accurately and efficiently, the input in discrete element models should be carefully selected, including the contact model and applied particle shape. To study the effects of the contact model and applied particle shape on the ballast performance (shear strength and deformation), the direct shear test (DST) model and the large-scale process simulation test (LPST) model were developed on the basis of two types of contact models, namely the rolling resistance linear (RRL) model and the linear contact (LC) model. Particle shapes are differentiated by clumps. A clump is a sphere assembly for one ballast particle. The results show that compared with the typical LC model, the RRL method is more efficient and realistic to predict shear strength results of ballast assemblies in DSTs. In addition, the RRL contact model can also provide accurate vertical and lateral ballast deformation under the cyclic loading in LPSTs.


Introduction
Railways play a significant role in the transportation system worldwide and work in many sectors (urban rail, highspeed railway, heavy haul, intercity and metro) [1,2]. Ballasted tracks, as the most widely used track type, consist of rails, sleepers and the ballast layer [3,4]. It possesses the advantages such as low construction cost, simple design and construction, and easy maintenance [5].
The ballast layer, a crucial component of ballasted track, provides resistances to sleepers, transmits and distributes the loads or impacts from sleepers to the subgrade, as well as allows rapid drainage [6]. Generally, it is composed of blasted (quarried) rock aggregate, which is required to meet certain characteristics such as narrow-graded (20-60 mm) and irregular particle shape, specific surface roughness, density, hardness, resistance to attrition and weathering [7]. Even though various railway ballast standards in terms of particle size distribution or particle shape have already been formulated [7][8][9], their influences on ballast performance (resilience, shearing strength, and settlement) have not been sufficiently studied [10,11].
Laboratory or field tests are of limited use in studying the ballast performance, because the test conditions cannot be kept the same and many characteristics (e.g. ballast density and sleeper type) cannot be controlled [12]. Additionally, due to the discrete nature of ballast, it is not accurate or realistic to use the finite element method, which simulates the ballast layer as continuous layer [13]. The ballast performance keeps changing due to the ballast degradation (abrasion and breakage) [12,[14][15][16]. In addition, the sliding and rolling of individual ballast particles also influence the performance of the ballast layer [17].
On the one hand, the computational cost is the most considerable limitation in developing DEM models that may have millions of spheres (e.g. full-scale track model) [21]. Larger number of particles means the increase in the total number of particle contacts, which results in great computational cost. This problem becomes more severe when non-spherical particles are present in the DEM models. The usage of the non-spherical particles can provide more realistic load-deformation response [18,33]. A non-spherical particle is generally made by a sphere assembly, named clump or cluster in the particle flow code (PFC, commercial DEM software) [27,34]. Using the nonspherical particles (sphere assembly) increases the spheres and the number of contacts (contact points between the particles). The contacts are updated with every cycle according to the force-displacement law, which finally increases the computation time considerably.
In most cases, the contact method used in the earlier models was elementary linear model (spring-damping model). By using the RRL, simple spheres can also be possible to attain similar ballast performance, which can save a large amount of computational time. For example, it was demonstrated that the linear rolling resistance contact model (using spheres as ballast particles) can obtain the same ballast lateral resistance results as those from field tests [35].
On the other hand, even if the sleeper-ballast model uses the simple spherical or less-spherical particles, it has to be developed in a large-scale manner (e.g. three-sleeper track model) due to the scale effect and the boundary condition. The scale effect means that the sample dimensions should be 4-6 times larger than the ballast particle (in laboratory tests), to ensure that the results are stable and unaffected [19]. The boundary condition means that when a DEM model represents only a part of the whole system (e.g. halfsleeper track model for the whole ballasted track), the model boundary normally provides different reactions (displacements, forces). For example, when building the half-sleeper track model by DEM, the boundary of the ballast layer is mostly restricted (no displacements) [36]. This will lead to false boundary-ballast reaction, since the boundary imposes larger forces to the contact ballast particles than in reality. When applying the dynamic loads, such boundary condition will result false results due to waves reflection effect.
To solve the issue of the boundary condition, the largescale process simulation test (LPST) model as described in [18] was developed, in contrast to small-scale track model (e.g. ballast box test model) [37]. It has five movable walls at one side to provide consistent pressure stress during the cyclic loading, and in this way the boundary condition is included by moving the lateral walls and providing the lateral deformation.
Therefore, to develop an efficient and accurate method for DEM simulation, this work explores the effects of the rolling resistance linear (RRL) model on the ballast performance of the direct shear test (DST) model and LPST model. The shear strength and settlement of the RRL model is compared with those of the LC model. Specifically, the contact model of the spheres is the RRL, whereas the LC model is used for the non-spherical particles.

Methodology
The DST model and LPST model are developed with the commercial DEM software called Particle Flow Code in 3D (PFC3D). The numerical results derived from these models are compared with those from Ref. [18]. The adoption of two models can fully describe ballast performance such as shear strength, resilience, settlement/permanent deformation. Of these indicators, the shear strength is most widely used and is measured generally by the DSTs [12,28,38]. The settlement/permanent deformation is another key characteristic concerning the performance of ballast assemblies (especially in the field), and is measured by the cutting-edge LPSTs [18]. More importantly, this test model is applicable to the lateral deformation of the ballast assemblies. Figure 1 presents the setup of the DST and the corresponding DEM models. The contact model parameters in Discrete element modelling of railway ballast performance considering particle shape and… 383 the DST model are calibrated using the DST results. Afterwards, we compare the results obtained from two different contact models, i.e. the RRL model and the LC model.

Experimental
In the DSTs, the ballast material is the commonly used aggregate of basalt rock produced in Quarry Pulandian, Dalian, China. The ballast particles have a uniformed shape, sufficient strength, and particle size distribution that follow the British standard [7]. The ballast density is 2530 kg/m 3 . The DST rig consists of three main parts: a steel square box, two hydraulic servo actuators and a computer control system (see Fig. 1a). The steel square consists of an upper steel square box (inner size: 400 mm 9 400 mm 9 200 mm), a lower steel square box (inner size: 400 mm 9 400 mm 9 200 mm) and a steel loading plate (size: 400 mm 9 400 mm 9 20 mm). There is a gap of 10 mm between the upper and lower boxes.
The vertical and lateral hydraulic servo actuators can create the maximum loading of 30 t and 10 t, respectively (Fig. 1a). The vertical actuator can apply the normal force on the steel plate placed on the top of the upper box. This is utilised to provide a constant normal stress in ballast samples. The lateral actuator is used to shift the lower box with a constant speed. The computer control system is utilised to measure vertical and lateral displacements through the linear variable differential transformers (LVDT). It also controls the application of the force or speed of the two hydraulic actuators and records the applied stress.
The ballast particles are placed in the shear box and experience three steps. After placing ballast particles each time, a vibratory compactor is used for compacting the layer. After the third time of compaction, the steel plate (weight of 25.64 kg) is placed on the top of the ballast sample. Then, the direct shear tests are performed at a shearing speed of 2 mm/min under three different normal stresses of 24, 54 and 104 kPa. The final horizontal Fig. 1 Schematic diagram of the applied methodology displacement of the lower DST box is 80 mm (20% shear strain), which is adequate to obtain the peak shear stress.

DST model description
The DST model (Fig. 1b) is utilised to measure the shear strength of the two contact models and four kinds of particle shapes. The porosity of the sample is 0.4, and the particle size distribution (PSD) is based on the above experimental tests. Note that the PSD of all the models remains the same. The model configuration is set as the experimental test configuration (Fig. 1a), including the box size and the applied normal stresses.
The basic contact mode of DEM is a kind of spheresphere contact interactions. Even though in some models the non-spherical particles (clumps) are used, the interaction in the contact areas is still based on the sphere-sphere contact model [39]. However, if the non-spherical particles are present, the number of contact points increases and particle interlocking occurs, finally restraining the particle rotation. On the other hand, if there are simple-shape particles (spheres) with certain rolling friction, it is also possible to result the same effect as the non-spherical particles [35]. Therefore, the rolling friction [39] is used in the DST model.

Contact model and particle shape
In order to determine whether the simple-shape (sphere) particles with the rolling friction can provide the same performance of the model as the complex-shape (clump) particles, two types of contact model are utilised in the model, namely, the LC model and the RRL model. The models with the spheres use the RRL model, while the models with the clumps use the LC model.
The RRL model has one more parameter (rolling friction) than the LC model. In other words, the only difference between the two contact models is the rolling friction. The rolling friction will resist the particle rolling when a force is acting on it. To be more specific, the rolling friction decides the maximum value that equals to the product of the rolling friction with the corresponding normal force. The restriction is defined as rolling stiffness that is assumed as the clockwork spring (Fig. 2), and it increases with the relative rotation.
The four types of the particle shape used in the models are a sphere, a 5-sphere clump, a 12-sphere clump and a 23-sphere clump. Note that one model corresponds to only one type of particle shape. The clump particles are created with the identical template that was obtained by scanning the real ballast particle [40].
In addition, the normal stiffness and shear stiffness (the springs in Fig. 2) are another two parameters in the two contact models that considerably influence the calculation time. Figure 2 describes the LC model.
The calculation time is decided by the timestep calculated based on the two types of stiffness. Specifically, a higher stiffness leads to a smaller timestep, causing more calculation time. The timestep is the smallest time period in simulation, in which the force-displacement law is applied to every updated contact. In other words, a particle moves at a speed in one timestep, and after the time is reached, the forces and displacements are updated. The specific introduction of the timestep can be found in [39]. For this, several values of these two parameters (shear and normal stiffnesses) are selected, and the results are compared for both efficient and accurate simulation.
The properties of the DST model and contact model parameters are listed in Table 1, where four DST models respectively use four types of particle shapes, i.e. the sphere, 5-sphere clump, 12-sphere clump and 23-sphere clump. For the DST model using spheres, the RRL model is utilised, and the particle-particle rolling friction coefficient and the values of the two stiffness (normal and shear) are calibrated. For the DST models using non-spherical particles, the LC model is used and its results are compared with that of the DST model using spheres. Discrete element modelling of railway ballast performance considering particle shape and… 385 The damping applied in the model is local damping (not damping at particle contacts), and the damping value is set according to Ref. [41]. Though there is no consensus in a universal damping value, it has little influence on comparison results. Moreover, it has been proved that different damping values have little influence on the shear strength results when the shear speed is very slow. High Particle-particle friction coefficient 0.5 Particle-particle rolling friction coefficient 0.3 Normal stiffness (N/m) 1 9 10 5 , 2 9 10 5 , 4 9 10 5 and 1 9 10 6 Shear stiffness (N/m) 1 9 10 5 , 2 9 10 5 , 4 9 10 5 and 1 9 10 6 Gravity (m/s 2 ) 9.81  Particle-particle friction coefficient 0. damping value tends to accelerate the formation of equilibrium state.

Model description of large-scale process simulation test
The development of the LPST model refers to the LPST apparatus. As shown in Fig. 3a, the LPST apparatus was designed by Indraratna to develop physical simulation of ''in situ'' railway track. It can contain specimens that are 800 mm long, 600 mm wide, and 600 mm high [42]. Most importantly, one side of the apparatus is made by five movable plates, which can provide consistent principal stresses in the cyclic loading. More explanations on the LPST apparatus can be found in Ref. [3].
The LPST model shown in Fig. 3b includes sleeper, ballast layer and test box. The dimension of the specimen is 800 mm 9 600 mm 9 475 mm, with the ballast thickness (under the sleeper) of 325 mm. The sleeper is constituted by the overlapped spheres (clump), and the ballast particles are simulated with spheres or clumps (same as DST model).
For the model with the spheres, the sample porosity is 0.354, which is larger than the one (0.338) in Ref. [18], whereas the models with three types of clumps have the same porosity (0.338). Even though the porosity is different, the results have shown that their performances can still be the same.
The model properties and parameters are listed in Table 2, including density, friction, stiffness, rolling friction, etc. The movable plates are simulated by walls that keep moving slightly to provide consistent principal stress (10 kPa). The maximum moving speed of the plates is set as 10 mm/s. Four developed LPST models use four different types of particle shapes, i.e. the sphere, 5-sphere clump, 12-sphere clump and 23-sphere clump. For the model with spheres, the RRL model is utilised and the values of   Fig. 4a and b, it can be seen that the peak shear stress increases with the stiffness; however, it increases slightly and stays at around 60-70 kPa after the stiffness is over 4 9 10 5 . Moreover, the peak shear stress is reached with shorter horizontal box displacement, when the contact stiffness is increased from 5 9 10 5 to 1 9 10 7 . Figure 4c and d presents the deformation results at different stiffness values. They illustrate that lower stiffness will cause significant shear contraction, and higher stiffness can lead to deformation results more similar to the experimental test, i.e. 4 9 10 5 , 5 9 10 5 , 1 9 10 6 and 1 9 10 7 . Based on the above results, it can be seen that the shear peak stress increases with the stiffness, but for the model with sphere and low rolling friction (0.1) it does not agree with the experimental shear peak stress. In addition, it is reasonable that the peak shear stress appears when the shear displacement is 30 mm.
Particularly, lower stiffness leads to less computation time. For example, using the spheres with the stiffness at 1 9 10 5 and 1 9 10 7 take the computation time of 433 and 1242 s, respectively. In the same test condition, using the sphere, 5-sphere clump, 12-sphere clump and 23-sphere clump take the computation time of 51, 80, 306 and 400 min, respectively. This means using the spheres is 8 times efficient at most.
In Fig. 4e-h and i-l, the shear stress and deformation results with increasing rolling friction are presented. From the shear stress results, it can be observed that the peak shear stress considerably increase with the rolling friction, and ballast assemblies with higher rolling friction needs larger shear displacement to reach the peak shear stress. Another fact is that with the higher stiffness, the peak shear stress increases at a faster rate than the rolling friction. From the deformation results, it can be seen that the deformation increases with the rolling friction, and high stiffness can cause large deformation change under the increasing rolling friction. Through comparing the experimental results with simulation ones in Fig. 4, we find that the stiffness of 4 9 10 5 can be chosen as the most suitable value for the DST model. Based on the above results, both normal and shear stiffness take the value of 4 9 10 5 . Nevertheless, we design the following simulation conditions to validate this value. The DST simulations under different normal stresses are performed and the results of shear stress and deformation are shown in Fig. 5.  Fig. 6a-c, it can be observed that for the RRL model (i.e. applying rolling friction), the simple sphere and complex shapes (clumps) have similar shear stress results. In addition, it can be observed that the shear stress of the 12-sphere clump is almost the same as that of the 23-sphere clump, but under the normal stresses of 54 and 104, their peak shear stress values are lower than that of the 5-sphere clump. From Fig. 6d-f, it can be observed that the deformation results for using the sphere can better accord with the experimental results than those for using the clump. The 5-sphere clump deformation is higher than the deformation of other two types of clump. The 23-sphere clump has the most realistic shape, but it provides the lowest deformation, which is much lower than the experimental ones.
It is indicated that spheres with rolling friction can replace complex-shaped particles (clumps). Interestingly, it is found that the 23-sphere clump sample provides lower shear stress than the 5-sphere clump sample. This means the interlocks of the 5-sphere clump is stronger than the 23-sphere clump, as some particles link each other to become one big particle. The rolling friction has the same effects as strengthening the contacts and acting as the interlocks.
For further testing the vertical settlement and lateral deformation of the sphere with rolling friction, the LPST model is developed. The simulation results are compared with the clumps and the results from Ref. [18].

Contact force analysis
The contact force analysis is crucial for observing the differences of different particle shapes and contact models at the mesoscopic level. Most of the earlier studies utilised the criterions at the macroscopic level, such as the shear strength in the DST (or triaxial test) [38], and the friction angle in the hopper discharge [37]. They compared the shear stress and strain or the repose angle to present that the parameters in the model can be confirmed. However, different parameters can similarly match the same test results. In other words, large difference in the parameters of the contact model may still lead to similar response. For this reason, the analysis at the mesoscopic level is necessary to perform with respect to the contact force chain, contact force distribution, and coordination number. ). This is because using the spheres reduce the numbers of particles and contacts, and then each contact contributes larger forces. For easy comparison, contact force anisotropy and their distribution are shown in Fig. 8g by the rose diagrams. In Fig. 8, the contact force under the normal stress of 24 kPa with 5 different shearing displacements (0 and 80 mm) is presented, and all the rose diagrams are given in Fig. 15.
In Fig. 8g, the average contact force is calculated from the projected forces. The contact forces are projected to the YZ plane, and the Y-axis directs the shearing direction, as shown in Fig. 1b. The YZ plane is chosen as the shearing direction has the most apparent contact force chain change during shearing.
The average contact force is calculated by averaging the forces within a certain angle range (every five degrees). Specifically, the forces have a direction vector that has an angle to the Y-axis. While 360°are divided every 5°into 72 ranges, the forces with the direction vectors in one range are averaged. The points in every ranges are connected to form one closed curve like the black curve in Fig. 8g. The red curve in Fig. 8g is obtained by smoothing the closed curve for observing the primary orientation more easily. The primary orientation is the purple line drawn by evaluating the direction of the red curve. Specifically, the purple line separates the area into two equal ones.
From this figure, it can be seen that with the increase of the shearing displacement, the primary orientation decreases from around 90°(0 mm) to the lowest (29.7°/ 29.9°/30.6°/36.3°); afterwards, the value slightly increases. Note that the lowest primary orientation with the spheres is approximately the same as that with the clumps except for the 12-sphere clump.
What is more, the contact force with the spheres is 2-2.5 times larger than that with the clumps. However, the average contact forces are approximately the same for the models with the clumps (Fig. 8e, f). For example, in Fig. 8b, the largest average contact force with the sphere is around 100 N, and the smallest is around 40 N. Correspondingly, the models with clumps produce the maximum and minimum values of 40-50 N and 15-25 N, respectively. It can also be seen that the average contact forces increase with the shearing displacement.
It is significant to find that the average contact force of the spheres is 2-2.5 times larger than that of the clumps. This is because the contacts of the spheres are approximately half of the clumps, and thus every contact bear more shearing stress. Alternatively, every contact of the spheres are strengthened. The contact number of each particle can be presented by the coordination number, which will be discussed in the following section. Fig. 9b. The coordination number is the average number of active contacts for each particle. The coordination number is calculated by the use of the particles that lie at the shearing zone within the four measurement spheres, as shown in Fig. 9a. As the shearing zone is the most important position to produce the shearing stress, the particles at the shearing zone have the most obvious movements. The radius of the measurement spheres is 0. From Fig. 9, it can be seen that the coordination number increases as the particle shape is more complex (from spheres to clumps). Moreover, with the increase of the  Discrete element modelling of railway ballast performance considering particle shape and… 395 Fig. 10 Particle rotation of using the spheres or clumps normal stress (24,54, and 104 kPa), the coordination number also increases, as the assemblies are more compacted. The 23-sphere clumps can produce approximately twice coordination number than the spheres. The coordination number results demonstrate that the contacts of the spheres are less than those of the clumps.

Particle rotation
In order to confirm the effects of the rolling resistance on the particle rotation, the sphere model with the RRL model is compared with the clump model with the LC model, as shown in Fig. 10.  Discrete element modelling of railway ballast performance considering particle shape and… 397 circle colour helps to distinguish the rotation magnitude (Fig. 10a). The figure shows part of the results that are obtained under the normal stress 104 kPa, and all the particle rotation results are given in Fig. 16 (''Appendix''). The particle rotation is calculated by Eq. (1) [39]. In the equation, the Euler angles are utilised to calculate the particle rotation (i.e. /, h and w), which present the precession rotation, nutation rotation and intrinsic rotation, respectively.
From the figures, it can be observed that the particle rotation of the sphere model is almost the same as that of the clump model. Specifically, the largest rotation (over 180°) appears at the similar positions, which are the left side of the upper shear box and the right side of the lower shear box, and both of them are near the shearing interface. In addition, most of the large circles (green, purple, and red) appear along the diagonal line of the shear box and the line is approximately perpendicular to the contact chain direction.

Large-scale process simulation test model
3.2.1 Stiffness and particle shape Figure 11 presents the applied stress vs vertical displacement with four kinds of particles and different normal and shear stiffnesses. From the figure, it can be observed that the elastic deformation and plastic deformation reduce as the stiffness increase. In addition, the results at the stiffness of 4 9 10 5 cannot accord with the results in Ref [18], where the elastic deformation and plastic deformation are within 0.5 mm (Fig. 11a). After the comparison, it is found that the results that correspond to the spheres or 5-sphere clumps with the stiffness of 1 9 10 7 or 1 9 10 8 N/m (Fig. 11c, e) can approximately accord with those in [18]. This proves that using one set of contact model parameters may not be fit to all the tests, despite testing on the same material. Even though the DST is a well-known method for confirming the   Additionally, it also demonstrates that the sphere model with rolling friction can have the same or even better performance compared with the clump model. Particularly, the hypothesis of the interlocks (Sect. 3.1.2) can be proved again by comparing the results in Fig. 11c, e, g and i. Specifically, the spheres with rolling friction has less vertical deformation than the clumps. Moreover, the 5-sphere clump has less deformation than the 12-sphere clump and 23-sphere clump.
The stiffness cannot be 1 9 10 7 or 1 9 10 8 due to two facts: (1) the initial stage (before applying loadings) of ballast particle in the numerical simulation plays an important role on the first a few cycles. This can be illustrated from Fig. 11, showing that the first cycle has the largest deformation than the other following cycles. In addition, there is large deformation in the first 5 cycles; afterwards, the deformation becomes small and stable, which means the contacts between the sleeper and the ballast particles become more and the ballast particles near the sleeper are rapidly compacted. (2) The elastic deformation and plastic deformation values have a large range due to the discrete nature of railway ballast [43,44]. In response to this, the lateral displacement and stress results of five movable plates are presented and compared with the results in [18].

Lateral displacement and stress
The lateral displacement results of the five movable plates are shown in Fig. 12, and their stress results are shown in Fig. 13. From Fig. 12a-d, it can be seen that the lateral displacements reduce with the stiffness increase. According to Fig. 12q, after 100 cycles, the lateral displacements of all five movable plates are within 5 mm. To match the test lateral displacements, the stiffness values of 4 9 10 5 and 1 9 10 6 are not suitable for the LPST model.
Particularly, the stress from the five movable plates become a constant (10 kPa) for the stiffness of 1 9 10 8 , as shown in Fig. 13d, h, l, and p. This demonstrates that the sphere with rolling friction can provide adequately reliable results. To be more specific, in Fig. 13d, the stresses become stable from 32 kPa to 10 kPa after one cycle, while the other particle shapes stabilize from much higher values, 75 kPa (5-sphere clump), 105 kPa (12-sphere clump) and 185 kPa (23-sphere clump). In addition, the 12-sphere clump and 23-sphere clump need more cycles to become stable, 11 and 9 cycles, respectively.

Conclusions
To increase the DEM simulation efficiency, the DEM models for the DST and the LPST are developed and applied to analyse the ballast performance in terms of shear strength and deformation. The efficiency of different contact model types and particle shapes are studied. The numerical results are compared with the experimental results and results from the literature. From the results and discussion, the following conclusions can be summarised: Nevertheless, these models have quite different particle scale performance; e.g. there are still large discrepancies among the particle performances (movements) in these models.
The LPST model was tested in 15 loading cycles, and it is necessary to observe the long-term deformation performance after thousands of cycles. Additionally, the spheres with the rolling friction can lead to the same results as those from the tests; however, the detailed reasons and mesoscopic mechanics at the particle contacts can be analysed deeper in a way of the DEM simulations. It should be emphasised that the contact model parameters need further investigations to settle the most suitable ones. Finally, the particle degradation has not been considered in this work, and further studies will be performed in this respect.
Universities (111 Project, Grant No. B16041). We want to thank the support during my work at the International Joint Laboratory on Railway Engineering System Dynamics in Southwest Jiaotong University. We would like to thank Xu Zhang from Guangdong University of Technology, China and Shunying Ji from Dalian University of Technology, China, for their contribution to this work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.

Fig. 14 continued
Discrete element modelling of railway ballast performance considering particle shape and… 403  Discrete element modelling of railway ballast performance considering particle shape and… 405