Advancements in Geo-Inclusions for Ballasted Track: Constitutive Modelling and Numerical Analysis

This paper reviewed some salient features evolving through mathematical and numerical modelling of ballasted track components incorporating recycled rubber products. Firstly, a constitutive model based on the bounding surface concept was introduced to simulate the shear stress-strain response of waste mixtures (i.e., recycled rubber crumbs, coal wash, and steel furnace slag) used for the capping layer placed below the ballast medium, whereby the energy absorbing property resulting from the inclusion of different amounts of rubber has been captured. Subsequently, key research findings concerning the inclusion of recycled rubber mats on ballasted tracks for reduced particle degradation under cyclic loading were examined and discussed. Discrete element modelling (DEM) coupled with Finite element modelling (FEM) to micro-mechanically characterise ballast behaviour with and without rubber mats offers invaluable insight into real-life track operations. In particular, this coupled DEM-FEM model facilitates the exploration of micromechanical aspects of particle breakage, contact force distributions within the granular assembly, and the orientation of contacts during cyclic loading.


Introduction
Rail tracks contribute to the largest worldwide transportation network. The increasing use of rail as one of the most sustainable and energy-efficient means of transport along with the associated increases in axle loads and train speeds have led to a rapid deterioration of track substructure components and exacerbated maintenance costs [1]. The need to maintain a competitive advantage implies there is pressure on the rail industry to improve its operational efficiency and passenger comfort, while minimising the maintenance and construction expenses [2,3]. Despite recent advances in the field of rail geotechnics, the rapid deterioration of track substructure (especially the ballast layer) is considered to be the most pronounced factor contributing to increased track maintenance costs [4,5]. To cope with this problem, there have been numerous studies into using artificial inclusions (rubber mats, geosynthetics, geocells, granulated rubber, etc.) to enhance rail ballast performance and reduce its degradation [6][7][8][9][10][11][12]. Field trials by Indraratna et al. [13] showed that rubber mats could significantly minimise ballast deformations, apart from reducing particle breakage. Previous laboratory studies using large-scale test facilities indicated that plastic under sleeper pads placed beneath sleepers could minimise ballast breakage by mitigating the transmission of impact loading effects [14].
In recent years, the Australian railway industry has trialled an innovation where waste and marginal materials are used as a sustainable replacement for the traditional capping layer. Energy absorbing recycled materials such as discarded rubber and ore wastes are blended into a new ductile material that will stabilise heavy haul tracks built on problematic terrain. The engineering properties of these new materials can be optimised such that they can be used during track construction instead of the traditional sub-ballast (capping) materials. Steel furnace slag (SFS) and coal wash (CW) are granular byproducts derived from steel manufacturing and coal mining, and since Australia is based on mining, it can produce several hundred million tons of SFS and CW per year [15,16]. In order to recycle these waste materials in civil engineering and curtail their adverse geotechnical properties (i.e., the swelling potential of SFS and breakage susceptibility of CW), engineers usually mix SFS with dredged materials, cement, fly ash, or asphalt for earthwork construction or pavement applications [17][18][19][20], while CW blended with SFS has been used very successfully for a Wollongong port reclamation project [21]. Not only are these recycled waste materials environmentally friendly and economically beneficial, but they can also help to solve some geotechnical problems. It has been found that adding rubber crumbs (RC) resulting from waste-car tyres to SFS-CW mixtures can significantly improve their energy absorbing capacity, which is associated with the high damping ratio of rubber products [12,[22][23][24][25][26]. This means that SFS+CW mixtures may also be used in projects involving dynamic loading by adding RC. Indraratna et al. [27] developed a selection system to optimise the SFS+CW+RC matrix for railway sub-ballast and found that a waste matrix consisting of SFS:CW = 7:3 and with the inclusion of 10% RC (by weight) is the optimal mixture for subballast materials. Qi and Indraratna [28] further confirmed this outcome via large-scale physical modelling and proved that using this optimal waste mixture (SFS+CW+RC) as subballast (or called capping layer which is the layer directly below ballast and above the subgrade, with a usual thickness of 150-200 mm) could efficiently absorb the energy and reduce ballast breakage. However, further investigation is needed to look into the mechanism of energy absorbing mixtures from a mathematical point of view.
Apart from traditional laboratory testing to understand the geotechnical properties of ballast and capping materials, mathematical and computational models are increasingly being demonstrated as alternative approaches with the possibility of implementing complex loadings subjected to varying boundary conditions [29][30][31][32]. Numerical simulations for real-scale track embankments subject to actual moving train loads (i.e., considering rest periods) could be more complex and time-consuming with either finite element or discrete element modelling approaches [33][34][35][36][37], so there is a need for co-simulations that can couple several computational schemes together to utilise computational resources and further advance railway research. This paper described mathematical and numerical approaches to investigate the improved performance of track substructure layers incorporating geo-inclusions (i.e., SFS+CW+RC mixtures and recycled rubber mats). A constitutive model adopting the bounding surface concept was developed for SFS+CW+RC mixtures, where the energy absorption capacity resulting from the RC inclusion is captured through an empirical relationship for the critical state parameter (M cs ) in relation to the total work input (W total ). Coupled discrete-continuum modelling was then used to examine the cyclic loading behaviour of a track specimen and the role of rubber mats in reducing ballast degradation. These mathematical and numerical models were calibrated and validated based on extensive laboratory testing.

Use of Waste Mixture for Subballast (Capping) Layer in Rail Tracks
Previous studies have examined the stress-strain behaviour of SFS+CW+RC blends subjected to monotonic loading [38] and found that the inclusion of rubber crumbs (RC) has a significant effect on the geotechnical properties of the mixture such that a higher RC content (R b ) makes the waste matrix more ductile and contractive, but with reduced shear strength. Qi et al. [39] found that with more than 20% RC, the SFS+CW+RC matrix could not attain a critical state, albeit there was a trend towards reaching a critical state with total axial strain over 25%, and thus extrapolation was adopted to estimate the critical state parameters for R b ≥ 20%.
In the following sections, the effect of the RC content on the critical state response is investigated, and the total work input is considered to be related to the critical state Sustainability 2021, 13,9048 3 of 20 parameters to capture the energy absorbing characteristics due to adding rubber. In this context, a constitutive model following the bounding surface theory and the critical state concept is then developed for the waste mixtures.
To calibrate the parameters of this constitutive model, the results of consolidated drained triaxial tests of the SFS+CW+RC matrix (SFS:CW = 7:3, R b = 0%, 10%, 20%, 30%, and 40%) from Indraratna et al. [27] were adopted. The waste matrix specimens were compacted to >95% of their maximum dry-unit weight, and before shearing they were saturated and then consolidated under varying effective confining pressures (σ 3 = 10, 40, or 70 kPa). The specimen and loading conditions for these tests were aimed at simulating the field conditions of subballast along the Australian coastline. The proposed model was then verified by comparing the prediction results with the laboratory data for SFS+CW+RC blends [27], as well as SFS+CW blends tested by Tasalloti et al. [40] and sand-RC blends examined by Youwai and Bergado [26].

Governing Equations
Conventional triaxial q − p notation was adopted in this study, where q = σ 1 − σ 3 denotes the deviator stress, p = (σ 1 + 2σ 3 )/3 is the mean effective stress, and σ 1 and σ 3 are the axial and radial stresses, respectively. Following the theory of bounding surface plasticity [41], the governing equation for the common stress-strain relationship can be expressed as: where m is the unit vector of plastic flow, n is the unit normal loading vector, H is the hardening modulus, ε v and ε q are the volumetric and deviator strains, respectively, and D e denotes the elastic compliance, which can be defined by: where K and G are the tangential bulk and shear modulus, respectively, as determined by: where κ refers to the gradient of the swelling line, e 0 refers to the initial void ratio, and v is the Poisson's ratio.

The Critical State
The critical state refers to a condition of a material where there is no further change in the mean effective stress and deviator stress with increasing axial strain, meanwhile, the dilatancy dε p v /dε q v is also approaching zero, i.e., dp dε q = dp Since rubber chips or tyre fibre continuously change shape towards the end of the triaxial test, it is difficult or even impossible to achieve a critical state. However, the rubber crumbs are usually smaller, so they are hindered from deformation by the surrounding stiff particles (e.g., SFS and CW). This implies that soil-RC mixtures can attain a critical state, albeit with larger amounts of RC the axial strain ε 1 can be higher than 25% [39]. Since most laboratory conditions are limited, ε 1 > 25% is difficult to achieve, and thus the parameters related to the critical state (e.g., the stress ratio at critical state M cs = q cs /p cs , and the critical void ratio e cs ) were obtained by extrapolation following Carrera et al. [42]. More details about the extrapolation for SFS+CW+RC mixtures can be found in Qi et al. [39]. The information associated with the critical state parameters for the materials investigated in this study is summarised in Table 1. The critical state line (CSL) for each SFS+CW+RC blend with certain R b (%) in the space of e − ln p is linear; this is identical to other granular materials and is represented by: e cs = Γ − λ ln p cs (6) where Γ is the void ratio at p cs = 1 kPa, and λ is the slope of the CSL in e − ln p space. By changing the amount of R b , the critical state parameters Γ and λ change with the amount of RC accordingly, and this effect can be captured by a linear relationship ( Figure 1a) as: where Γ 2 , λ 1 , and λ 2 are the calibration parameters computed on the basis of laboratory data for SFS+CW+RC mixtures (SFS:CW = 7:3 with different R b ) obtained by [39], and the test data for SFS75+CW25 (SFS:CW ≈ 7 : 3, R b = 0%) reported by Tasalloti et al. [40] is also shown in Figure 1a for validation purposes. The slight divergence between the proposed empirical relationship and the validation data may be due to the minor difference in the SFS:CW ratio. Similarly, a linear relationship is also found for sand-RC mixtures from [26], as shown in Figure 1b. This linear relationship is based on sand-RC mixtures that can reach or almost reach a critical state under laboratory conditions (test data provided in Table 1). Substituting Equations (7) and (8) into Equation (6), the CSL for SFS+CW+RC blends is now a function of p cs and R b , indicating a critical state surface can be formed in the space of e − ln p − R b , as illustrated in Figure 1c, that can be expressed by: where e * cs is the modified critical void ratio related to the RC content (R b ). Substituting Equations (7) and (8) into Equation (6), the CSL for SFS+CW+RC blends is now a function of ′ and , indicating a critical state surface can be formed in the space of − ln ′ − , as illustrated in Figure 1c, that can be expressed by: * = Γ + Γ − ( + ) ln ′ where * is the modified critical void ratio related to the RC content ( ). for each SFS+CW+RC blend with a certain amount of RC varies with ′ and , as seen in Table 1. It is also aware that the total work input is also influenced by the ′ and , which indicates that is the potential parameter that will reflect the changes of , and also help to incorporate the energy-absorbing property of the M cs for each SFS+CW+RC blend with a certain amount of RC varies with σ 3 and R b , as seen in Table 1. It is also aware that the total work input W total is also influenced by the σ 3 and R b , which indicates that W total is the potential parameter that will reflect the changes of M cs , and also help to incorporate the energy-absorbing property of the blends. This is in agreement with Chavez and Alonso [43], who also proposed that the work input is a good indicator to capture the variations of M cs that are influenced by parameters such as σ 3 and suction of rockfill. The increment of W total is the total incremental work input by the applied stress q and p , as shown below: where dW total is the total incremental work input and dε v and dε q are the incremental volumetric and deviator strains, respectively. Here, the total work input is calculated up to failure, which is taken as the peak deviator stress introduced by Zornberg et al. [44]. By connecting W total to M cs in Figure 2, an empirical relationship is established: where M * cs is the modified critical stress ratio by relating to the total work input, M 0 is the critical stress ratio when W total = 1 kPa, α is a regression coefficient, and W 0 = 1 kPa is used to keep the same units on both sides of the equation.
blends. This is in agreement with Chavez and Alonso [43], who also proposed that the work input is a good indicator to capture the variations of that are influenced by parameters such as ′ and suction of rockfill. The increment of is the total incremental work input by the applied stress and ′, as shown below: where is the total incremental work input and and are the incremental volumetric and deviator strains, respectively. Here, the total work input is calculated up to failure, which is taken as the peak deviator stress introduced by Zornberg et al. [44].
By connecting to in Figure 2, an empirical relationship is established: where * is the modified critical stress ratio by relating to the total work input, is the critical stress ratio when = 1 , is a regression coefficient, and = 1 is used to keep the same units on both sides of the equation. For validation purposes, data for SFS75+RC25 from Tasalloti et al. [40] are also shown in Figure 2, as is the good correlation with the proposed empirical relationship. Moreover, in Figure 2 a similar power relationship is established between and for sand-RC mixtures based on test data for mixtures involving 0%, 20%, and 30% RC [26]. This empirical equation will help to estimate the critical state parameters for soil-RC mixtures with larger quantities of RC, which cannot reach a critical state due to the limited laboratory conditions.

Bounding Surface and Loading Surface
This study used bounding surface plasticity theory because, unlike conventional elastic-plasticity theory, it has a versatile formulation for the loading and bounding surfaces, which eliminates the abrupt change between elastic and elastic-plastic response [45][46][47]. For validation purposes, data for SFS75+RC25 from Tasalloti et al. [40] are also shown in Figure 2, as is the good correlation with the proposed empirical relationship. Moreover, in Figure 2 a similar power relationship is established between M cs and W total for sand-RC mixtures based on test data for mixtures involving 0%, 20%, and 30% RC [26]. This empirical equation will help to estimate the critical state parameters for soil-RC mixtures with larger quantities of RC, which cannot reach a critical state due to the limited laboratory conditions.

Bounding Surface and Loading Surface
This study used bounding surface plasticity theory because, unlike conventional elastic-plasticity theory, it has a versatile formulation for the loading and bounding surfaces, which eliminates the abrupt change between elastic and elastic-plastic response [45][46][47].
The shape of the bounding surface was selected experimentally based on the undrained behaviour of granular materials. The bounding surface is shaped like a half water drop which aims to cover the triaxial compaction behaviour [48,49]. To facilitate further analysis, the same shape was adopted for the loading surface [45,50]; in this instance, the bounding Sustainability 2021, 13, 9048 7 of 20 surface F p , q, p c = 0 and the loading surface f (p , q, p C ) = 0 for the waste matrix are represented by [46]: where p c and p C refer to the interception points of the bounding surface and loading surface with the q = 0 axis, respectively, and they are the key factors controlling the size of the bounding surface and the loading surface ( Figure 3). N ≥ 1 is a constant relating to the material properties and it governs the curvature of the bounding surface.
is another material constant to express the ratio between p and p C , and the ratio between p and p C . According to the radial mapping rule [51,52], the stress ratio η is expressed as . By combining Equations (12) and (13), and then substituting the stress ratio η function, the material constant R can be rewritten as ility 2021, 13, x FOR PEER REVIEW 7 of 21 The shape of the bounding surface was selected experimentally based on the undrained behaviour of granular materials. The bounding surface is shaped like a half water drop which aims to cover the triaxial compaction behaviour [48,49]. To facilitate further analysis, the same shape was adopted for the loading surface [45,50]; in this instance, the bounding surface ( ̅ ′, , ̅ ′ ) = 0 and the loading surface ( ′, , ′ ) = 0 for the waste matrix are represented by [46]: where ̅ ′ and ′ refer to the interception points of the bounding surface and loading surface with the = 0 axis, respectively, and they are the key factors controlling the size of the bounding surface and the loading surface ( Figure 3). ≥ 1 is a constant relating to the material properties and it governs the curvature of the bounding surface. = = is another material constant to express the ratio between ′ and ′ , and the ratio between ′ and ′ . According to the radial mapping rule [51,52], the stress ratio is expressed as = = . By combining Equations (12) and (13), and then substituting the stress ratio function, the material constant R can be rewritten as * = − 1 ( * ) (14) The value of ̅ ′ controls the size change of the bounding surface which is connected to the change of the volumetric strain; therefore, to locate the position of ̅ ′ , the swelling line is combined with Equations (6)- (8) and Equation (14), and then the position of p c is obtained where e κ0 refers to the void ratio when p = 1 in Equation (15), p r denotes the unit pressure, and Γ(R b ) and λ(R b ) are critical state parameters modified by the amount of RC. The unit normal loading vector n at the image point on the bounding surface is then obtained through Equation (17): where σ represents the effective stress on the bounding surface, and n p and n q are the loading direction vectors.

Plastic Potential
The plastic potential, also regarded as the dilatancy relationship d = dε p v dε p q , is related to the state of soil corresponding to its compaction status and the applied pressure. Been and Jefferies [53] proposed a state parameter ψ after Worth and Bassett [54] to include the influence of the soil density and applied pressure on the deformation behaviour of soil, and ψ can be defined by Equation (18) as the difference between the current void ratio and the critical void ratio under the same loading conditions: By substituting Equation (9), ψ can be re-written as: According to Li and Dafalias [55], the dilatancy (d) is related to ψ, and may be represented by: where g denotes the plastic potential, and d 0 and m are constants related to the material properties. The plastic potential g = 0 can then be obtained by integrating Equation (20). The unit vector of plastic flow (m) at σ (the effective stress on the loading surface) can generally be expressed as: where m p and m q refer to the plastic flow direction vectors.

Hardening Law
According to the bounding surface concept, H (the hardening modulus) is expressed with two components: where H b refers to the plastic modulus atσ , and H δ is the arbitrary modulus at σ . H b is then attained by adopting an isotropic hardening rule associated with ε v p by using the following equation: Based on the bounding surface theory, H δ represents a decreasing function related to the distance between σ andσ [56], and by taking an arbitrary form it can be expressed as: where h 0 refers to a scaling parameter that controls the steepness of deformation behaviour in the ε v − ε q space, so it may be considered as a material constant. δ max is the distance between the stress origin and the image stress point, and δ is the distance between the current stress point to the image stress point (Figure 3). Following the radial mapping rule, remains positive, H δ is always positive, and only when p c = p C , H δ reaches zero, at which H = H b , and in this case the material is purely plastic. When δ max ≤ δ, H δ = +∞, H tends to be very large and the soil behaviour becomes purely elastic, but when H = −H b , H = 0, the soil sample changes from strain hardening to strain softening.

Parameters Calibration
This proposed constitutive model needs 12 parameters; here, there are two elastic parameters, i.e., Poisson's ratio v which is a constant related to materials properties, and κ which can be attained by unloading and reloading in isotropic compression.
The parameters reflecting the critical state are M 0 , α, Γ 1 , Γ 2 , λ 1 and λ 2 . α and M 0 are attained based on the relationship between the work input and critical stress ratio shown in Figure 2, and Γ 1 , Γ 2 , λ 1 , and λ 2 are calculated by curve fitting (Figure 1a,b). The parameter N is used to control the shape of the bounding surface and can be defined using a normalised q ∼ p plot of the yield points from undrained shearing behaviour. Due to the unavailability of the undrained test results for this study, N = 1 was set for simplicity, in light of previous studies [46,57,58] for granular soils.
d 0 and m are related to the soil dilatancy; m is obtained through Equation (20) at the changing point where the material turns from contraction to dilation; therefore, and the parameter d 0 can be determined by Equation (20) at the failure point which achieves the peak deviator stress [27], so d = d peak , ψ * = ψ * peak , and η = η peak . Therefore, h 0 is the parameter associated with the hardening modulus and it is obtained by fitting ε v − ε q curves. All the required parameters are listed in Tables 1 and 2.

Model Validation
In order to validate the proposed constitutive model, a comparison was established between the model predictions and the test results for the above-mentioned mixtures, including SFS+CW+RC mixtures [38] (Figure 4), SFS75+CW25 mixtures [40] (Figure 5a), and sand-RC mixtures [26] (Figure 5b). Note that the developed model captured the overall stress-strain and volumetric responses of the mixtures with adequate accuracy. Figure 4 shows that by increasing R b in the waste matrix, (i) the peak deviator stress decreases, (ii) the specimen becomes more ductile, (iii) the volumetric response is more contractive, and (iv) the behaviour changes from strain-softening to strain-hardening. Moreover, in view of Figures 4b and 5, when σ 3 increases, (i) the peak deviator stress increases, (ii) the strainhardening response becomes more pronounced, and (iii) the specimen tends to exhibit a more significant contractive behaviour.  Without exception, the model developed for the waste matrix has some limitations. For instance, the bounding surface is defined for > 0, which means that this bounding surface model can only predict compressive loading conditions, and the proposed relationship for relating to is specified for soil-RC mixtures with varying . It is not appropriate for soil-rubber mixtures with larger particles of rubber, such as chips or tyre fibre because they may not reach a critical state even after large axial strain (>25%) due to the continuous deformation experienced by larger rubber particles towards the end of the tests.

Discrete Element Modelling (DEM) of Ballasted Tracks Stabilised with Rubber Mats
The DEM was firstly proposed by Cundall and Strack [60] and is now widely used to explore the micromechanical response of ballast grains [5,34,[61][62][63]. The pioneering work in applying the DEM to simulate ballast was conducted by McDowell and Bolton [64] at the Cambridge University; it then continued with researchers at Nottingham University [65], the University of Illinois [61], the University of Wollongong Australia [66], and at the Southwest Jiaotong University [34], among others. It is noteworthy that the macroscopic strength and load-deformation responses of ballast are associated with the internal arrangement of grains (i.e., the granular fabric). Yimsiri and Soga [67] found that this granular fabric influences the shear behaviour of the granular assembly such that once the material shears in the same direction as the major contact normals, it becomes stiffer, stronger, and more dilative [64,68,69]; however, this material will not be as stiff or dilative if it is sheared in an orthogonal direction. Cui and O'Sullivan [70] investigated the microscale parameters of granular materials and confirmed there are significant stress and strain non-uniformities within specimens subjected to shearing.

Modelling of Ballast Grains in DEM
It is well known that the size, shape, and angularity of particles influence the mechanical properties and stress-strain responses of ballast grains [64,[71][72][73]. In this study, a 3D laser scanner was adopted to capture the actual shape, size and angularity of ballast. A number of representative ballast particles were chosen and scanned through the laser lights to form the polygonal-mesh, as shown in Figure 6a. Sub-routines were then programmed to build ballast particles in DEM by clumping many spheres together to fill up the polygonal-mesh (i.e., clump), as illustrated in Figure 6b. The parameters adopted for the DEM simulation of ballast were calibrated using the results obtained from laboratory tests for similar types of ballast grains. The input parameters were modified accordingly Without exception, the model developed for the waste matrix has some limitations. For instance, the bounding surface is defined for q > 0, which means that this bounding surface model can only predict compressive loading conditions, and the proposed relationship for M cs relating to W total is specified for soil-RC mixtures with varying M cs . It is not appropriate for soil-rubber mixtures with larger particles of rubber, such as chips or tyre fibre because they may not reach a critical state even after large axial strain (>25%) due to the continuous deformation experienced by larger rubber particles towards the end of the tests.

Discrete Element Modelling (DEM) of Ballasted Tracks Stabilised with Rubber Mats
The DEM was firstly proposed by Cundall and Strack [60] and is now widely used to explore the micromechanical response of ballast grains [5,34,[61][62][63]. The pioneering work in applying the DEM to simulate ballast was conducted by McDowell and Bolton [64] at the Cambridge University; it then continued with researchers at Nottingham University [65], the University of Illinois [61], the University of Wollongong Australia [66], and at the Southwest Jiaotong University [34], among others. It is noteworthy that the macroscopic strength and load-deformation responses of ballast are associated with the internal arrangement of grains (i.e., the granular fabric). Yimsiri and Soga [67] found that this granular fabric influences the shear behaviour of the granular assembly such that once the material shears in the same direction as the major contact normals, it becomes stiffer, stronger, and more dilative [64,68,69]; however, this material will not be as stiff or dilative if it is sheared in an orthogonal direction. Cui and O'Sullivan [70] investigated the microscale parameters of granular materials and confirmed there are significant stress and strain non-uniformities within specimens subjected to shearing.

Modelling of Ballast Grains in DEM
It is well known that the size, shape, and angularity of particles influence the mechanical properties and stress-strain responses of ballast grains [64,[71][72][73]. In this study, a 3D laser scanner was adopted to capture the actual shape, size and angularity of ballast. A number of representative ballast particles were chosen and scanned through the laser lights to form the polygonal-mesh, as shown in Figure 6a. Sub-routines were then programmed to build ballast particles in DEM by clumping many spheres together to fill up the polygonal-mesh (i.e., clump), as illustrated in Figure 6b. The parameters adopted for the DEM simulation of ballast were calibrated using the results obtained from laboratory tests for similar types of ballast grains. The input parameters were modified accordingly until the shear stress-strain behaviour predicted from the simulations agreed adequately with the experimental results.

Coupled Discrete-Continuum Modelling (DEM-FEM) for Ballasted Tracks with Rubber Mats
Numerical modelling via the coupled DEM-FEM approach enables an insightful investigation of the mechanical response of ballast and how rubber mats reduce ballast breakage under cyclic loading, as also observed from laboratory tests using a large-scale cyclic triaxial equipment (Figure 7a) [74]. Ballast particles are modelled by bonding many disks with predetermined sizes together (Figure 7b), and ballast breakage is assumed to occur when those bonds are broken. Figure 7c illustrates a diagram of a combined DEM-FEM model of the whole test specimen, where the discrete grains are modelled in DEM and the layers of subgrade soil, capping, and rubber mat are modelled following the continuum approach. The thickness of the ballast layer was 300mm, and is typically used in Australian railway lines. The rubber mat was modelled as a linear-elastic material and the capping was simulated as an elastoplastic material adopting the Drucker-Prager model. The rubber mat used in this study was locally made by an Australian recycling company that has a thickness of 10mm. The performance of this recycled rubber mat was justified by previous laboratory tests [10]. The dynamic bedding modulus of the rubber

Coupled Discrete-Continuum Modelling (DEM-FEM) for Ballasted Tracks with Rubber Mats
Numerical modelling via the coupled DEM-FEM approach enables an insightful investigation of the mechanical response of ballast and how rubber mats reduce ballast breakage under cyclic loading, as also observed from laboratory tests using a large-scale cyclic triaxial equipment (Figure 7a) [74]. Ballast particles are modelled by bonding many disks with predetermined sizes together (Figure 7b), and ballast breakage is assumed to occur when those bonds are broken. Figure 7c illustrates a diagram of a combined DEM-FEM model of the whole test specimen, where the discrete grains are modelled in DEM and the layers of subgrade soil, capping, and rubber mat are modelled following the continuum approach. The thickness of the ballast layer was 300mm, and is typically used in Australian railway lines. The rubber mat was modelled as a linear-elastic material and the capping was simulated as an elastoplastic material adopting the Drucker-Prager model. The rubber mat used in this study was locally made by an Australian recycling company that has a thickness of 10mm. The performance of this recycled rubber mat was justified by previous laboratory tests [10]. The dynamic bedding modulus of the rubber mat used in the test was C dyn = 0.107 (N/mm 3 ). Material properties were determined from the laboratory tests on materials sourced from a local quarry [75]. It is noted that there is only one type of ballast (fresh basalt quarried in Bombo, NSW) following the Australian standard AS 2758. 7-2015 [76] and was used in laboratory tests; therefore the effect of changes in ballast type on the performance of rubber mat was not studied. A series of large-scale laboratory tests were carried out under a stress-controlled mode where the maximum cyclic stress applied on the ballast sample was 230 kPa, which is equivalent to 16.25 kN vertical force applied by a dynamic actuator. It is noted that this force may not be exactly the same as the wheel-rail-sleeper-ballast interaction, as it was determined to generate an equivalent vertical stress in the model test. The laboratory test carried out in the laboratory was only considered vertical forces where lateral forces that were caused by a wheel impact on a locally bent rail and longitudinal force (e.g., braking force) have not been considered. The subgrade is simulated by a conventional Mohr Coulomb model whose parameters were predetermined in the laboratory. Note that the concept of coupled discrete-continuum has been established earlier by many studies [30,36,77] and adopted in many engineering problems, including stone columns, cellular geostructures, soil-structure interaction and railway engineering. As the main aim of this section is to review recent research work carried out by the authors in the context of ballasted track stabilisation using rubber mats, only the associated results were considered herein. large-scale laboratory tests were carried out under a stress-controlled mode where the maximum cyclic stress applied on the ballast sample was 230 kPa, which is equivalent to 16.25 kN vertical force applied by a dynamic actuator. It is noted that this force may not be exactly the same as the wheel-rail-sleeper-ballast interaction, as it was determined to generate an equivalent vertical stress in the model test. The laboratory test carried out in the laboratory was only considered vertical forces where lateral forces that were caused by a wheel impact on a locally bent rail and longitudinal force (e.g., braking force) have not been considered. The subgrade is simulated by a conventional Mohr Coulomb model whose parameters were predetermined in the laboratory. Note that the concept of coupled discrete-continuum has been established earlier by many studies [30,36,77] and adopted in many engineering problems, including stone columns, cellular geostructures, soilstructure interaction and railway engineering. As the main aim of this section is to review recent research work carried out by the authors in the context of ballasted track stabilisation using rubber mats, only the associated results were considered herein. The coupled model was implemented by transferring the forces and displacements between the two domains. Figure 7d shows how this exchange of forces (Fn, Fs) and displacements ( ) at their interfaces is facilitated. This figure shows that at a certain time The coupled model was implemented by transferring the forces and displacements between the two domains. Figure 7d shows how this exchange of forces (F n , F s ) and displacements ( . X [E] i ) at their interfaces is facilitated. This figure shows that at a certain time ∆t, contact forces in the DEM (Zone 1) were transferred to Zone 2 (continuum domain), and the displacements were transferred back to the DEM. This exchange of force and displacement can be performed using a mathematical framework and interface elements to distribute the forces and moments to nodal points (Figure 7e and Appendix A). This DEM-FEM model was adopted to simulate cyclic tests under different load frequencies (f = 10-40 Hz) and throughout N = 10,000 loading cycles.

Particle Breakage
This coupled model can be used to predict the accumulated broken bonds (B r ) as the load cycles N increase under varying frequencies (f = 10-40 Hz), as shown in Figure 8a.
Here, B r increased as the frequency f increased; this trend agrees with the ballast breakage measured in the laboratory [78]. B r increased remarkably in the initial N = 5000 cycles, and remained nearly constant towards the end of the simulations. Placing a rubber mat under the ballast layer led to a decrease in B r , an observation that reinforces the energy-absorbing capacity of the rubber mat.
∆t, contact forces in the DEM (Zone 1) were transferred to Zone 2 (continuum domain), and the displacements were transferred back to the DEM. This exchange of force and displacement can be performed using a mathematical framework and interface elements to distribute the forces and moments to nodal points (Figure 7e and Appendix A). This DEM-FEM model was adopted to simulate cyclic tests under different load frequencies (f = 10-40 Hz) and throughout N = 10,000 loading cycles.

Particle Breakage
This coupled model can be used to predict the accumulated broken bonds (Br) as the load cycles N increase under varying frequencies (f = 10-40 Hz), as shown in Figure 8a. Here, Br increased as the frequency f increased; this trend agrees with the ballast breakage measured in the laboratory [78]. Br increased remarkably in the initial N = 5000 cycles, and remained nearly constant towards the end of the simulations. Placing a rubber mat under the ballast layer led to a decrease in Br, an observation that reinforces the energy-absorbing capacity of the rubber mat. Figure 8b presents snapshots of broken bonds obtained at given numbers of load cycles (N) for f = 20 Hz. The number of broken bonds increased with N, however it reduced when a rubber mat is provided. It is seen that many bonds break directly under the top loading plate within the first N = 1000 cycles. Zhang et al. [34] studied the dynamic response of ballast and observed that ballast grains approximately 200 mm beneath the sleepers are more likely to experience higher stress levels and therefore undergo more pronounced breakage.  Hz. The number of broken bonds increased with N, however it reduced when a rubber mat is provided. It is seen that many bonds break directly under the top loading plate within the first N = 1000 cycles. Zhang et al. [34] studied the dynamic response of ballast and observed that ballast grains approximately 200 mm beneath the sleepers are more likely to experience higher stress levels and therefore undergo more pronounced breakage. Figure 9 shows the changes of contact forces in ballast and the contours of predicted stress in the subgrade and capping layers captured at N = 100-10,000 cycles, and for f = 20 Hz. Note here that the imposed loads were transmitted to the ballast particles as chains of contact forces whose thickness is proportional to their magnitude. High contact forces were concentrated below the top-loading wall, and a large portion of the load was transferred vertically to the underlying layer of the subgrade. Increasing the number of load cycles N led to an increasing number of contacts, which may be related to the rearrangement and breakage of particles. The inclusion of recycled rubber mats increases the contact area between the ballast and capping layer, reduces the maximum stresses, and also helps to attenuate the load. This, in turn, attenuates excessive contact forces and reduces the deformation and breakage of the ballast grains, which in practice would result in extended track longevity with associated technical, financial, and environmental benefits.  Figure 9 shows the changes of contact forces in ballast and the contours of predicted stress in the subgrade and capping layers captured at N = 100-10,000 cycles, and for f = 20 Hz. Note here that the imposed loads were transmitted to the ballast particles as chains of contact forces whose thickness is proportional to their magnitude. High contact forces were concentrated below the top-loading wall, and a large portion of the load was transferred vertically to the underlying layer of the subgrade. Increasing the number of load cycles N led to an increasing number of contacts, which may be related to the re-arrangement and breakage of particles. The inclusion of recycled rubber mats increases the contact area between the ballast and capping layer, reduces the maximum stresses, and also helps to attenuate the load. This, in turn, attenuates excessive contact forces and reduces the deformation and breakage of the ballast grains, which in practice would result in extended track longevity with associated technical, financial, and environmental benefits.

Conclusions
This paper presented salient outcomes based on recent mathematical and computational modelling of rail track foundation layers (i.e., subballast and ballast) incorporating marginal and recycled products, such as steel furnace slag (SFS), coal wash (CW), recycled rubber crumbs (RC) and rubber mats. A bounding surface constitutive model was established to predict the stress-strain behaviour of a blended waste matrix (SFS+CW+RC) expected to be adopted as an alternative sustainable capping layer. Computational modelling was then implemented using the discrete element method (DEM) and the combined discrete-continuum modelling approach (i.e., coupled DEM-FEM). The effect of recycled rubber mats (energy absorbing layer) placed beneath a layer of ballast was simulated under different cyclic loading frequencies (f = 10-20 Hz). The main conclusions from this study are summarised below: • A bounding surface model for SFS+CW+RC mixtures as a form of potential capping layer was introduced and validated based on laboratory data. The energy-absorbing characteristics of the mixture with different rubber contents were captured through an empirical relationship between the critical stress ratio (M cs ) and the total work input (W total ). The proposed model was also validated by considering other rubber-soil mixtures tested during previous studies [26,40]. The predicted stress-strain behaviour of the mixture was comparable with the laboratory observations. The observed properties of this compacted waste product could be considered favourable as an alternative capping medium to be adopted in railways, in lieu of traditional (geologic) materials.

•
The changes of broken bonds (representing ballast breakage) and the re-distribution and re-orientation of contact forces to support the external loads were captured via the coupled DEM-FEM modelling. The predicted results were in satisfactory agreement with the large-scale laboratory data. The inclusion of rubber mats reduced the broken bonds and the magnitude of maximum contact forces, hence reducing ballast deformation and degradation.

•
The constitutive and numerical modelling facilitated the analysis and investigation of track performance incorporating resilient recycled rubber materials (i.e., rubber crumbs and rubber shock mats). The proposed critical state surface capturing the energy-absorbing characteristics of the waste blends including RC is suitable for and can be extended to relevant studies on rubber-soil mixtures. The established DEM-FEM modelling provides a novel approach for further studies to look into the combined micro and macro-mechanical response of ballasted tracks considering ballast degradation. From an environmental sustainability perspective, the use of recycled waste products in track construction and rehabilitation will contribute not only to reduce waste disposal volumes to landfills, but also to more effective use of non-renewable natural resources, thereby reducing the carbon footprint of future rail infrastructure projects.
of the CONSTRUCT-Instituto de I&D em Estruturas e Construções-funded by national funds through the FCT/MCTES (PIDDAC).

Data Availability Statement:
The data presented in this study are openly available in reference number [38,74].

Acknowledgments:
The authors wish to acknowledge the financial assistance provided by the Australian Research Council (ARC) Discovery Project (ARC-DP180101916), ARC Industry Linkage Project (LP200200915) and ARC Industry Transformation Training Centre for Advanced Rail Track Technologies (ARC-ITTC-Rail). The financial and technical assistance provided by various in-dustry partners over the years including Australasian Centre for Rail Innovation (ACRI), RMS and Sydney Trains, c/o Transport for NSW, SMEC, Bestech, Douglas Partners, ASMS, South 32, Ecoflex Australia, Bridgestone, and Tyre Crumbs Australia) is gratefully acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Interface elements only receive forces at their end nodes (F XB , F YB and F XA , F YA ,). Thus, it is required to transfer moments (M) and forces (F X , F Y ) from a grain to an element, as illustrated earlier in Figure 7e. By taking the equilibrium of forces, we can derive: The equilibrium of moment is derived as: Equations (A1) and (A2) can be extended by: in which, the parameter β is given by: Replacing the Equations (A6) and (A7) to Equation (A3) results in: Rearranging Equation (A8) gives: Equation (A9) can be used together with Equations (A6) and (A7) to exchange moment and forces between the DEM and FDM via user-defined subroutines developed in FISH language (ITASCA).