1. Introduction
Many of the world’s railways run on the ballasted track, which has for nearly 200 years provided stable support for train operation. However, with trafficking the geometry of the track deteriorates, mainly as a result of the development of differential settlement of the track-bed. When the geometry defects become too severe, maintenance is needed to realign the track to enable the continued safe running of trains [1]. Ballast failure is the reduction of the original quality due to various influences. The most significant factor contributing to the deterioration is the dynamic load. The dynamic load is directly related to the axle load and track geometry [2]. The factor that contributes to the degradation of railway infrastructure is usage (wear by physical contact and dynamic load), environmental conditions (climatic influence and water) and failures (faulty components, bad construction) [2]. And the process of track deterioration is mainly wearing and fatigue then final settlement [3].
Considerable maintenance work is needed for ballasted tracks due to the track deterioration. Therefore, it is very important to understand the mechanism of track deterioration and to predict the track settlement or track irregularity growth rate in order to reduce track maintenance costs and enable new track structures to be designed. Explanations on how the speed, load and repetition of traffic influence the long-term settlement of ballast in ballasted track are very scarce. Having in mind that tracks subjected to the same load show different settlement behaviors, explanations of track settlements in accordance with the speed, load and repetition are needed.
The study of the track structures and their response to applied traffic conditions are easier to study using software models. As the railway industry improves, the modeling of the railway track is becoming more realistic and representative of the real conditions. Track modeling is used to predict the performance of railway track components with the application of train movement and environmental conditions. Based on the complexity, railway track modeling techniques can be classified as two dimensional (2D), three dimensional (3D) and two point five dimensional (2.5D) [4]. From the very beginning, track has been modeled with the concept of Beam on Elastic Foundation (BOEF) which was developed by Winkler’s theory or Zimmermann method. These are categorized as two-dimensional analyses of a railway structure and are simplification of a beam (rail) laid on continuous support (soil’s subgrade or foundation) [4]. Through time the modeling technique becomes more realistic to three-dimensional finite element modeling (FEM) [5]. For granular materials like ballast and sub-ballast, discrete element modeling (DEM) is better representative [6].
This research sheds light on how traffic affects the long-term settlement of ballast and how the traffic parameters (speed, load and repetition) contribute in the process that has not been noticed by other researchers. This will be helpful in monitoring the conditions of existing railway tracks and designing new ones.
2. Materials and Methods
In this study, the effect of traffic was assessed based on an 8.08 m three-dimensional finite element track model loaded with a moving axel. Drucker Prager plastic model is applied to represent the granular behavior of ballast and sub-ballast. Parameters for the Drucker Prager plastic model are generated from a triaxial test model using PFC (3D) discrete element modeling software. On the model the sleeper spacing, rail gauge and characteristics of ballast material, sleeper and rail (steel) are adopted from CREC Addis Ababa-Mieso railway project, testing report for stone ballast No.1 and 2nd testing report for stone ballast [7], and Chinese Rail Standard (National Standards of the People’s Republic of China GB 2585-2007) [8].
2.1. Simulation of Triaxial Compression Test
Due to the need for plastic behavior of ballast material, for the study of the effect of traffic parameters on the settlement of ballast, triaxial compression test is chosen. This test is a promising way of evaluating the compression strength of a sample due to axial stress [9]. Axial and radial strains may be monitored during this test, to determine basic elastic constants (Young’s Modulus, E, and Poisson’s ratio, ν) [10]. Since it is more representative way to simulate the actual field stress conditions, many researchers used the simulation of triaxial test to assess the behavior of ballast system [11]. Parallel bonded clusters of spheres are created to simulate three types of ballast particle shapes [12] (bonding 13 and 14 pebbles to represent smaller particles and 28 pebbles for large sized particles) as shown on the Figure 1. This is achieved by generating clump templates of the three particle shapes and changing the pebble-pebble contact logic of clumps to parallel bonded pebbles to include breakage [13] in addition to the shape and angularity of ballast particles [14].
The particle breakage in PFC can be simulated by modeling particles as cluster of bonded spherical elements and breakage will be assumed to occur when the
![]()
(a) ![]()
(b) ![]()
(c)
Figure 1. Particle shapes used to represent ballast particles in the triaxial test simulation (a) With 14 pebbles; (b) With 13 pebbles; (c) With 28 pebbles.
bond strength between the elements exceeded by the internal forces due to applied load [15]. In this thesis, the contact model between clumped elements (called pebbles) is changed to parallel bond. With strength parameters mentioned below:
➢ Parallel bond normal and shear stiffness (kn and ks) = 6 × 1010 N/m.
➢ Parallel bond normal and shear strength = 5 × 106 N/m2.
The samples of the particles are inserted in a cylindrical shaped vessel and compressed to achieve a density of 1800 kg/m3. The test apparatus includes cylindrical vessel and two top and bottom horizontal plates [16]. Particle compaction and axial and confining pressure are applied by the movement of the plates and the cylindrical vessel [17]. The geometry of cylindrical vessel is composed of 20 triangular segments called facets. The top and bottom plates of the apparatus are composed of 8 segments. The test is performed under seven different confining pressures (30, 40, 80, 150, 200, 500, 1000 kPa). Since the size of the apparatus is set to 300 × 150 mm based on the H to D ratio of 2 (standard of triaxial compression test) [11]. The test apparatus and the particles assembly are shown on Figure 2.
![]()
Figure 2. The triaxial compression test simulation using PFC3D.
2.1.1. Virtual Triaxial Test Results
From this test, the relations of Deviatoric stress vs. Axial strain for different confining pressure (see Figure 3 for 40 kPa confinement) are used to plot p-t graph in Drucker Prager plastic model and then to calculate the friction angle (β) of the ballast material.
And relation between Volumetric strain vs. Axial strain (see Figure 4) is used to calculate the dilation (ψ) angle in linear Drucker Prager plastic model of ballast material for the FEM [18].
![]()
Figure 3. Deviatoric stress vs. axial strain plot from PFC3D simulation (40 kPa confinement).
![]()
Figure 4. Volumetric strain vs. Axial strain plot from PFC3D simulation (40 kPa confinement).
2.1.2. Validation of the Virtual Triaxial Test
Triaxial test is simulated using different DEM software (including PFC (3D)) in many former studies. Some of them also compared the result of this simulation with actual laboratory test results [19].
[20] has used PFC3D to simulate the triaxial test of sand material applying clumped particles. The results from this simulation have an error of about 0.46% to 5% (See Figure 5).
Clumps of PFC are used by [21] to investigate the behavior of fouled ballast material. The triaxial simulation at 3 different confining pressure resulted in precise results of volumetric strain vs. axial strain curve, with an error of about 1% - 5% (See Figure 6).
Based on the above discussion, the use of PFC3D software simulation for triaxial test of course grained material will give very precise “stress—strain” and “volumetric strain—axial strain” relation.
![]()
Figure 5. Deviatoric stress vs. axial strain diagram of DEM simulation and laboratory triaxial test [20].
![]()
Figure 6. Volumetric strain vs. axial strain curve of DEM simulation and experiment of triaxial test [21].
2.2. Finite Element Modeling (FEM) of Ballasted Railway Track
And the precise representation of these behaviors using software simulation is difficult. However, under essentially monotonic loading conditions rather simple constitutive models provide useful information. These constitutive models are essentially pressure-dependent plasticity models that have historically been popular in the geotechnical engineering field.
The model used in this research (The Extended Drucker-Prager model) is extension of the original Drucker-Prager model. In the context of geotechnical materials, the extensions of interest include the use of curved yield surfaces in the meridional plane, the use of noncircular yield surfaces in the deviatoric stress plane, the use of non-associated flow laws and creep [22].
Drucker Prager model along with its optimized models are successfully been employed to simulate the material response behavior of pressure-dependent soil, rocks and concrete [10]. As described on [23], the Extended Drucker Prager model in Abaqus is used to model frictional materials, which are typically granular-like soils and rock, and exhibit pressure-dependent yield (the material becomes stronger as the pressure increases). Ballast material is one of those.
The input parameters (dilation angle [ψ], friction angle [β], flow ratio [K]) for the linear Drucker-Prager analysis are calculated based on the triaxial test simulation results. These parameters’ description and their respective calculation for this thesis are discussed in the following sub sections.
The stress invariants of the model can be written in terms triaxial compression principal stresses as [22]:
(1)
(2)
(3)
where:
and
are the principal axial and lateral stresses respectively. Both are negative (compression).
The triaxial compression results can, thus, be plotted in the meridional plane.
Based on the equations above and rules the three input parameters of ABAQUS Linear Drucker-Prager model are calculated as follows.
Friction angle (β)
This parameter is calculated by plotting and fitting the best straight line through the t-p drawn from triaxial compression results [22].
According to [10], for triaxial compressiont = q. Therefore, using Equation (1) and Equation (2), t-p plot is generated from triaxial compression stress-strain plots of different confining pressure and best straight line fitted on the points as shown in Figure 7.
Then the value of β is calculated from the equation of the fitted line and is equal to 54˚.
![]()
Figure 7. t-p plot generated from the triaxial compression test simulation results.
Dilation angle (ψ): It is the volume change observed in granular materials when they are subjected to shear deformations. This effect was first described scientifically by Osborne Reynolds in 1885/1886 and is also known as Reynolds dilatancy [18].
Unlike most other solid materials, the tendency of compacted granular material is to dilate (expand in volume) as it is sheared. This occurs because the grains in a compacted state are interlocking and therefore do not have the freedom to move around one another. When stressed, a lever motion occurs between neighboring grains, which produces a bulk expansion of the material. On the other hand, when a granular material starts in a very loose state it may initially compact instead of dilating under shear. A sample of a material is called dilative if its volume increases with increasing shear and contractive if the volume decreases with increasing shear [18].
This parameter is calculated by using the formula recommended by [18], from the plot of volumetric strain-axial strain.
(4)
where: ψ is dilation angle in degrees
changes in volumetric strain;
changes in axial strain.
The value of the dilation angle of the ballast are calculated for each confining pressure (30, 40, 80, 150, 200, 500, 1000 kPa) and are equal to 12.84˚, 14.22˚, 15.9˚, 16.2˚, 16.24˚, 16.36˚, 16.7˚ respectively. As written in different researches the confining pressure of real ballast material is 30 to 40 kPa. So, the dilation angle for 40 kPa confining pressure which is equal to 14.22˚ is chosen.
Flow stress ratio (K) is the ratio of the yield stress in triaxial tension to the yield stress in triaxial compression. To ensure that the yield surface of a Drucker-Prager model remains convex requires 0.778 ≤ K ≤ 1.0.
This parameter is estimated in range to be 0.8, because the compression strength of ballast is much greater than its tensile strength.
Elastic parameters: Modulus of Elasticity (E) = 320 Mpa and Poisson’s ratio (ν) = 0.3 calculated from the test results under 40 kPa confining pressure.
Track Model Using ABAQUS
After the components of the track are modeled (see Figure 8 and Figure 9), twelve 8.08 m long three-dimensional track models are created (i.e., four for each three parameters of traffic load: speed, load and number of repetition). For speed parameter four models are created making number of pass (2 cycles) and load (15 tons) constant; with four speed variations of 80 km/hr., 100 km/hr., 120 km/hr. and 160 km/hr. And four models with varying load of 15, 20, 25 and 30 tons, making number of pass (2 cycles) and speed (80 km/hr.) constant. Four models also made varying number of passes of 2, 4, 6, and 8 repetitions, with constant speed (80 km/hr.) and load (15 tons). To assess the long-term behavior of ballast material due to cyclic traffic load, the model is subjected to an axel load of 25 tons with 120 km/hr. speed for about 8000 cycles. The properties of components of the track are summarized in Table 1.
![]()
Table 1. Summary of properties of track components.
• Step module
In this module two steps are created one for the application of the axel load. Its total time is 0.1 sec, and the other for the cyclic movement of the axel. Both the steps Dynamic-Explicit to avoid the convergence problem.
• Constraints
Kinematic coupling of the middle rod with wheels is performed defining reference point at the center of axel (see Figure 10).
Two contact interaction properties are created. One is for the wheel-rail contact and the other for other contacts in the model.
Wheel-rail contact: Frictionless tangential behavior and Hertz normal contact behavior is created based on [25].
Other contacts: a contact behavior with penalty frictional tangential behavior of friction coefficient equal to 0.7 and normal contact pressure-over closure relationship of linear type with stiffness of 400,000 N/m is created for the other contacts in the model.
• Springs/Dashpots
The rail pads are represented by spring/dashpot elements connecting the sleepers and the rails. Three spring elements with spring stiffness of 1.2 × 108 N/m2 and dashpot coefficient of 75,000 Ns/m2 are created per sleeper.
• Load module
In this module, two types of loads are created. The first is a concentrated axel load which is applied at the reference point created in the previous part. And the other is uniformly distributed (pressure) which is the sum of the dead load of the components of the track is applied on ballast, sub ballast and subgrade. The loads are applied based on the amplitudes defined.
• Boundary conditions (BC)
Three types of boundary conditions are created. The first one is a velocity/angular velocity type boundary condition which is applied to the axel to constrain the axel from any movement during loading step and to allow it to move cyclically with a specified velocity at amplitude created during movement step (see Figure 11). And the second is to constrain the lateral and longitudinal
![]()
Figure 11. Constant amplitude (Amp 2) used for cyclic movement of axel.
movement of track components. In the third BC the bottom of the subgrade is constrained from all displacements and rotations.
• Amplitudes
Amplitudes are created to define the repetition of load or boundary conditions in cyclic actions. Here two amplitudes are created; one is for the axel load (Amp 1) and the other is for cyclic movement of the loaded axel on the track (Amp 2).
• Mesh
The components assembled above are meshed separately based on their size to a total no of mesh elements equal to 9111 (See Figure 12).
3. Results and Discussion
The results used to meet the objectives of the research are stress-strain relationship of ballast, strain-time relationship of ballast, and stress-time relationship of subgrade and contours showing the distribution of stress and strain (See Figure 13 and Figure 14).
3.1. The Effect of Traffic Parameters on Ballast Strain
The cumulative plastic strain of ballast for each model is plotted against each corresponding speed, axel load and a number of cycles.
![]()
Figure 12. Meshed ballasted railway track model.
![]()
Figure 13. Stress contour for ballast, sub-ballast and subgrade.
![]()
Figure 14. Strain contour for ballast, sub-ballast and subgrade.
As it can be seen from the above plots (see Figures 15-17) the speed-strain relation has a relatively steepest slope at the beginning. The increase in the axel speed by 20 km/hr. increases the cumulative plastic strain of the ballast drastically in the first stage and gradually the slope became gentle; this is because in the first stages, the strain increases due to voids between particles but, gradually the amount of void between the particles will decrease and the strain becomes smaller. And on the speed-strain plot, significant amounts of total plastic strain are noticed. That means the variation of axel speed on ballasted railway track affects the total plastic strain of ballast material more than the variation of a number of cycles and axel loads.
![]()
Figure 15. Axel speed vs. total plastic strain.
![]()
Figure 16. Axel load vs. total plastic strain.
![]()
Figure 17. Number of cycles vs. total plastic strain.
3.2. The Effect of Traffic Parameters on the Stress Transfer to Subgrade
One of the most important performance criteria of ballast material is to resist load and to distribute it to wider area to minimize the stress transferred to the subgrade [26]. Here the stress transferred to the subgrade vs. time graphs of each model were extracted from the software. The maximum stress values are plotted against the corresponding speed, load and number of cycles. (See Figures 18-20).
From the above plots, it can be easily noticed that the variation of speed affects the load transferred to subgrade more than the two parameters.
![]()
Figure 18. Axel speed vs. maximum stress transferred to subgrade.
![]()
Figure 19. Maximum stress on subgrade vs. axel load.
![]()
Figure 20. Maximum stress on subgrade vs. number of cycles.
3.3. Strain-Time Relationship of Ballast
To assess the Strain-time relation of ballast, 25 tons of axel load moving with a speed of 120 km/hr. is applied for about 8000 cycles. According to the result, in the first about 25 to 30 cycles the ballast shows a rapid increasing of plastic strain that is due to reduction of depth as particles begin to interlock (See Figure 21). Then up to cycle, about 1200 the strain growth rate fluctuates (almost zero strain and then rapid growth). This is because of the particle resistance to stress due to enough interlock, and then particle breakage. In the next cycles up to 5000 plastic strain grows slowly due to particle final remaining breakage and compaction. Then after the strain rate decreases to have an almost constant value. That shows the ballast particles broken and well interlocked again. In this stage, the ballast’s drainage and damping performance will be lost since its particles are broken and it is compacted enough.
![]()
Figure 21. Plastic strain-time plot of ballast.
3.4. Stress-Time Relationship of Subgrade
In a track structure, a layer on the top is stronger than the one below it; so, it is the responsibility of the above layer to protect the layers under it. Here in ballasted track, the increase in stress on the subgrade is due to the degradation of the performance of ballast and sub-ballast to distribute load to wide area. As the thickness of ballast and damping performance decreases the pressure to be transferred to subgrade becomes higher. According to the finite element analysis result (see Figure 22), the value of stress on the subgrade fluctuates between the values 0.5 MPa and 2 MPa and tends to converge to a constant value around 1.5 MPa. This is because of the ballast settlement and reduced damping performance.
Validation of the finite element modeling
In many former studies, similar finite element modeling of railway track and simulation of train-track interactions are performed successfully. On some of them, it is managed to compare the results with experiments. [24] compared the stress transfer in similar track model using Abaqus, with a full-scale half-track structure experiment done by Jerry Rose from University of Kentucky. According to study stress transferred to the layers from the Abaqus model gives a precise result to the experiment with a difference of about 2% to 5% (See Figure 23).
![]()
Figure 23. Comparison stress between experiment and analytical model [24].
4. Conclusions
The effects of train speed and a number of cycles of train movement on ballast settlement are higher than that of the axel load. According to the results the increase in the speed of train movement by 20 km/hr. will increase the stress transferred to the subgrade by up to about 1000 kPa. Increasing the number of cycles of train movement by four affects the total plastic strain of ballast more than increasing of tonnage by 5 tons.
Ballast maintenance is mainly applied to correct two major deteriorations (failure) of ballast material due to cyclic traffic load and environmental conditions. These two failures are differential settlement of the ballast layer and the fouling or contamination of ballast material with fine materials. Tamping or stone blowing is the maintenance action that can be taken to correct differential settlement and another geometrical failure whereas ballast fouling is maintained by ballast cleaning action. Usually, ballast differential settlement happens before fouling because the main cause of ballast fouling is particle breakage and wear. According to the results and discussions in the above chapters, the performance reduction due to much settlement is reached after about 6000 cycles of moving axel load. Therefore, periodic monitoring of the performance of ballast material is recommended after 6000 cycles of axel load.
Acknowledgements
The authors would like to thank the Ethiopian Railway Corporation (ERC) and Addis Ababa Institute of Technology (AAiT) for the M.Sc. program and financial support given. Also owe their gratitude to the contractor of the Ethio-Djibouti railway line, China Railway Group Limited (CREC) especially Mr. Chen for their cooperation in giving valuable test data.