A numerical framework for modelling settlements of railway ballast layers

Permanent deformation in ballast layers is a major contributing factor to the railway track geometry deterioration. In spite of a considerable amount of research on understanding and predicting performance of ballast layers, accurately capturing their settlements remains a challenge. In order to contribute to solving this important issue, a new numerical method for predicting ballast settlements is presented in this paper. This method is based on the finite element (FE) method combined with a constitutive model that captures permanent deformation accumulation in unbound materials under cyclic loading. This allows predicting permanent deformations of large structures and at large number of load cycles in a computationally efficient manner. The developed constitutive model is validated based on triaxial test measurements over wide range of loading conditions. Stress state in ballast layers has been examined with a 3D FE model, for several embankment structures and traffic load magnitudes. The determined stress distributions and loading frequencies were used as an input of the constitutive model to evaluate permanent strains and settlements of ballast layer. The influence of embankment structural designs and traffic loading magnitudes on the ballast layers settlements is examined and the results obtained are compared with the existing empirical performance models.


Introduction
Granular ballast is a major component of the railway track substructure controlling largely the track geometry deterioration.Correct predictions of ballast layer requirements are crucial, as under dimensioning implies a risk for structural failure and large maintenance costs and over dimensioning implies unnecessary large costs and environmental impact.
As observed in field or lab acceleration tests, the settlement in ballasted tracks occurs in two phases: (1) the initial settlement with a high rate, which corresponds to the process of reducing and consolidating the gaps between ballast particles; (2) track deterioration, at a lower rate that is linearly proportional with respect to time or tonnage [1].Generally, there are several mechanisms involved in track degradation, mainly the track settlement, fatigue of track elements and rail wear [2].In that, deterioration of unbound granular layers, caused by particle rearrangement, breakdown, and abrasive wear, plays a great role.Performance predictions of unbound granular layers is mainly addressed through empirical approaches, due to the complexity of the problem.Typical examples of those models include: the Shenton [3,4], the Sato [5] and the ORE [6] models.While those empirical models allow capturing the track settlements in an average sense and in a computationally straightforward manner, they have two important limitations.Namely, the model's validity is limited to their calibration interval which precludes evaluation of settlements for new material types, new structural designs, and new traffic load levels.Furthermore, possibilities of using those models to gain insights regarding the relative contributions to ballast settlements of different influencing factors, i.e., load/ stress amplitudes, porosity, subgrade stiffness, etc. are somewhat limited.Those shortcomings may, at least partially, be addressed by establishing a performance prediction model based on fundamental granular mechanics principles.The intention of the present study is to contribute to this important topic.
The need for accurate predictions of ballast layers performance with respect to permanent deformation accumulation motivated a considerable amount of experimental, theoretical and numerical research, cf.e.g [7][8][9][10].Semi-empirical models are available in the literature for prediction of permanent deformation as function of number of cycles and of tri-axial and deviatoric stresses, cf.e.g.[10,11].To address the problem from a granular mechanics point of view, several numerical modelling approaches for evaluation of permanent deformation in unbound materials have recently been proposed in the literature relying on finite element (FEM) and discrete element (DEM) methods [8,12].
Numerical modelling based on the DEM or FEM is a promising tool for ballast layer performance prediction, as it allows incorporating fundamental material properties into the analysis.Over past three decades, a considerable amount of recent research has been focused on development of FE models of railway track structures at different levels of complexity, considering the ballast layers as linear or non-linear elastic material, cf.e.g.Nguyen et al. [13], Gallego et al. [14] and Kalliainen et al. [15].With respect to ballast layer settlements, Mohr-Coulomb plasticity and hardening soil models have been used to model explicitly permanent deformation accumulation in ballast layers, by e.g.Jing et al. [16] and Kalliainen et al. [15].In particular, as discussed by Kalliainen et al. [15] a hardening soil model allows to capture the effect of non-linear ballast properties on the load-bearing capacity of the railway track as well as on local stress and strain distributions in different structural layers.Computational results presented by Jing et al. [16] indicate that incorporation of Mohr-Coulomb plasticity allows accurately capturing the effect of steel slag on ballast layers' dynamic properties, in terms of their damping ratios and dominant frequencies.At the same time, non-linear FEM modelling of railway structures is associated with long computational times, which precludes its applications for analysis of larger structures and at longer time scales.The latter aspect is of particular importance with respect to ballast layer settlement prediction during the track service life.There are also numerous studies of prediction of railway embankment settlements using DEM presented in the literature.Among others, Huang and Tutumluer [17] studied the effect of fouling, and Liu et al. [18] studied how the dynamic behaviour of the track is affected by ballast degradation and investigation on how the stress state is affecting the settlements was presented by Bian et al. [19].While DEM certainly is a suitable tool for studying these effects in detail as it directly provides a link between the ballast particles and the settlements, the studies are rather limited to the studied geometrical size and time duration due to the computational expense of DEM simulations.
In a recent study, Grossoni et al. [20] proposed a simplified semianalytical approach allowing to evaluate the progression of permanent deformation accumulation in railway track beds and their effect on the dynamic loads arising at vehicle-track interaction.In their approach, short-term analysis of train-track interaction during one train passage has been combined with a semi-analytical model to predict ballast layer settlements over number of cycles.The semi-analytical model for ballast layer settlements proposed by the authors calculates the permanent strain accumulation of the ballast layer as a function of the vertical pressure exerted on the ballast layer.As shown by Grossoni et al. [20], combining the short-term and long-term analysis allows to capture the long-term settlements in the ballast layer in a computationally efficient manner.At the same time, Grossoni et al. [20] study a relatively simple short-term analysis approach has been used to evaluate the loads applied to the ballast from slippers, with a railway track structure represented as an elastic beam on discrete supports, characterized by their dynamic stiffnesses and damping ratios.This simplified approach precludes evaluation of 3D effects, such as e.g.change of confinement level with depth in the ballast layer, angle of side slopes of the ballast layer as well as variation of the ballast layer properties in the lateral direction, on permanent deformation accumulation.Furthermore, the permanent deformation accumulation model proposed by Grossoni et al. [20] has not been evaluated fully with respect to capturing the permanent deformation behaviour of the granular materials over wide range of triaxial loading conditions and loading frequencies.The intention of the present research is to remedy these shortcomings, through proposing an alternative computational approach for ballast layer settlements, combining short term 3D analysis of the railway structure based on the FEM and new constitutive model capturing the accumulation of permanent deformation in ballast materials under cyclic loading.
The aim with the present research is to address the gap of having a computationally efficient tool to predict long term permanent deformation accumulation in ballast layer materials.An interesting approach, that serves as inspiration for this work, was presented by Suiker and de Brost [21].They used a rather advanced viscoplastic constitutive model to predict settlements using cyclic triaxial test data.However, even though their approach allowed to study a full track structure in 2D up to 50 000 load cycles, an even more computational efficient approach would be beneficial.To address this gap of having a computational efficient tool, a new modelling framework for unbound granular materials (UGM) is proposed, capable of capturing the accumulation of permanent deformation under cyclic loading in computationally efficient manner as the stress state only needs to be calculated at one representative load cycle.The model's capability to capture UGMs behaviour over wide range of loading conditions is evaluated based on triaxial test data from the literature [9].To investigate the feasibility of using the developed model for ballast layers settlements predictions, the stress state induced in ballast layers by traffic loads has been examined with a 3D finite element model of railway embankments.The determined stress distributions and representative loading frequencies were used as an input of the constitutive model to evaluate permanent strains and settlements induced in a ballast layer by traffic load cycles.The computational framework developed has been applied to perform a brief parametric study with respect to the influence of embankment structural designs and traffic loading magnitudes on the ballast layers settlements.The frameworks capability to capture the effect of those parameters on ballast settlements is evaluated and compared with the existing empirical performance prediction models.A clear future use of the framework is to investigate how different structural components, like crossings and switches, affect the local settlements.

Modelling of permanent deformations in unbound materials
A model capturing the accumulation of permanent deformation in unbound granular materials (UGMs) as a function of fundamental material properties and of the triaxial stress state acting on the material volume is a principal component of ballast layer settlement performance prediction.The reliance on fundamental material properties is of particular importance, for the performance predictions to be valid for new structural and material designs as well as at new traffic loads.The model and material parameter identification procedure for the permanent strains in UGMs are described in detail in sections 2.1 and 2.2 respectively.

Constitutive model for permanent deformation in UGM
Cyclic stress states cause different types of deformations in UGM depending on the relationship between the static and cyclic stress states.A sketch at which combinations of cyclic effective stress and static pressure these types of deformation occur is drawn in Fig. 1.For predicting settlements in ballast layers, the regime with cyclic permanent strains is the main interest and especially how this region changes with the number of load cycles.
The proposed model is based on the model developed by Suiker and de Borst [21] and the experimental work by Sun et al. [9].The model has a similar structure as the Perzyna model for viscoplasticity giving it a strong theoretical foundation.The main difference from the approach in [18] is that the outlined model is only used for calculating the permanent strain in a post-processing step of the, later discussed, FE results of the stress state in the railway embankment.This leads to a computationally efficient approach and makes it possible to use any suitable model for calculating the stress state in the ballast from one train passage.The derivate of the permanent strains with respect to loading cycles, N, is decomposed into a part accounting for frictional sliding of the stones and a volumetric part accounting for compaction where κ f and κ v are history parameters accounting for these types of deformations and will be discussed below, δ ij is the Kronecker delta and n f ij giving the directions of the permanent strains at frictional sliding.
Assuming that the deviatoric shape change is co-linear with the devia- where p c is the hydrostatic pressure of the cyclic part of the stress tensor and q is the cyclic von Mises stress.The factor f ( ) determines the volume change due to frictional sliding.At low pressures, a volume expansion is expected due to the rotation of the individual stones and at higher pressures a compaction behaviour is expected.It is also found in the experimental data that the dependency of the static pressure is nonlinear with a decreasing effect of pressure at increasing pressure.A simple model for this behaviour is where b 1 , b 2 and b 3 are positive material parameters and p s is the static hydrostatic pressure.In [21], the following expression for the evolution of k f was proposed where α f and γ f are material parameters, f f (σ ij ) is a function accounting for the severity of the stress state that will be discussed below and h f (κ f ) is a stabilizing term.The notation 〈x〉 is the Macaulay brackets , no further strains due to frictional sliding develop.In [21], f f ( σ ij ) = p/q but presently a model based on the KST-model [22,23] is used to fully capture the pressure dependency seen experimentally which was not achieved with initial attempts using the linear model originally proposed by Suiker and de Borst [21].The pressure dependency function reads where a 0 , a 1 and a 2 are material parameters.The expression for h f ( κ f ) same in the present work as in (Suiker and de Borst, 2003) with h 0 , h m and η f being material parameters.The parameter κ f,0 accounts for permanent deformations prior to the analysis like deformations induced during construction.With Eqs. ( 6) and ( 7), Eq. ( 4) takes the form where the new parameters A, A 1 , A 2 , H 1 and H 2 becomes The cyclic loading of the ballast material also leads to a volume compaction accounted by the state variable κ v in Eq. ( 1).This is modelled similar to the permanent strains due to frictional sliding.However, as will be shown later, these permanent strains are significantly smaller than the strains at frictional sliding, it is deemed that a linear function is sufficient.
where the use of Macaulay brackets results in only a positive compaction contribution.However, the volume expansion can still be negative due to volume changes from frictional sliding, the dilatation.A comment is that this type of permanent strain model also can account for mixed types of traffic if N instead of cycles is the number of days.Then, the right-hand side of Eqs.(8), and (10) will consist of several terms, one for each train passage.

Material model parameter identification
To describe the behaviour of UGM under cyclic loading at general stress states, many material parameters must be determined.In the constitutive model presented above, there are seven parameters to describe the evolution of deviatoric permanent strains due to frictional sliding, A f , A 1 , A 2 , H 0 , H 1 , η f and γ f .Furthermore, there are in total seven parameters for describing the volumetric strains, three, for capturing volumetric changes at frictional sliding b 1 , b 2 and b 3 and four, c 1 , c 2 , A v and γ v to describe the volumetric compaction.A parameter calibration is performed for each testing frequency.An attempt was made to find parameters that are common for all loading frequencies, but such model did not accurately enough describe the test data.Further development should be based on more fundamental studies involving movement of individual particles, preferably using the Discrete Element Method (DEM).However, this is beyond the scope of the present work.
The model for the development of permanent strains in the ballast material has been calibrated against the large triaxial dataset in [9] including different stress state and loading frequencies.The material was lattite basalt which is a common aggregate for railway ballast layers in Australia with gradation specified according to [24].A cylindrical sample was used having a diameter of 300 mm and a height of 600 mm.This size is sufficiently large to avoid that local packing configurations affect the response.
By calculating the deviatoric axial strain e 11 = ε ij − ε kk /3δ ij the parameters determining the frictional sliding behaviour can be determined.Numerical experimentation shows that the effect of H 1 is very small and can therefore be set to zero.Furthermore, it is assumed that no initial compaction has occurred and hence κ f,0 is set to zero.The results of this calibration are shown as dashed lines in Fig. 2.
It is seen that the model captures the experimental results well for all load cases and frequencies.The strain levels at the highest frequency, 40 Hz, are well predicted but the instability of the sample at deviatoric stresses of 276 kPa and 460 kPa is not fully captured and will occur at higher deviatoric stress levels in the model.However, this instability can be strongly dependent on the behaviour of individual aggregates leading to a large spread in the strains where this instability occurs.Also, the instability occurs at a combination of high frequency and high deviatoric stresses.It will be shown later when discussing the stress state in the embankment that this combination does not occur in the ballast during traffic loading.It should be noted that at each loading condition only one experiment has been performed and hence the scatter in the data is unknown.
The two contributions to the volumetric strains cannot be separated experimentally and hence all parameters governing the volumetric behaviour must be determined simultaneously.The outcome of this calibration is shown with dashed lines in Fig. 3 where it is seen that the model could describe the experiments well, including the dilatation at low confining pressures and loading frequencies.This is an clear advancement form the original model by Suiker and de Brost [21] where only a monotonic accumulation of permanent strain could be predicted.The parameters determined in all calibrations are presented in Appendix A.
The calibration of the model to different loading frequencies enables to investigate how the settlements are affected by the train speed.In [25] it was concluded that rest periods between the loads do not affect the evolution of permanent strains.Thus, is concluded that it is the loading rate, increasing with increasing frequency that has a major influence on the accumulation of the permanent deformation instead of the loading frequency itself.This can be further illustrated that it is the load rate that governs the dynamics of individual stones in the ballast and thus affect the local movements causing settlements.This stress rate is varying spatially in the embankment structure.However, for simplicity, in the present analysis a constant frequency will be assumed for the cyclic stresses in the ballast.The effect of train speeds can still be studied as the loading rate will be proportional to the loading frequency.This will be commented more upon in Section 4.
To further illustrate the predictability capabilities of the permanent strain model, the of the deviatoric axial strain at uniaxial conditions after 500 000 cycles as function of the effective stress is shown in Fig. 4 (a) for different static pressure and frequencies and compared with the Fig. 2. Comparison between the deviatoric axial strain obtained experimentally in [9] and the predictions from the calibrated model.experimental final strains from Fig. 2. A decent agreement with the experimental results is seen and the instabilities at 40 Hz is predicted at q = 300 kPa and q = 420kPa for pressures of 30 kPa and 60 kPa respectively which is close what is seen in the experiments.Furthermore, combinations of different deviatoric stresses and hydrostatic pressures needed to obtain a given deviatoric axial strain ε eff after 1 000 000 load cycles is illustrated in Fig. 4 (b).It will be shown later that pressures approximately between 5 kPa and 30 kPa and deviatoric stresses up to 150 kPa are the most relevant for settlement predictions.
To further show the applicability of the permanent strain model, the model has been calibrated to an additional dataset taken from literature.The dataset, from Tennakoon and Indraratna [26], present permanent deformation of ballast showing different degree of clay fouling.It uses the same sample size as in [9] with a loading frequency of 20 Hz, static pressure of 10 kPa and cyclic deviatoric stress of 230 kPa.It is evident from the results in Fig. 5 where the experimental results are compared with the calibrated model that the permanent strain model also can capture this dataset and thus be used for capturing different degree of ballast fouling.It can also be noted that this material is stronger than the material used for Figs.2-4 as the base material (VCI = 0 %) shows lower permanent strains than the corresponding data in Fig. 2 for 20 Hz.

Finite element modelling of railway embankments
The constitutive model developed in section 2 may be used to predict settlements in ballast layers, provided that the stress state induced by passing trains is known.Presently, the stresses induced in ballast layers are evaluated with the 3D FE model described in section 3.1.To gain insight into the approach capabilities with respect to capturing the effects of different structural designs and traffic load levels, several analysis cases are defined as follows.Models of embankments having different heights are created and the loading can either be applied using a continuous concrete slab or sleepers.The models have been loaded using three different axle load levels, 17.5 tonnes to simulate high-speed trains, 22.5 tonnes pertinent to normal track loading and 30 tonnes for heavy freight traffic.
The procedure used to calculate the settlements from the calculated permanent deformations is described in section 3.2.

Railway embankment model
The FE software Abaqus [27] has been used for calculating stress state in the embankment.Full 3D models of embankment structures are created, and they consist of layers with different densities and mechanical properties as follows.The embankment geometry and material properties of the layers used in this study are chosen based Swedish Transport administration norms [28].The analysed structures are however somewhat simplified as compared to the type sections, in Swedish Transport Administration guidelines.In particular, as the triaxial test data is only available for ballast material, the ballast and sub-ballast layers are modelled with the same constitutive model, both with respect to their elastic properties and their resistance to permanent deformations.Furthermore, for the cases of slab tracks, the slabs are placed on a thin layer of cement stabilized material.As the focus of this study is permanent deformation in unbound layers, the cementstabilized layer is not included explicitly in the model.The subgrade is modelled as an elastic material to obtain a realistic flexibility of the ground.The aim with including the subgrade is to provide a bottom layer with realistic flexibility for the stress calculation in the ballast.The rails and the concrete for ballast and slabs are also assumed to be elastic with standard values.Using a purely elastic material for the ballast results in significant tensile stresses in the ballast and to avoid this issue a Fig. 3. Comparison between the volumetric strain obtained experimentally in [9] and the predictions from the calibrated model.

E. Olsson et al.
tension cut-of 10 kPa was set being a compromise between having a realistic tension cut-off value and achieving numerical convergence in a reasonable time.The elastic modulus, E, the Poisson's ratio ν and the density ρ of each material is given in Table 1.The presented simulation framework allows for changing the constitutive model for calculating the stresses incorporating different non-linear effects in the ballast layer but for the present case, a simple model is chosen for brevity.
Due to symmetry as only straight embankments are modelled, only half of the embankment is modelled.A cross-section of a generic

Table 1
Material properties used for the different layers in the embankment model, for the concrete used for the slab and sleepers, and the steel used for the rails.embankment geometry is shown in Fig. 6 along with the boundary conditions used.The width of the layers beneath the ballast W, the height of the subgrade H subgrade and the length of the embankment are chosen so that any effect from the remote boundaries is negligible.With a length of 5.5 m together with symmetry conditions the studied track length is 11 m representing an infinitely long straight section.The displacement on all remote boundaries is set to zero.Two different geometries are studied to investigate the effect of different embankment heights.The models are named "High embankment" and "Low embankment" with dimensions given in Fig. 6.The loading is applied by explicitly modelling rails mounted on sleepers or on a concrete slab.The sleepers are modelled as rectangular blocks made of concrete with a dimensions of 2.5 m × 0.175 m × 0.265 m (LxHxW).The concrete slab extents over the whole length of the railway embankment and has a height of 0.5 m and a width of 2.5 m.The material parameters for concrete are given in Table 1.

Material E[MPa] ν[ − ] ρ[kg/m
The FE model has been discretized with second order brick elements, the maximum element length was set to 0.1 m for the ballast and the concrete slab, 0.5 m for the subgrade, 0.05 m for the rail and 0.025 m for the sleepers.Using this discretization, the "High embankment" and the "Low embankment" models consist of 126 830 and 94 695 elements respectively.Two examples of meshed models are shown in Fig. 7, the "High embankment" with sleepers in Fig. 7 (a) and (b) and the low embankment with a concrete slab in Fig. 7 (c) and (d).The mesh has been checked for convergence by having a coarser mesh which gave a negligible difference in the predicted stress state.
The model is loaded in two steps; firstly, by gravity to obtain the static stresses and thereafter by traffic load.Currently, only loading from one axle is studied and that load is applied as a pressure on the upper surface of the rail element closest to the symmetry plane in the length direction as shown in Fig. 7 (b).It is assumed that the accumulated deformations do not affect the geometry significantly and hence the stresses calculated at the initial geometry will be valid throughout the analysis.However, it is possible with this modelling approach to update the geometry using the calculated permanent deformations in a computational efficient manner by displacing the nodes using the nodal permanent deformations calculated in the next section followed by a new stress calculation.

Settlement calculations
The permanent strain and settlement calculations are performed as a post-processing step, using the static stresses and the cyclic stresses, defined as the difference between the stresses at the end of the loading and gravity steps.The permanent strain tensor is evaluated at all integration points of each element by integrating the rate equation Eq. ( 1) numerically.
The calculation of the permanent deformations from permanent strains is based on how strains are calculated from deformations in the FEM.Given nodal displacements in the x-y-and z-directions, u, v and w of an element with n nodes, u e = [u 1 , v 1 , w 1 , u 2 , v 2 , w 2 ⋯u n , v n , w n ] T , the strains within that element are calculated using the matrix B as ε = Bu e (11) where ε is a vector of the six strain components and the matrix B has six rows and 3n columns.The strains at all integration points in the whole model can thus be calculated simultaneously from the displacements from all nodes using ε G = B G u G (12) with ε G being a vector with all strain components in the model organised as follows 13) giving in total 6N e m number of strain components where N e is the number of elements and m the number of integration points per element.
The vector u G contains all 3n n nodal displacements u G = [u 1 , v 1 , w 1 , ⋯w Nn ] T with N n being the number of nodes in the model.The global strain-displacement matrix B G will have 6N e m rows and 3N n columns and will be sparse and its non-zero components are by the assembly algorithm in Algorithm 1.
Presently, the strains ε G are known and the permanent nodal dis- The boundary conditions used in the finite element model are used as boundary conditions in the evaluation of the permanent deformations and as only settlements of the ballast layers are considered, the vertical movement of the nodes connected to the subgrade layer is set to zero.The permanent deformations are determined by solving the system of equations in Equation ( 14) and each equation, pertinent to one strain component in a specific integration point, is weighted with the volume belonging to this integration point.One benefit with this method is that it will regularise the deformations as the combination of low confining pressures and large deviatoric stresses in some integration points could cause locally unrealistically large strains.

Algorithm 1: Calculation of the non-zero components of BG
For each element i = 1 : Ne: Collect the numbers of all nodes belonging to the element in a list nodes For each integration point j in element i, j = 1 : m: Calculate strain-displacement matrix B at integration point j Set l = 1 For each node number k in nodes:

Stresses in the ballast layer
Frist it seems appropriate to discuss the stresses states obtained from the FE models as they govern the final settlement calculations.From Eq. ( 6) it is given that the combination of low static pressures and large deviatoric stresses causes large permanent strains.An example of the static pressure field, for the model "High Embankment" with sleepers, is shown in Fig. 8 (a).At the edge of the sleepers, large pressures occur due to singularities in the stress field at contact edges.The cyclic von Mises stress for the same model loaded with 22.5 tonnes axle load is shown in Fig. 8 (b) where the largest von Mises stress, except for the singularity at the sleeper edge, are found slightly below the surface as expected from contact mechanics considerations, cf.[29].
A quantitative analysis of the stresses is found in Fig. 9.The static pressure is visualized Fig. 9 (a) along a vertical line under the rail, shown in Fig. 8, for the "High embankment" models.The results for the low embankment are almost identical.It is seen that the concrete slab creates slightly higher pressures than sleepers in the upper part of the embankment due to its weight.The difference between the slabs and sleepers is much more pronounced when studying the cyclic von Mises stress, Fig. 9 (b), as it becomes approximately five times higher in the case of sleepers than for a concrete slab.The results in Fig. 9 also indicate that the stress levels used for calibrating the model are relevant for the application.
To connect loading frequency and train speed, the rate of the cyclic stress is examined by considering the stress variation along the embankment.The loading rate at a given train velocity v can be calculated from and will be compared with the rates used in [9].The cyclic von Mises stress below the rail, at the depth where the maximum von Mises stress is found in Fig. 9 (b), is visualized in Fig. 10 (a) and in Fig. 10 (b) for the embankments with a concrete slab and sleepers respectively.Only results from the high embankment model are presented.In the case of a slab track, the increase in cyclic deviatoric stress occurs over approximately 5 m whereas in the case of sleepers the increase in stress occurs over approximately 2 m and the estimated stress gradients are presented in Table 3.As most of the cyclic deviatoric von Mises stresses used in [9] was 230 kPa, the relevant frequency can be estimated by comparing the loading rate in the finite element model and used in [9].It is seen that for the concrete slab that the relevant frequencies are well below the lowest frequency used for calibrating the model, 5 Hz and thus the model results for 5 Hz will be seen as an upper bound for settlements when using a concrete slab.For an embankment with sleepers, the relevant loading frequencies are between 5 Hz and 10 Hz depending on the axle load and the train velocity as seen in Table 3.A final remark is that the size of the region with highest cyclic stresses is similar to the size of the sample in the cyclic testing which limits the risk of size effects when calculating permanent strains using data obtained from the cyclic tests.,

Permanent deformations
The permanent strains, along the lines in Fig. 8 are shown in Fig. 11.The permanent strains are evaluated using the model outlined in section 2.1 with parameters determined by calibrating the model against the dataset from Sun et al [9].The volumetric strains, ε v = ε 11 + ε 33 + ε 33 , is drawn in Fig. 11 (a) and the effective deviatoric strain, defined in the same way as the von Mises stress, is depicted in Fig. 11 (b).It is seen that the volumetric strains have only a minor contribution to the total strain field.A large difference is seen as expected when comparing the deviatoric strains for embankments with slabs and sleepers.Deviatoric strains below 1 % are seen for all axle loads for the concrete slab whereas the region with strains over 1 % is substantial for the case of sleepers.
The permanent strains are integrated to permanent deformations according to the procedure in Section 3.2.The vertical permanent deformations at 22.5 tonnes axle load, loading frequency of 5 Hz and after 1 000 000 cycles, are shown in Fig. 12 (a) and (b) for models with concrete slab and sleepers respectively.As expected, a significant difference is seen between the two models with about 4 times larger permanent deformations in the case of sleepers.
The permanent vertical deformations along the lines in Fig. 8, at different number of load cycles and axle loads, are shown Fig. 13, for the loading frequency of 5 Hz.In addition to the permanent deformations' magnitudes, two significant differences are seen between the sleeper and slab models.Firstly, the accumulation of permanent deformation is approximately uniform over the height of the embankment in the model

Table 3
Estimated stress gradients dq/dx and from them estimated relevant triaxial test frequencies by comparing stress rate in the triaxial testing in [9] for the investigated load cases at different train velocities.with a concrete slab whereas in the case of sleepers the majority of the deformation is accumulated in a region extending 1 m below the top of the embankment.This is explained by the stress state visualized in Fig. 9 where large cyclic stresses occur in the same region.The second difference is that for the model with a concrete slab there is a small but continuous accumulation of the permanent deformation over the load cycles whereas for the sleeper model, most of the permanent deformation is accumulated the first 10 000 cycles.The evaluated settlements as function of load cycles, loading frequency and axle load are shown in Fig. 14.The settlement is defined as the vertical permanent deformation at the ballast top surface.The results include also loading frequencies of 20 Hz and 40 Hz which is much larger than the estimated relevant frequencies in Table 3 as those might indicate how the settlements are affected by inferior materials.
In order to show the usefulness of the present Finite Element approach, a detailed visualization of the settlements around individual sleepers is shown in Fig. 15 where the settlement along the track is drawn with the sleeper locations indicated in grey.While the loading in the present case is simplified as just one axle at one position is simplified it is clear that the present approach can provide detailed information about settlements at local features like sleepers.

Discussion
The constitutive model, presented in Section 2, can capture the accumulation of permanent deformation in unbound granular materials (UGM) for a wide stress range and loading frequencies as discussed in detail in Section 2.2.In the present study, UGMs material parameters are identified from triaxial tests data by [9] and [26].As discussed in the  previous studies by the authors, [30,31] one potential way to reduce the volume of triaxial testing for UGM characterization is to employ modelling with the Discrete Element Method (DEM) for studying the effects of stress states, gradation and loading rates on the response.
The applicability of the developed model for ballast settlements predictions has been examined for several embankment structures and traffic loads.The maximum levels of confining and deviatoric stresses in ballast layers were in the same range as the ones in the experimental dataset used for model calibration.
As seen from Fig. 14, the loading frequency has a profound influence on the permanent deformation of UGMs and should be accounted for.A simplified analysis has been performed presently based on evaluation of deviatoric stress gradients in longitudinal direction, cf.Fig. 10.In addition to train speed, embankment design was found to affect the loading rates in ballast layers.As seen in Table 3, 5 Hz loading rate may be considered as an upper boundary for the case of slab track while rates up to 12 Hz are relevant for the tracks with sleepers.A more accurate analysis of the local loading rates in ballast would require a full dynamic analysis including wave propagation and material damping, which is beyond the scope of this study.
Once the stress distributions and loading frequencies are determined, equations ( 1)-( 10) may be used to evaluate permanent strains in a ballast layer.As illustrated with Fig. 11, the obtained permanent strain distributions are reasonable, at least in qualitatively.Namely, using concrete slab instead of sleepers significantly reduces the accumulated strains and results in more uniform strain distribution.The load magnitude also has a profound effect on the permanent strain accumulation.For the case embankment with sleepers, increasing the axial load with 70% results in approximately 120% higher maximum permanent strains in the ballast layer.
Settlements in ballast layers may be calculated from permanent strains, by employing the procedure in section 3.2.The obtained settlements are presented in Fig. 14, for all the embankment structures and axle loads examined.With respect to quantitative results, it may first be noticed that the results presented in Fig. 14, are reasonable, with the settlements in the range 5 to 40 mm for the loading frequencies ≤ 20 Hz.This is an encouraging result as no additional calibrations of the model have been performed.The relative effects of structural designs and axial load levels on the settlements accumulations are also reasonable, with e. g. 30 t axial loads resulting in 40-100% higher settlements at 10 6 load cycles compared to 17.5 tones.
With respect to quantitative predictions of ballast settlements, field calibration of the developed computational tool is of crucial importance.In this study, it is assumed that in the beginning of embankments service life, UGMs state in ballast layer and in triaxial test specimen is identical with respect to compaction degree, gradation, etc.This is not necessarily the case, especially regarding compaction.One possible way to correct for this, is to adjust the number of load cycles, attributing some initial load cycles to pre-compaction.
In absence of full-scale measurements, a rough quantitative evaluation of the developed model is performed by comparing the predicted settlements with empirical models from the literature.In Fig. 16, settlements predicted with 4 models [3,[32][33][34] are compared to the developed model predictions.It must be emphasized that the empirical models are intended to capture the total settlements in embankments, in contrast to the ballast layers settlements captured by the present model.Furthermore, the empirical models include the effect of pre-compaction of the ballast, which presently is not accounted for.Nevertheless, the results in Fig. 14 seem to compare favourably to the ones in Fig. 16.Namely, the settlements predicted after 1e6 load cycles, generally differ less than 7 mm from the ones obtained with empirical models, which is in line with the differences between empirical models.Furthermore, the computationally predicted 50-80% effects of changing axle loads (depending on the frequency chosen) is comparable with 60% and 90% obtained from the [32] and [3] models respectively.

Conclusion
This paper addresses the gap found in the literature of having a computational efficient tool for calculating settlements in large embankment structures after long-time traffic loading.In the framework, the stress state in the embankment is calculated by the finite element method and thereafter used for calculating the permanent strain state due to train traffic.The permanent strains are in a last computation step translated into permanent deformations which provide the settlement.
The modelling results clearly show that the frequency largely affects the magnitude of the permanent deformations.It is suggested in the present work that the governing parameter for the development of permanent strain except from the stress magnitude is not only the frequency itself but also the rate of the cyclic stress state.This should be investigated further to account for train speed in a fully quantitative way.The presented results concerning the ballast layer settlements are reasonably close to the settlements predicted by empirical models.Given that no additional calibration on a structural level has been performed, this result is encouraging and indicates that the developed framework captures basic mechanics of permanent deformation accumulation in ballast layers well.
Finally, the pre-compaction of the ballast material is identified as a major factor affecting the predicted settlement levels as most of the deformation occurs during the first 1000 load cycles for the relevant load cases.This effect could explicitly be modelled in the framework if

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

E
.Olsson et al.

Fig. 1 .
Fig. 1.A sketch of the different types of deformation occurring in ballast materials depending on the static and cyclic stress state.From [21].
E.Olsson et al.

Fig. 4 .
Fig. 4. (a) The predicted axial deviatoric strain after 500 000 cycles as function of deviatoric von Mises stress for two different pressure levels, 30 kPa and 60 kPa.(b) Different combinations of deviatoric stresses and hydrostatic pressure to obtain a given deviatoric axial permanent strain after 500 000 load cycles.

Fig. 5 .
Fig. 5. Predicted and experimentally obtained permanent strains for different degree of ballast fouling represented by void contaminant index (VCI).The experimental data is taken from [26].
placements u G are the unknown.To calculate displacements uniquely from strains, boundary conditions are needed.A reduced global strain-displacement matrix is created B red G containing only the columns pertinent to the non-prescribed components of u G .Prescribed nodal displacements are considered by collecting these values in a new column vector u 0 and multiplying this vector by a matrix B 0 G containing the columns of the matrix B G pertinent to those nodal displacements.The remaining unknown nodal displacements are collected in a reduced vector u red G .Finally, an overly determined system of equations for u red G can be formulated as

Fig. 7 .
Fig. 7.The meshed finite element models of the embankments.The high embankment with sleepers is shown in (a) and a magnified view of the load point of this model is shown in (b).The "Low embankment" model with a concrete slab is pictured in (c) with the load point magnified in (d).

Fig. 8 .
Fig. 8. (a) The static pressure field due to gravity in the model "High Embankment" with sleepers.(b) The cyclic von Mises stress field for the same model loaded with 22.5 tonnes axle load.

Fig. 9 .
Fig. 9. (a) Graph of the static pressures along the red dotted line in Fig. 8 for the different investigated models.(b) The cyclic von Mises stress along the red dotted line in Fig. 8 (b) for the different investigated models and axle loads.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.The distribution of the cyclic von Mises stress beneath the rail at the depth where the maximum stress occurs in Fig. 9 (b).Results for the embankment with a concrete slab is shown in (a) and results for the embankment with sleepers in (b) where the position of the sleepers is indicated as grey areas.

Fig. 11 .
Fig. 11.(a) The volumetric strain and (b) the effective deviatoric strain, along the red lines defined in Fig. 8.The results are evaluated at 1 000 000 load cycles at a frequency of 5 Hz.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. Visualization of the permanent vertical deformations after 1 000 000 load cycles of an axle load of 22.5 tonnes at 5 Hz.Results for the "High Embankment" with a concrete slab is shown in (a) and results for the same model with sleepers in (b).

Fig. 13 .
Fig. 13.Distribution of the vertical permanent deformation for the different investigated models with different axle loads and at different number of load cycles.The loading frequency is 5 Hz for all cases.

Fig. 14 .
Fig. 14.Evaluation of the numerically predicted settlement, defined as the vertical permanent deformation at the top of the ballast layer under the rail, as function of load cycles, axle load and loading frequency.

Fig. 15 .
Fig. 15.Predicted settlements along the track or the case for 22.5 tonnes axle load and 10 Hz frequency.The sleeper locations are indicated with grey colour.

Fig. 16 .
Fig. 16.Comparison of the prediction of different empirical models with the predictions at 1 000 000 load cycles from the present modelling work where different axle load levels are shown with different line types for the models accounting for such effects.

Table 2
Geometrical dimensions of the studied ballast geometries.
Fig.6.Sketch of the cross-section of the railway embankment model.The dimensions for the studied geometries are given in Table2.Boundary conditions for the FE-model are indicated.E.Olsson et al.