Discrete Element Modelling of Fractal Behavior of Particle Size Distribution and Breakage of Ballast under Monotonic Loading

: Breakage of particles has a great inﬂuence on the particle size distribution (PSD) and the associated mechanical behavior of ballast under train loads. A discrete element method (DEM) simulation of triaxial testing under monotonic loading was carried out using FRM (fragment replacement method) breakable particles as ballast and a ﬂexible shell model as membrane. The coupled model was validated by comparing the load-deformation responses with those measured in previous experiments and was then used to analyze the contact orientations and the distribution of particle breakage from a micromechanical perspective. The simulation results show that higher conﬁning pressure and larger axial strain may increase the grain breakage ( B g ) and the fractal dimension ( D ) of ballast. It was observed that most breakage was ﬁrst-generation breakage, and that the proportions of the second- to ﬁfth-generation breakage decreased successively. Moreover, as the axial strain or conﬁning pressure increased, the percentage of small particle fragments increased in correspondence with the PSD curves that remained concave upwards, as the fractal dimension D of PSD increased. In addition, the evolution of D exhibited a linear correlation with grain breakage B g . Contrarily, a quadratic curve relation between D and volumetric strain was exhibited under different axial strain stages. Therefore, D has the potential to be a key indicator to evaluate the degree of ballast crushing and PSD degradation, which may contribute to better decision making concerning track bed maintenance.


Introduction
Ballast is a granular material used for construction of railway beds. It can be characterized as having an extremely irregular particle shape with high porosity and strength. Ballast helps to transmit and distribute axle loads from the sleepers to the sub-ballast layer [1]. The dynamic stress of a passing train degrades and fouls the ballast layer which usually experiences continuous breakage, directly resulting in track settlement, densification, and the need for frequent maintenance. Ballast breakage may lead to a sharp increase in the number of fragments, weakening interlocking and impacting mechanical performances [2]. In practice, ballast needs to be maintained and replaced when its quality deteriorates to a certain extent.
Recently, fractal theory has been widely used in describing and studying the particle size distribution (PSD) of granular material [3]. Compared with traditional methods based on Euclidean geometry, fractal theory describes PSD well because it can characterize the local shape, size, and structure of particles, and any similarities among them [4]. It is beneficial to analyze the pore structure of materials by combining more than one method, considering such factors as the relationship between fractal characteristics, or the physical or chemical attributes of materials [5]. Early studies revealed that fractal theory could be used to study the fractal characteristics of soil [4][5][6], sand [7], rockfill [8][9][10] and other material [11,12]. Compared with natural soil or quartz sand of smaller size and wider grading distribution, layers of ballast are usually composed of medium-to-coarse gravelsized particles (10-60 mm) whose uniformity coefficients vary from 1.5 to 3.0 according to the ballast specifications [13] which can be regarded as uniformly sized granular material. As more and more small fragments were produced by ballast breakage under train loads, the PSD curve of ballast also changed, with a fractal dimension which was directly proportional to the fragment content. Therefore, the fractal theory could be used to study PSD evolution and ballast breakage. However, studies on the effect of the fractal distribution of particle size on the characteristics of ballast are still not enough in themselves.
The discrete element method (DEM) has been extensively used in railway engineering and also pavement engineering [14][15][16]. In order to analyze the evolution of fractal PSD using the fragment replacement method (FRM), the particle is represented by a group of smaller ones without bonding as implemented in [16][17][18][19][20][21]. McDowell et al. [16] adopted octahedral shear stress as the fracture criterion and replaced broken spheres with smaller ones with overlap in their successful investigation of the evolution of fractal PSD for silica sand. Ciantia et al. [17] used a multigenerational DEM approach to simulate soil breakage in a one-dimensional compression test, and obtained a limit fractal dimension of 2.47. Zhou et al. [18] found that both particle size and particle contact force were distributed in a fractal manner by one-dimensional compression test simulation, which is consistent with the findings of Vallejo and Lobo-Guerrero [19]. De Bono et al. [20] conducted one-dimensional normal compression test simulation with different fracture criteria, and pointed out that both octahedral shear stress and the largest normal contact stress can lead to the correct normal compression lines and fractal PSDs. Li et al. [21] quantified the evolution of microstructures of carbonate sands under mechanical loading and found that the fractal dimension of crack networks increased with external loading due to crack branching via cleavage.
Although the FRMs used in the above studies can accurately simulate particle breakage based on DEM, most of the parent breakable particles involved were spherical and so could hardly simulate the interlocking between particles [16][17][18][19][20][21]. However, particle shape plays a key role in the fractal behavior of ballast with irregular shapes [22,23]. Therefore, in order to reveal the fractal characteristics of ballast, a DEM-FDM (finite differential method) coupled simulation of triaxial testing under monotonic loading using an irregular-shaped breakable ballast model is carried out in this study. The effects of axial strain and confining pressure on ballast breakage and the fractal behavior of ballast PSD are investigated. The relationships between fractal dimensions, grain breakage and the stress states of ballast are then established.

Ballast Breakage Model Using FRM
Firstly, the DEM study using FRM to simulate particle breakage should include a reasonable fracture criterion and a fragment replacement mode. They should be defined to check (1) whether a particle should break and (2) how any individual particle does break.
It is widely accepted that the failure of spherical particles under compression is in fact a tensile failure. Jaeger [24] suggested that the 'tensile strength' of rock grains can be measured by diametral compression between flat platens. For a grain of diameter d under a diametral force F, a characteristic tensile stress σ induced within it can be defined as For railway ballast with typical sizes ranging from 10 to 60 mm, Lim et al. [25] conducted indirect diametral compression and measured the tensile strengths of ballast between flat platens pointing out that Weibull strength σ 0 is a function of the particle size d. This is usually expressed in the form where σ 0 is the value of σ such that 37% of the blocks survive the test intact; while b is a constant and a function of the statistical distribution of flaws in the material, describing the size-hardening law. Typical values of b range from −0.69 to −0.40 for different types of ballasts. Following the single-particle crushing test results [25], the characteristic tensile strength σ 0 of ballast size d 0 (19.7 mm) is 14.7 MPa, and the value of b is −0.56, which is used to attribute random strengths to new fragments and other initial particles.
In order to determine the correct behavior of ballast under normal compression, De Bono and McDowell [20] investigated various possibilities for the fracture criterion, including octahedral shear stress, major principal stress within a particle, mean stress, and stress calculated from the maximum contact force on a particle. The criterion based on the maximum contact force (F max ) was shown to better match the experiments and simulations, such that According to Raymond et al. [1], the process of ballast degradation due to loading can occur in three ways: by the grinding of small-scale asperities, by breakage of angular projections, and by the splitting of particles into approximately equal parts. In addition, an irregular clump produces more contact compared to a spherical particle, leading to a lower maximum contact force. Hence, these particles will be replaced by new fragments when the F max stress exceeds than half of the characteristic strength.
From the perspective of computational efficiency, in this study, each particle was allowed to split into several smaller particles, maintaining conservation of mass when the value of σ Fm was greater than or equal to its Weibull strength. Previous studies have used spheres or disks to simulate particle breakage within the bonded parent sphere, where new sphere fragments to be contained overlap sufficiently [16,18,20] or tangentially to each other [17,19]. As shown in Figure 1a, new sphere fragments overlap sufficiently to be contained within the bounding parent sphere, with the axes joining their centers aligned along the direction of the minor principal stress. This process produces undesirable local pressure spikes during breakage and continues until internal forces are completely released and the system reaches equilibrium [16].

Simulations of Monotonic Triaxial Test
Based on the experimental monotonic triaxial tests conducted by Indraratna et al. [26], DEM-FDM coupled simulations of monotonic triaxial tests were carried out under different confining pressures ranging from 60 kPa to 240 kPa. Figure 2 shows the DEM-FDM model of a triaxial test with dimensions of 300 mm diameter and 600 mm height. It  It is therefore acknowledged that it is impossible to simulate a perfectly realistic ballast fracture mechanism using disks or spheres [16]; therefore, eight-pebble clumps of irregular shape were adopted for this study. The fragment replacement mode is depicted in Figure 1b, which shows that irregular particles were replaced by smaller fragments within the bounding parent particles without overlaps. During the first-generation breakage, the eight-pebble clump might split into a four-pebble clump and four single spheres. All the fragments could continue to break and generate smaller particles. Except for single spheres, most of the particles could break and generate fragments without overlap, which means internal forces could be releases and the system could reach equilibrium immediately. Hence this fragment replacement mode can greatly improve calculation efficiency, and is therefore used in this study.

Simulations of Monotonic Triaxial Test
Based on the experimental monotonic triaxial tests conducted by Indraratna et al. [26], DEM-FDM coupled simulations of monotonic triaxial tests were carried out under different confining pressures ranging from 60 kPa to 240 kPa. Figure 2 shows the DEM-FDM model of a triaxial test with dimensions of 300 mm diameter and 600 mm height. It consists of 971 eight-pebble clumps as breakable ballast using FRM in PFC3D (particle flow code in three dimensions) with a flexible boundary as membrane using shell model in FLAC3D (fast Lagrangian analysis of continua in three dimensions). The top and bottom plates using disk walls have the same dimensions as the triaxial test, to maintain a consistent area of contact with the ballast sample. The eight-pebble clump has an irregular shape with an aspect ratio of 1.83, as expressed by the ratio of the particle length to the particle width. In accordance with the PSD of ballast in triaxial tests, the PSD used in these simulations is shown in Table 1.   During sample generation, the particles generated at random orientation by th Distribute command. Firstly, the sample was divided into five equal layers, each was generated in turn and allowed to fall under gravity, then compacted as perform the lab. Once the sample attained equilibrium, a flexible membrane model was gene to replace the cylinder wall. The membrane model was linked with the top and bo walls by each node. The shell model as a flexible membrane consisted of 1440 facet which were assumed to be triangles of uniform thickness lying between three n points. Each shell structural element was defined by its geometric and material prope The parameters of FRM ballast particles and shell elements are listed in Table 2.  During sample generation, the particles generated at random orientation by the Ball Distribute command. Firstly, the sample was divided into five equal layers, each layer was generated in turn and allowed to fall under gravity, then compacted as performed in the lab. Once the sample attained equilibrium, a flexible membrane model was generated to replace the cylinder wall. The membrane model was linked with the top and bottom walls by each node. The shell model as a flexible membrane consisted of 1440 facet units which were assumed to be triangles of uniform thickness lying between three nodal points. Each shell structural element was defined by its geometric and material properties. The parameters of FRM ballast particles and shell elements are listed in Table 2. Table 2. Simulation parameters in PFC3D and FLAC3D.

Input Parameters Value
Normal and shear stiffness of particles: N/m 2 × 10 6 Friction coefficient of clumps 0. In order to reduce the influence of the initial arrangement of particles on the simulation results, the same sample was used for different loading conditions. Once the assembly generation procedure was completed, artificial isotropic compression was applied to the plates and membrane to achieve the required stress state. The assigned uniform pressure, ranging from 60 to 240 kPa, was applied to all shell elements of the membrane model, so that the required confining pressure was obtained. Monotonic loading was applied to both the top and bottom disk walls with a constant speed of 0.02 m/s, and the monotonic load test continued until the axial strain reached 20%. During the loading process, data relating to porosity, the stress-strain curve, volumetric strain and particle breakage were all recorded. Figure 3 shows the simulations and experimental results, in terms of deviator stress and of volumetric strain against axial strain under different confining pressures. Comparing the experimental results with the simulation results, a trend of stress-strain relationship is seen to be well simulated. With the development of the axial strain, the deviator stress gradually increases and attenuates after reaching a relatively stable value, showing the strain softening behavior and the maximum deviator stress increases with the increase in confining pressure. However, the shear strength of the simulated assembly is overestimated in comparison with the experimental one. This could be caused by the irregular shape of particles with an aspect ratio of 1.83, which were recognized as flaky particles. Similarly, after having conducted a set of triaxial ballast tests to investigate the effects of ballast shape on ballast performance, Roner [27] found that randomly placed flaky material exhibited higher deviator stress and higher angles of internal friction than non-flaky material at the same void ratio. This offers a possible explanation as to why the deviator stress was higher than expected. In addition, an orientation parallel to the failure plane will cause a substantial strength reduction when a significant proportion of the particles are flaky.

Stress-Strain Response
ballast shape on ballast performance, Roner [27] found that randomly placed flaky material exhibited higher deviator stress and higher angles of internal friction than non-flaky material at the same void ratio. This offers a possible explanation as to why the deviator stress was higher than expected. In addition, an orientation parallel to the failure plane will cause a substantial strength reduction when a significant proportion of the particles are flaky.
(a) (b)   However, the purpose of this study is not to resolve the particle shape problem, but rather to adopt the best approach to investigate the fragment replacement mode and particle crushing under shear. Therefore, this study applies simplified and irregular shaped particles, and focus on the fractal fragmentation and gradation evolution of particles. Figure 3b shows that the sample exhibits shear shrinkage in the initial stage, but with the increase in axial deformation, an obvious shear expansion is observed, especially under low confining pressure. As the confining pressure decreases from 240 kPa to 60 kPa, the dilatancy effect becomes more obvious. Volumetric strain agrees reasonably well with the experimental results under a high confining pressure of 240 kPa where breakage is more pronounced. This is probably attributable to the generation of smaller size fragments from the split particles. Using more and smaller fragments to form the origin clumps may provide more realistic volumetric strains at lower confining pressures. Therefore, further investigation of the effects of fragment replacement mode is needed. Figure 4 shows the effect of different confining pressures on the peak friction angle ϕ p of the ballast sample, and both the simulation and experimental results show that the peak friction angle decreases with the increase in confining pressure. Moreover, the simulation using uniform flat shape particles shows a higher angle of internal friction than the experimental results. This is due to the fact that the FRM particle only undergoes particle splitting at higher stress thresholds and is unable to undergo corner breakage and surface abrasion to release internal forces at the initial strain stage. Figure 4 shows the effect of different confining pressures on the peak friction angle φp of the ballast sample, and both the simulation and experimental results show that the peak friction angle decreases with the increase in confining pressure. Moreover, the simulation using uniform flat shape particles shows a higher angle of internal friction than the experimental results. This is due to the fact that the FRM particle only undergoes particle splitting at higher stress thresholds and is unable to undergo corner breakage and surface abrasion to release internal forces at the initial strain stage. It can be seen from Figure 5a that the deformation of the sample is allowed to be uneven, and the expansion of the middle part is the most obvious after using the flexible membrane model. This verifies that the flexible membrane model is better than the servocontrolled cylinder walls for simulating the triaxial test. In order to study the development of porosity during loading, three measurement spheres with radii of 85, 120, and 85 mm were arranged in the sample, as shown in Figure 5b. Figure 5b shows the evolution of porosity with axial strain in the upper, middle and lower regions under 240 kPa confining pressure. As the axial strain increases to 12%, the porosity of each part exhibits the characteristics of first decreasing and then increasing. Subsequently, the porosities in the upper and lower parts are stable, but the porosity in the middle region continues to increase. It can be seen from Figure 5a that the deformation of the sample is allowed to be uneven, and the expansion of the middle part is the most obvious after using the flexible membrane model. This verifies that the flexible membrane model is better than the servocontrolled cylinder walls for simulating the triaxial test. In order to study the development of porosity during loading, three measurement spheres with radii of 85, 120, and 85 mm were arranged in the sample, as shown in Figure 5b. Figure 5b shows the evolution of porosity with axial strain in the upper, middle and lower regions under 240 kPa confining pressure. As the axial strain increases to 12%, the porosity of each part exhibits the characteristics of first decreasing and then increasing. Subsequently, the porosities in the upper and lower parts are stable, but the porosity in the middle region continues to increase. The evolution process of the particle displacement vector of the sample under the confining pressure of 240 kPa is shown in Figure 6. During the loading process, the particles within the top and bottom half move toward the horizontal plane in the middle. When axial strain reaches 4%, due to mutual resistance, the intermediate particles start to move outward horizontally to produce a shear sliding surface, and the angle between the sliding surface and the horizontal direction gradually increases. The evolution process of the particle displacement vector of the sample under the confining pressure of 240 kPa is shown in Figure 6. During the loading process, the particles within the top and bottom half move toward the horizontal plane in the middle. When axial strain reaches 4%, due to mutual resistance, the intermediate particles start to move outward horizontally to produce a shear sliding surface, and the angle between the sliding surface and the horizontal direction gradually increases.
Under a confining pressure of 240 kPa, the polar coordinate distribution of contact numbers and average contact forces within the sample is plotted as shown in Figure 7. With the increase in axial strain, the average contact force in the sample first increases and then decreases between 12% and 20% axial strain, which is similar to the substantial strength reduction in the deviator stress after peaking shown in Figure 3b  The evolution process of the particle displacement vector of the sample under the confining pressure of 240 kPa is shown in Figure 6. During the loading process, the particles within the top and bottom half move toward the horizontal plane in the middle. When axial strain reaches 4%, due to mutual resistance, the intermediate particles start to move outward horizontally to produce a shear sliding surface, and the angle between the sliding surface and the horizontal direction gradually increases. Under a confining pressure of 240 kPa, the polar coordinate distribution of contact numbers and average contact forces within the sample is plotted as shown in Figure 7. With the increase in axial strain, the average contact force in the sample first increases and then decreases between 12% and 20% axial strain, which is similar to the substantial strength reduction in the deviator stress after peaking shown in Figure 3b. Most of the internal contacts of the sample are distributed in the range of ± 30° from the vertical plane in the initial stage; the contact number in the vertical direction then decreases with the loading, while the contact number in the horizontal direction increases. At the end of loading, the average contact force and contact number in the vertical direction are almost 3.36 and 2.49 times larger than those in the horizontal direction.

Particle Breakage
As shown in Figure 8, the number of fragments could be ignored at the initial compaction stage (4% axial strain) before increasing linearly above 8% axial strain. This indicates that the sample becomes denser in the compaction stage due to particle rearrangement and a reduction in the number of pores. When the porosity reaches the minimum value, the particle interlocking phenomenon is further exacerbated and a large number of particles are broken off by shearing. In addition, more particle breakage occurs under higher confining pressure.

Particle Breakage
As shown in Figure 8, the number of fragments could be ignored at the initial compaction stage (4% axial strain) before increasing linearly above 8% axial strain. This indicates that the sample becomes denser in the compaction stage due to particle rearrangement and a reduction in the number of pores. When the porosity reaches the minimum value, the particle interlocking phenomenon is further exacerbated and a large number of particles are broken off by shearing. In addition, more particle breakage occurs under higher confining pressure.
In accordance with experiment [26], the differences in the percentage by weight of each particle size fraction before and after simulation (∆W k ) are plotted against the aperture of the lower sieve corresponding to that fraction, as shown in Figure 9. The equivalent diameters and weights of each particle size are calculated by tracing all particles and extracting the volumes in simulation. For aperture sizes between 18.9 mm and 29 mm, the percentage change in particle size decreased and formed a main trough. Meanwhile, For aperture sizes between 18.9 mm and 53 mm, the percentage change of particle size increased and formed a main peak. Figure 8, the number of fragments could be ignored at the initial compaction stage (4% axial strain) before increasing linearly above 8% axial strain. This indicates that the sample becomes denser in the compaction stage due to particle rearrangement and a reduction in the number of pores. When the porosity reaches the minimum value, the particle interlocking phenomenon is further exacerbated and a large number of particles are broken off by shearing. In addition, more particle breakage occurs under higher confining pressure. In accordance with experiment [26], the differences in the percentage by weight of each particle size fraction before and after simulation (∆Wk) are plotted against the aperture of the lower sieve corresponding to that fraction, as shown in Figure 9. The equivalent diameters and weights of each particle size are calculated by tracing all particles and extracting the volumes in simulation. For aperture sizes between 18.9 mm and 29 mm, the percentage change in particle size decreased and formed a main trough. Meanwhile, For aperture sizes between 18.9 mm and 53 mm, the percentage change of particle size increased and formed a main peak. A measure of grain breakage (Bg) equal to the sum of positive values of ∆Wk is now plotted and expressed as a percentage, as introduced by Marsal [28] for granular materials. The Bg under 240 kPa confining pressure increases with increasing axial strain. In particular, the largest increment of Bg can be found the stage from 8% to 12% axial strain, which also corresponds to the peak in the growth rate of crushing, as shown in Figure 8. Moreover, with a large number of particles broken in this phase, the peak deviator stress is gradually deceased when the axial strain reaches 8% as shown in Figure 3a. Figure 10b shows levels of grain breakage Bg after loading under different confining pressures. It can be observed that more particles are broken off with increasing confining pressure in both the experimental test and the simulations, which proves the acceptability of the fracture criterion and fragment replacement mode used in the simulation study. It should be noted that the breakage in simulation is underestimated compared with the experimental test. This could be explained by the fact that the simple splitting model does not take particle abrasion and small corner breakage into consideration. A measure of grain breakage (B g ) equal to the sum of positive values of ∆W k is now plotted and expressed as a percentage, as introduced by Marsal [28] for granular materials. The B g under 240 kPa confining pressure increases with increasing axial strain. In particular, the largest increment of B g can be found the stage from 8% to 12% axial strain, which also corresponds to the peak in the growth rate of crushing, as shown in Figure 8. Moreover, with a large number of particles broken in this phase, the peak deviator stress is gradually deceased when the axial strain reaches 8% as shown in Figure 3a. Figure 10b shows levels of grain breakage B g after loading under different confining pressures. It can be observed that more particles are broken off with increasing confining pressure in both the experimental test and the simulations, which proves the acceptability of the fracture criterion and fragment replacement mode used in the simulation study. It should be noted that the breakage in simulation is underestimated compared with the experimental test. This could be explained by the fact that the simple splitting model does not take particle abrasion and small corner breakage into consideration.

As shown in
For increasing levels of axial strain, the quantities and positions of fragments under 240 kPa confining pressure are plotted in Figure 11. The generations of breakage are defined from first generation to fifth generation as shown in Figure 1. At the initial stage of 8% axial strain, almost all of the fragments belong to the first generation, and the lateral bulging of the sample is not obvious. When ε a = 12%, the fragments coming from the second and third generations begin to generate at locations near the top and bottom plates. As ε a increases to 16%, and to 20%, new fragments coming from the fourth and fifth generations, respectively, start to appear in the sample. Apart from the area near the top and bottom plates, most of these fragments (second-third generation) are distributed along the shear surface as shown in Figure 11e, which is consistent with the distribution of particle displacement in Figure 6e. In summary, 160 of 971 eight-ball clumps were split and a total number of 1341 fragments were generated. Among these, 1237 fragments were of the first generation, accounting for 92.2% of all breakage and located all over the sample, while 104 fragments of the second to fourth generation accounted for 7.8% of all breakage and are more concentrated in the upper half of the sample. Therefore, the first-generation breakage is dominant, and the proportions of breakage decrease from the second to the fifth generations, successively A measure of grain breakage (Bg) equal to the sum of positive values of ∆Wk is now plotted and expressed as a percentage, as introduced by Marsal [28] for granular materials. The Bg under 240 kPa confining pressure increases with increasing axial strain. In particular, the largest increment of Bg can be found the stage from 8% to 12% axial strain, which also corresponds to the peak in the growth rate of crushing, as shown in Figure 8. Moreover, with a large number of particles broken in this phase, the peak deviator stress is gradually deceased when the axial strain reaches 8% as shown in Figure 3a. Figure 10b shows levels of grain breakage Bg after loading under different confining pressures. It can be observed that more particles are broken off with increasing confining pressure in both the experimental test and the simulations, which proves the acceptability of the fracture criterion and fragment replacement mode used in the simulation study. It should be noted that the breakage in simulation is underestimated compared with the experimental test. This could be explained by the fact that the simple splitting model does not take particle abrasion and small corner breakage into consideration. For increasing levels of axial strain, the quantities and positions of fragments under 240 kPa confining pressure are plotted in Figure 11. The generations of breakage are defined from first generation to fifth generation as shown in Figure 1. At the initial stage of 8% axial strain, almost all of the fragments belong to the first generation, and the lateral bulging of the sample is not obvious. When εa = 12%, the fragments coming from the second and third generations begin to generate at locations near the top and bottom plates. As εa increases to 16%, and to 20%, new fragments coming from the fourth and fifth generations, respectively, start to appear in the sample. Apart from the area near the top and bottom plates, most of these fragments (second-third generation) are distributed along the shear surface as shown in Figure 11e, which is consistent with the distribution of particle displacement in Figure 6e. In summary, 160 of 971 eight-ball clumps were split and a total number of 1341 fragments were generated. Among these, 1237 fragments were of the first generation, accounting for 92.2% of all breakage and located all over the sample, while 104 fragments of the second to fourth generation accounted for 7.8% of all breakage and are more concentrated in the upper half of the sample. Therefore, the first-generation breakage is dominant, and the proportions of breakage decrease from the second to the fifth generations, successively

Gradation Evolution and Fractal Distribution
Turcotte et al. [3] demonstrated that any initial distribution of soil particles would tend towards a self-similar fractal distribution with increasing particle fragmentation. Railway track bed typically consists of ballast of size 10-60 mm in different ballast grading specifications, generally larger than soil or calcareous sand. Considering the ballast breakage results in the previous section, the fractal evolution of the PSD against the axial strain and confining pressure were also studied. Figure 12a shows the evolution of PSD under monotonic loading and a confining pressure of 240 kPa. It can be seen that, as the axial strain increases, the percentage of small particle fragments ranging from 10 mm to 35 mm increases, with the grading curves remaining concave upwards. This can be attributed to the supplement of aggregates induced by particle breakage from the upper size ranges. Figure 12b shows the comparison of PSDs under different confining pressures, ranging from 60 kPa to 240 kPa. It can be seen that all the sample PSDs after test are shifted left where smaller aggregates were generated, when compared with the initial grading. With increases in confining pressure, the sample grading shifts more towards the left, especially under the highest confining pressure of 240 kPa.

Gradation Evolution and Fractal Distribution
Turcotte et al. [3] demonstrated that any initial distribution of soil particles would tend towards a self-similar fractal distribution with increasing particle fragmentation. Railway track bed typically consists of ballast of size 10-60 mm in different ballast grading specifications, generally larger than soil or calcareous sand. Considering the ballast breakage results in the previous section, the fractal evolution of the PSD against the axial strain and confining pressure were also studied. Figure 12a shows the evolution of PSD under monotonic loading and a confining pressure of 240 kPa. It can be seen that, as the axial strain increases, the percentage of small particle fragments ranging from 10 mm to 35 mm increases, with the grading curves remaining concave upwards. This can be attributed to the supplement of aggregates induced by particle breakage from the upper size ranges. Figure 12b shows the comparison of PSDs under different confining pressures, ranging from 60 kPa to 240 kPa. It can be seen that all the sample PSDs after test are shifted left where smaller aggregates were generated, when compared with the initial grading. With increases in confining pressure, the sample grading shifts more towards the left, especially under the highest confining pressure of 240 kPa. According to the study by Tyler et al. [11], the fractal grading can be described as where Fu(d) is the percentage finer, dmax is the maximum diameter for the given grain assembly, MT is the total mass of particles, δ is the particle size, and M(δ < di) is the cumulative mass of the particles whose grain size is smaller than di. The fractal dimension D can be determined by: The PSD curve presents the correlations between the percentages passing and normalized particle diameters which gradually close to a linear line with the log-log scales. Figure 13a shows that the fitting lines shift upward, and the grading curve gradually becomes linear as the axial strain increases. Meanwhile, the linear portion of the curve from which the fractal dimension can be obtained increases in length, suggesting a more reliable value. Moreover, the slope becomes steeper (i.e., the fractal dimension D increases) with increasing confining pressure, with the slope appearing to become constant as shown in Figure 13b. The size ratio of the largest and smallest particle is much smaller than that of soil or quartz sand, with values of Cu ranging from 1.5 to 3.0, which can be regarded as granular material of uniform size. This is quite different from granular materials in natural state because the grading ballast used in railway is qualified by the ballast standard, leading to a much smaller initial value of D (0.2-0.3). Predictably, a continued compression will cause a further crushing and increasing of fine particles, leading to a further fractal.
With the faster and heavier of the passing train, the ballast layer are subjetced to settle down and crushing under the long-term cyclic loading, and the fractal dimension gradually increases. Ballast will need to be maintained or replaced once the fractal dimension reaches a constant value. Further experimental research is required to predict the critical value of the fractal dimension which indicates when ballast needs to be maintained. Figure 14 shows the evolution of fractal dimension against grain breakage under different confining pressures. It shows that the fractal dimension of the sample increases with rising confining pressure and also axial strain. The relationship between D and Bg is almost linear. Therefore, the fractal dimension can be well used to evaluate the degree of crushing and deterioration of ballast. According to the study by Tyler et al. [11], the fractal grading can be described as where F u (d) is the percentage finer, d max is the maximum diameter for the given grain assembly, M T is the total mass of particles, δ is the particle size, and M(δ < d i ) is the cumulative mass of the particles whose grain size is smaller than d i . The fractal dimension D can be determined by: The PSD curve presents the correlations between the percentages passing and normalized particle diameters which gradually close to a linear line with the log-log scales. Figure 13a shows that the fitting lines shift upward, and the grading curve gradually becomes linear as the axial strain increases. Meanwhile, the linear portion of the curve from which the fractal dimension can be obtained increases in length, suggesting a more reliable value. Moreover, the slope becomes steeper (i.e., the fractal dimension D increases) with increasing confining pressure, with the slope appearing to become constant as shown in Figure 13b. The size ratio of the largest and smallest particle is much smaller than that of soil or quartz sand, with values of C u ranging from 1.5 to 3.0, which can be regarded as granular material of uniform size. This is quite different from granular materials in natural state because the grading ballast used in railway is qualified by the ballast standard, leading to a much smaller initial value of D (0.2-0.3). Predictably, a continued compression will cause a further crushing and increasing of fine particles, leading to a further fractal.
With the faster and heavier of the passing train, the ballast layer are subjetced to settle down and crushing under the long-term cyclic loading, and the fractal dimension gradually increases. Ballast will need to be maintained or replaced once the fractal dimension reaches a constant value. Further experimental research is required to predict the critical value of the fractal dimension which indicates when ballast needs to be maintained.     Figure 14 shows the evolution of fractal dimension against grain breakage under different confining pressures. It shows that the fractal dimension of the sample increases with rising confining pressure and also axial strain. The relationship between D and B g is almost linear. Therefore, the fractal dimension can be well used to evaluate the degree of crushing and deterioration of ballast.
The relationship between D and the volumetric strain ε v for these four samples under different confining pressures can be plotted in Figure 15. At the initial stage (4% axial strain), the samples are compressed with an almost constant value of D; a linear relationship between volumetric strain and D can then be observed above 8% axial strain under different confining pressures, which corresponds with the maximum sample dilation angle as shown in Figure 3b. A quadratic curve relation between the volumetric strain and fractal dimension can be obtained for different levels of axial strain, as shown in Figure 15   The relationship between D and the volumetric strain εv for these four samples under different confining pressures can be plotted in Figure 15. At the initial stage (4% axial strain), the samples are compressed with an almost constant value of D ; a linear relationship between volumetric strain and D can then be observed above 8% axial strain under different confining pressures, which corresponds with the maximum sample dilation angle as shown in Figure 3b. A quadratic curve relation between the volumetric strain and fractal dimension can be obtained for different levels of axial strain, as shown in Figure  15. At 20% axial strain, the values of a, b and c are 1.48 × 10 −3 , 5.25 × 10 −4 and 0.26, respectively.

Conclusions
A DEM model of irregular-shaped ballast particles using the fractal replacement method (FRM) and a FDM flexible membrane model using shell elements are presented in this paper. A series of DEM-FDM coupled simulations of triaxial testing under monotonic loading and with different confining pressures were performed to study the ballast breakage and fractal behavior of PSD. We draw the following main conclusions: (1) The proposed FRM model can simulate the ballast breakage and capture the essential features of the stress-dilatancy in the monotonic triaxial test. During the loading process, the particles within the top and bottom half moved toward the horizontal plane in the middle. With axial strains above 4%, due to mutual resistance, the intermediate particles started to move outward to produce a shear sliding surface, where the angle between the sliding surface and the horizontal direction gradually increased. (2) More particle breakage occurred under higher confining pressure or larger axial strain and most such breakage was in the second generation or earlier. The number of fragments could be ignored at the initial compaction stage and the majority of particle breakage occurred with axial strain of 8% or above. Compared with the experimental results, the simulation underestimated the degree of particle breakage. This could be explained by that the simple splitting model of FRM did not take particle abrasion and small corner breakage into consideration. (3) With increasing axial strain or confining pressure, the percentage of small particle fragments increased in correspondence with the PSD curves which remained concave upwards, leading to an increasing fractal dimension. However, a different situation applies to granular materials in natural state because the grading ballast used in railway is usually qualified by the ballast specification, leading to much smaller values of D.

Conclusions
A DEM model of irregular-shaped ballast particles using the fractal replacement method (FRM) and a FDM flexible membrane model using shell elements are presented in this paper. A series of DEM-FDM coupled simulations of triaxial testing under monotonic loading and with different confining pressures were performed to study the ballast breakage and fractal behavior of PSD. We draw the following main conclusions: (1) The proposed FRM model can simulate the ballast breakage and capture the essential features of the stress-dilatancy in the monotonic triaxial test. During the loading process, the particles within the top and bottom half moved toward the horizontal plane in the middle. With axial strains above 4%, due to mutual resistance, the intermediate particles started to move outward to produce a shear sliding surface, where the angle between the sliding surface and the horizontal direction gradually increased. (2) More particle breakage occurred under higher confining pressure or larger axial strain and most such breakage was in the second generation or earlier. The number of fragments could be ignored at the initial compaction stage and the majority of particle breakage occurred with axial strain of 8% or above. Compared with the experimental results, the simulation underestimated the degree of particle breakage. This could be explained by that the simple splitting model of FRM did not take particle abrasion and small corner breakage into consideration. (3) With increasing axial strain or confining pressure, the percentage of small particle fragments increased in correspondence with the PSD curves which remained concave upwards, leading to an increasing fractal dimension. However, a different situation applies to granular materials in natural state because the grading ballast used in railway is usually qualified by the ballast specification, leading to much smaller values of D. (4) The evolution of fractal dimension against grain breakage showed a linear relation, Conversely, a quadratic curve relation was obtained between the fractal dimension and volumetric strain under different axial strain stages. Therefore, D has the potential to be a key indicator to evaluate the degree of ballast crushing and PSD degradation, which can contribute to better decision making concerning railway track bed maintenance in practice.