Applicability of a Three-Layer Model for the Dynamic Analysis of Ballasted Railway Tracks

: In this paper, the three-layer model of ballasted railway track with discrete supports is analyzed to access its applicability. The model is referred as the discrete support model and abbreviated by DSM. For calibration, a 3D ﬁnite element (FE) model is created and validated by experiments. Formulas available in the literature are analyzed and new formulas for identifying parameters of the DSM are derived and validated over the range of typical track properties. These formulas are determined by ﬁtting the results of the DSM to the 3D FE model using metaheuristic optimization. In addition, the range of applicability of the DSM is established. The new formulas are presented as a simple computational engineering tool, allowing one to calculate all the data needed for the DSM by adopting the geometrical and basic mechanical properties of the track. It is demonstrated that the currently available formulas have to be adapted to include inertial effects of the dynamically activated part of the foundation and that the contribution of the shear stiffness, being determined by ballast and foundation properties, is essential. Based on this conclusion, all similar models that neglect the shear resistance of the model and inertial properties of the foundation are unable to reproduce the deﬂection shape of the rail in a general way.


Introduction
Numerical models of the railway track are fundamental tools for the study of their dynamic behavior, with implications for the safety and comfort of rail transport. The importance of these models has increased alongside the speed of the railway vehicles and capacity of the network over the last decades.
The use of 3D FE models is common practice, but reduced models are still relevant, due to simplicity of implementation and results interpretation, and low computational cost. Despite the wide use of the reduced models, there are still relatively few studies about their overall applicability and how to select appropriate parameters based on the properties of the railway track. Most works determine the necessary parameters in a way to match the response of the model to experimental measurements from the track. This means that the range of applicability of the methods developed and the conclusions reached in these works are always, to some extent, limited.
This paper provides a detailed analysis of the viability and applicability of the threelayer discrete support model [1], for simplicity referred to as the discrete support model and abbreviated by DSM, expanding on the work presented in [2] and analyzing the formulas proposed in [3]. The DSM is a 2D model composed of a beam, representing the rails, connected to sets of springs, dampers and masses that represent the support elements of the track: the fastening system, sleepers, ballast, and subgrade. The DSM has significant allowing the reader to calculate all the data needed for the DSM simply by adopting the geometrical and basic mechanical properties of the track constituents.
This paper is organized in the following way: Section 2 describes the DSM and reviews the relevant literature, Section 3 presents the 3D FE model and its validation, Section 4 is dedicated to the calibration of the DSM based on the results of the 3D FE model and derivation of the theoretical expressions and Section 5 summarizes the main conclusions.

DSM
Regarding the evolution of the reduced models, it is worthwhile to mention that discrete supports were introduced in [49] to remove the limitations of the Winkler and Pasternak models, which consist of a beam continuously supported by an elastic or viscoelastic foundation. Newton and Clark [50] were the first to use the abbreviation DSM and to consider the sleepers as discrete mass elements, which separated the stiffness and damping introduced by the fastening system from the one due to the ballast bed and subgrade, creating two layers. In [51], another set of masses that represent the ballast were introduced, and therefore, the stiffness and damping of the ballast and of the subgrade/foundation got also separated, creating three layers. Zhai and Sun [4] then introduced springs and dampers connecting consecutive ballast masses to model the shear behavior of the ballast, which already cover all the DSM components, as depicted in Figure  1. The rail is modelled as a Timoshenko beam, for which the bending ( EI ) and shear stiffness ( * GA ), and the mass per unit length ( m ) are well-known. Likewise, the fastening system stiffness and damping coefficient ( pad K and pad C ) can be determined experimentally. The sleepers experience very little deformation during the normal use of the track, so they are modelled as rigid elements that only contribute with their mass ( sleep M ), which is well-defined. The remaining unknown parameters are the vertical stiffness and damping coefficient of the ballast ( b K and b C ), the vertical stiffness and damping coefficient of the subgrade ( f K and f C ), the shear stiffness and damping coefficient of the ballast ( w K and w C ) and the dynamically activated ballast mass ( b M ). In [3,4], formulas for calculation of the vertical stiffness of the ballast are proposed as an adaptation of the stress cone (pyramid) method, which assumes that the stresses are transmitted from the effective load-bearing area under the sleepers to the ballast, and then dissipate at a constant angle, as shown in Figure 2a). The effective load-bearing area under the sleepers is determined in [45], based on experimental observations and simple geometrical considerations, requiring only the sleeper's dimensions and the track gauge. Then, the vertical stiffness b K can be determined as a function of a distribution angle  b , the Young modulus of the ballast and its depth. In [43], a reduction factor for the vertical stiffness due to superposition of stress cones of 0.5 was proposed based on experimental The rail is modelled as a Timoshenko beam, for which the bending (EI) and shear stiffness (GA * ), and the mass per unit length (m) are well-known. Likewise, the fastening system stiffness and damping coefficient (K pad and C pad ) can be determined experimentally. The sleepers experience very little deformation during the normal use of the track, so they are modelled as rigid elements that only contribute with their mass (M sleep ), which is well-defined. The remaining unknown parameters are the vertical stiffness and damping coefficient of the ballast (K b and C b ), the vertical stiffness and damping coefficient of the subgrade (K f and C f ), the shear stiffness and damping coefficient of the ballast (K w and C w ) and the dynamically activated ballast mass (M b ).
In [3,4], formulas for calculation of the vertical stiffness of the ballast are proposed as an adaptation of the stress cone (pyramid) method, which assumes that the stresses are transmitted from the effective load-bearing area under the sleepers to the ballast, and then dissipate at a constant angle, as shown in Figure 2a. The effective load-bearing area under the sleepers is determined in [45], based on experimental observations and simple geometrical considerations, requiring only the sleeper's dimensions and the track gauge. Then, the vertical stiffness K b can be determined as a function of a distribution angle α b , the Young modulus of the ballast and its depth. In [43], a reduction factor for the vertical stiffness due to superposition of stress cones of 0.5 was proposed based on experimental observations. In [3], superposition of the adjacent stress cones in the longitudinal direction is introduced exactly. The cone volume such defined is then used for the estimation of the dynamically activated ballast mass, M b . The area of the base of the stress cone is further used to calculate the foundation vertical stiffness, K f , but for this, the experimentally determined subgrade reaction modulus must be known. Vibration 2020, 3 FOR PEER REVIEW 2 observations. In [3], superposition of the adjacent stress cones in the longitudinal direction is introduced exactly. The cone volume such defined is then used for the estimation of the dynamically activated ballast mass, b M . The area of the base of the stress cone is further used to calculate the foundation vertical stiffness, f K , but for this, the experimentally determined subgrade reaction modulus must be known. No functional dependence is proposed for all damping coefficients ( b C , f C and w C ), for the shear stiffness of the ballast ( w K ) and, in fact, also for the foundation stiffness f K . There is also no justification for the value of the distribution angle  b . It is therefore the objective of the present work to propose general formulas for all components of the DSM without the necessity to call on experimental data. Alongside this, the range of applicability of the DSM is obtained. Contrary to the arrangement given in Figure 1, it is proven that the dynamically activated part of the foundation must be added to the ballast mass and that the shear stiffness of the ballast must be increased by the contribution from the foundation soils. If this is not respected, then the DSM is unable to reproduce the deflection shape of the rail in a general way.

Three-Dimensional Model
To obtain the reference results necessary to calibrate the DSM, a 3D FE model of the railway track is developed. This model has linear elastic behavior and can be used to obtain static and dynamic vertical displacements in the rail. It is implemented using the AN-SYS parametric design language [51], so the model can be generated in a fully automatic way for each set of parameters.
Static, modal and transient analyses are used to define the necessary discretization and which type of boundary conditions (BCs) are suitable to the problem in study. The static and modal analyses are necessary to validate the stiffness component of the BCs. The 3D FE model is validated by modelling an existing railway track for which the rail displacements due to a train passage are available in [46].
The sleepers, ballast and subgrade are modelled by brick and wedge solid elements, and beam elements are used for the rail and discrete spring-damper elements in the three orthogonal directions simulate the fastening system. The Timoshenko beam has the crosssection properties of a 60E1 rail profile [52] and shear factor of 0.4 as in [5]. The sleepers are spaced by 0.6m and their geometry is based on the DW post-tensioned mono-block concrete sleepers used by the Portuguese railway manager, I.P. S.A. [53,54], but altered to achieve a regular mesh (the following parameters in m are used: length 2.6, bottom width 0.3, top width 0.2 and height 0.22). The moment of inertia is approximately the same as the real section, and the change in volume is corrected by applying a factor to the mass of No functional dependence is proposed for all damping coefficients (C b , C f and C w ), for the shear stiffness of the ballast (K w ) and, in fact, also for the foundation stiffness K f . There is also no justification for the value of the distribution angle α b . It is therefore the objective of the present work to propose general formulas for all components of the DSM without the necessity to call on experimental data. Alongside this, the range of applicability of the DSM is obtained. Contrary to the arrangement given in Figure 1, it is proven that the dynamically activated part of the foundation must be added to the ballast mass and that the shear stiffness of the ballast must be increased by the contribution from the foundation soils. If this is not respected, then the DSM is unable to reproduce the deflection shape of the rail in a general way.

Three-Dimensional Model
To obtain the reference results necessary to calibrate the DSM, a 3D FE model of the railway track is developed. This model has linear elastic behavior and can be used to obtain static and dynamic vertical displacements in the rail. It is implemented using the ANSYS parametric design language [51], so the model can be generated in a fully automatic way for each set of parameters.
Static, modal and transient analyses are used to define the necessary discretization and which type of boundary conditions (BCs) are suitable to the problem in study. The static and modal analyses are necessary to validate the stiffness component of the BCs. The 3D FE model is validated by modelling an existing railway track for which the rail displacements due to a train passage are available in [46].
The sleepers, ballast and subgrade are modelled by brick and wedge solid elements, and beam elements are used for the rail and discrete spring-damper elements in the three orthogonal directions simulate the fastening system. The Timoshenko beam has the crosssection properties of a 60E1 rail profile [52] and shear factor of 0.4 as in [5]. The sleepers are spaced by 0.6 m and their geometry is based on the DW post-tensioned mono-block concrete sleepers used by the Portuguese railway manager, I.P. S.A. [53,54], but altered to achieve a regular mesh (the following parameters in m are used: length 2.6, bottom width 0.3, top width 0.2 and height 0.22). The moment of inertia is approximately the same as the real section, and the change in volume is corrected by applying a factor to the mass of the material of 0.86. Another simplification is that no loss of contact is possible between distinct materials, however, in the absence of sudden changes in stiffness or other irregularities, no loss of contact occurs anyway. Symmetry along the vertical plane parallel to the rails (the xy-plane) is assumed, reducing the size of the model to half, as seen in Figure 3. The shoulder width of 0.5 m and slope 1:2 are assumed for the ballast layer. Two thicknesses will be used of 0.3 and 0.6 m, leading to the bottom width of 4.8 and 6 m. The subgrade is discretized using a regular orthogonal mesh. The regularity of the mesh simplifies the process of implementing BCs and avoids spurious reflections between elements of different sizes [55]. Preliminary tests have also shown that the BCs perform worse at the interface between different element sizes. Convergence studies for static, dynamic and transient analyses led to the choice of the mesh: rail (x) 0.05 m; longitudinal direction (x) 0.075 m, lateral (z) and vertical direction except for sleepers (y) 0.1 m, along sleeper height (y) 0.055 m.
Vibration 2020, 3 FOR PEER REVIEW 2 the material of 0.86. Another simplification is that no loss of contact is possible between distinct materials, however, in the absence of sudden changes in stiffness or other irregularities, no loss of contact occurs anyway. Symmetry along the vertical plane parallel to the rails (the xy-plane) is assumed, reducing the size of the model to half, as seen in Figure  3. The shoulder width of 0.5m and slope 1:2 are assumed for the ballast layer. Two thicknesses will be used of 0.3 and 0.6m, leading to the bottom width of 4.8 and 6m. The subgrade is discretized using a regular orthogonal mesh. The regularity of the mesh simplifies the process of implementing BCs and avoids spurious reflections between elements of different sizes [55]. Preliminary tests have also shown that the BCs perform worse at the interface between different element sizes. Convergence studies for static, dynamic and transient analyses led to the choice of the mesh: rail (x) 0.05m; longitudinal direction (x) 0.075m, lateral (z) and vertical direction except for sleepers (y) 0.1m, along sleeper height (y) 0.055m. The subgrade depth can be defined as the depth at which a significantly rigid substrate is found, or as a reasonable depth after which the deformations due to surface loads are negligible (the so-called active depth of the soil [56,57] [65]. Given this variability, five different depths of the subgrade ( s h ) are studied: 3, 6 and 9 m (relatively shallow subgrade), 25 m (an average depth) and 50 m (a high depth). Not all the depths are considered for all the different analyses, as will be discussed in the relevant section.
All materials are assumed to have linear elastic behavior, since the purpose of the model is to analyze short-term behavior due to vehicle passage. Besides the elastic components, there are discrete linear viscous dampers in the fastening system. No material damping was considered in the ballast and subgrade, to avoid damping the dynamic response. The material properties of the rail and sleepers are well defined, since they are manufactured components that follow the relevant standards [53,54]. The properties of the ballast and subgrade, on the other hand, show great variability across the literature, particularly the Young modulus (see Rodrigues [66]), which generally has the greatest influence on the short-term response of the track. Three values of the Young modulus were chosen for the ballast and subgrade, which cover the relevant range found in the literature while excluding extreme values that are not representative of typical railway The subgrade depth can be defined as the depth at which a significantly rigid substrate is found, or as a reasonable depth after which the deformations due to surface loads are negligible (the so-called active depth of the soil [56,57] [65]. Given this variability, five different depths of the subgrade (h s ) are studied: 3, 6 and 9 m (relatively shallow subgrade), 25 m (an average depth) and 50 m (a high depth). Not all the depths are considered for all the different analyses, as will be discussed in the relevant section.
All materials are assumed to have linear elastic behavior, since the purpose of the model is to analyze short-term behavior due to vehicle passage. Besides the elastic components, there are discrete linear viscous dampers in the fastening system. No material damping was considered in the ballast and subgrade, to avoid damping the dynamic response. The material properties of the rail and sleepers are well defined, since they are manufactured components that follow the relevant standards [53,54]. The properties of the ballast and subgrade, on the other hand, show great variability across the literature, particularly the Young modulus (see Rodrigues [66]), which generally has the greatest influence on the short-term response of the track. Three values of the Young modulus were chosen for the ballast and subgrade, which cover the relevant range found in the literature while excluding extreme values that are not representative of typical railway tracks in good condition. Table 1 summarizes these values, the most typical values of Poisson's ratio and mass density and the properties of the sleepers and rail. The range considered for Young's modulus of the ballast may be considered large, but it was justified by extensive literature search and by summarizing values from more than 40 sources. The histogram for ballast Young's modulus, Poisson's ratio and density are presented in Figure 4.
Vibration 2020, 3 FOR PEER REVIEW 2 tracks in good condition. Table 1 summarizes these values, the most typical values of Poisson's ratio and mass density and the properties of the sleepers and rail. The range considered for Young's modulus of the ballast may be considered large, but it was justified by extensive literature search and by summarizing values from more than 40 sources. The histogram for ballast Young's modulus, Poisson's ratio and density are presented in Figure 4.  For the properties of the fastening system, the Vossloh Zw687a EVA rail pad with direct fastening is used, to match IP S.A. [54]. In agreement with [67], the static stiffness For the properties of the fastening system, the Vossloh Zw687a EVA rail pad with direct fastening is used, to match IP S.A. [54]. In agreement with [67], the static stiffness of 1300 MN/m and 100 MN/m, dynamic stiffness of 3550 MN/m and 280 MN/m and damping coefficient of 36 kNs/m and 10 kNs/m were used, with the first value reflecting the vertical and the second value the lateral and longitudinal directions, respectively (see [66] for more details).
In regard to the elastic part of the BCs, the normal stiffness of 4G/r and tangential stiffness of 2G/r is used at the lateral surfaces, as derived for spherical waves in [68,69], because they were found to work better than the ones derived for the cylindrical waves that were used in [70]. The bottom boundary is considered either fixed or E oed /h and G/h are used for the normal and tangential stiffness, respectively, as in [70]. Here, G stands for the shear modulus, r for the distance from the source, E oed is the oedometric modulus and h is the depth of the foundation not being modelled. At all (non-fixed) boundaries damping coefficients of ρc p for normal and ρc s for tangential direction, as in [68][69][70][71] are used. Here, ρ is the density and c p , c s , are the velocity of propagation of pressure and shear waves, respectively.
The model was tested for different combinations of properties and showed good convergence for the mesh specified for static displacement, the 10 first natural frequencies, the transient response to a pulse load and to a moving load at velocities of 50 and 100 m/s. Finally, the model was validated by comparison with results published in [46], where the measured vertical vibrations of a railway line in Portugal due to the passage of a train at 220 km/h are compared to the results of a 2D FE model. 66 kN per wheel is used in the front bogie. By adapting the properties of the 3D FE model developed for the purpose of this paper, the vertical displacements seen in Figure 5 were obtained. Vibration 2020, 3 FOR PEER REVIEW 2 of 1300MN/m and 100MN/m, dynamic stiffness of 3550MN/m and 280MN/m and damping coefficient of 36kNs/m and 10kNs/m were used, with the first value reflecting the vertical and the second value the lateral and longitudinal directions, respectively (see [66] for more details).
In regard to the elastic part of the BCs, the normal stiffness of 4/ Gr and tangential stiffness of 2/ Gr is used at the lateral surfaces, as derived for spherical waves in [68,69], because they were found to work better than the ones derived for the cylindrical waves that were used in [70]. The bottom boundary is considered either fixed or / oed Eh and / Gh are used for the normal and tangential stiffness, respectively, as in [70]. Here, G stands for the shear modulus, r for the distance from the source, oed E is the oedometric modulus and h is the depth of the foundation not being modelled. At all (non-fixed) boundaries damping coefficients of  p c for normal and  s c for tangential direction, as in [68][69][70][71] are used. Here,  is the density and p c , s c , are the velocity of propagation of pressure and shear waves, respectively. The model was tested for different combinations of properties and showed good convergence for the mesh specified for static displacement, the 10 first natural frequencies, the transient response to a pulse load and to a moving load at velocities of 50 and 100 m/s. Finally, the model was validated by comparison with results published in [46], where the measured vertical vibrations of a railway line in Portugal due to the passage of a train at 220 km/h are compared to the results of a 2D FE model. 66kN per wheel is used in the front bogie. By adapting the properties of the 3D FE model developed for the purpose of this paper, the vertical displacements seen in Figure 5 were obtained.   It can be concluded that the 3D FE model provides a good approximation to the experimental results, particularly when taking into account the variability of the track properties and the random nature of the dynamic response for high frequencies due to the short-wave irregularities of the rails and wheels. These results also confirm that the simplifications introduced in the model do not significantly impact its applicability, namely, the assumption of linear elasticity and the use of moving forces.
The input data for the 3D FE model are taken from [46]. Besides the typical properties related to the rail, rail pads and sleepers, there are other materials characterized by its thickness, Youngs's modulus, density and Poisson's ratio: the ballast layer (0. . The depth of the embankment soils and the stiffness of the natural foundation in relation to that of the former were deemed high enough that the latter can be considered as a rigid substrate, therefore subgrade was implemented with (7 m, 80 MPa, 2040 kg/m 3 , 0.3) because the capping layer, having low thickness, does not influence the equivalent stiffness much. As far as the ballast layer, as the generic model used for optimization has only one layer, (0.6 m, 150 MPa, 1530 kg/m 3 , 0.2) was implemented, where 150 MPa corresponds to an equivalent modulus obtained by considering the two layers characterized by equivalent springs connected in series.

Optimization
The optimization is performed using genetic algorithms (GA), which, being a heuristic method, provides no guarantee that the solution found is the global optimum. Therefore, several independent runs are necessary to guarantee the stability of the results obtained.
In what follows, only stable results will be shown for each set of optimization runs, for the sake of brevity. For all optimization runs, the following parameters are used: the population consists of 100 individuals, the number of generations is 20, the mutation rate is 1% and an elite of two individuals is kept between generations (see [72] for definitions). The objective function Φ is equal to the L 2 -norm of the difference between the rail displacements resulting from the DSM and from the 3D FE model, divided by the L 2 -norm of the displacement obtained by the 3D FE model. Therefore, Φ can theoretically reach zero when both deflection shapes are identical.
The philosophy behind the optimization is the following: firstly, dependence on basic mechanical properties is determined, identifying the proportionality factors; secondly, mechanistic expressions are used to specify these proportionality factors in terms of some geometric or similar specifications of the model. This is very important because the solution of the optimization problem is not unique, i.e., several combinations of parameters can lead to practically identical rail deflection shapes, but for some of those solutions, the parameters have no physical meaning. Expressing the proportionality factors in terms of some geometric parameters impose natural limits on the design space and consequently guarantee the obtainment of physically admissible values.
Besides this, the optimization is divided in two steps. First, the stiffness parameters characterizing the DSM are determined from the static analysis (v = 0 m/s). Then, a dynamic analysis (v = {50, 100} m/s) is used to determine damping and mass parameters. In addition, it is also important to distinguish the individual and combined optimizations. For the individual optimization, a unique set of the discrete values of the parameters is selected and the objective function Φ is used as specified above. For the combined optimization, some parameters are unique (fixed) while others assume all possible combinations. Then the objective function is defined as Φ = max i Φ i , where i identifies the particular combination.
As specified in previous section, the discrete values selected for variation are: • Φ is small for all combinations (1%-4%), so the DSM can adequately approximate the 3D FE model solution for a static load.

•
The optimum value of K f varies linearly with E s and has the greatest influence on Φ.

•
The optimum value of K w varies significantly with both E b and E s (and therefore with the respective shear moduli, G b and G s ), with the latter having a stronger influence, thus the resulting optimum value of K w can be suitably approximated using a linear combination of G b and G s . Although no clear conclusions can be drawn as regarding K b , further analyses will show that the formula proposed in [3] provides a good approximation. Therefore, taking also into account the observations above, the following relations are proposed and will be used in the combined optimization: where E oed s is the oedometric modulus of the subgrade. The proportionality factors ( f K,b , f K,s , f K,w,b , f K,w,s ) are in meters.

Stiffness Parameters (Static Solution), Combined Optimization
The static solution is optimized for all possible combinations of E b and E s simultaneously, while h b is considered unique and h s are grouped in two sets: a typical range h s = {6, 25, 50} m and a shallow range h s = {3, 6, 9} m. Some results can be presented for all three depths from the set, some have to be separated, as shown in Table 2. For the former set, the reference results of the 3D FE model were obtained using the elastic BCs at the bottom, so that only a fraction of the subgrade depth had to be modelled, while for the latter, the full subgrade depth was modelled with a fixed bottom boundary. The results of four proportionality coefficients from Equation (1) are summarized in Table 2. Some qualitative observations are listed below: • f K,b decreases as h b increases, but not as much as for a purely inversely proportional relation, which supports the stress cone theory; • f K,s does not change significantly with h b , which suggests that superposition of the stress cones occurs at a relatively shallow depth (see Figure 2); • f K,s decreases when h s increases, but is stays nearly constant from 25 to 50 m, which suggests that the relation between K f and h s is asymptotic in nature; • f K,w,b increases with h b at a rate higher than direct proportionality; • f K,w,s increases with h s , particularly in the shallower range, but not as much from 25 to 50, again suggesting an asymptotic relation; • f K,w,s also increases with h b , in some cases very significantly, while in others not so much; • f K,b and f K,s,6 are very similar for the two different optimization sets; • Despite differences in f K,w,b and f K,w,s,6 between the two sets, the optimum value of K w is similar.
The differences observed between the two sets for h s = 6 m are because for the second set, the complete subgrade depth was modelled instead of using elastic BCs at the bottom, which introduces some differences in the results. The impact on Φ, however, is not significant.

Damping and Mass Parameters (Dynamic Solution), Individual Optimization
Admitting K b , K w and K f obtained for the static case, the remaining parameters according to Figure 1  • For v = 50 m/s, Φ ranges from 8% to 13%, while for v = 100 m/s, from 8% to 18%; • The worst results for v = 100 m/s occur when E s = 50 MPa. When these combinations are omitted, the maximum value of Φ is 13%, the same that was observed for v = 50 m/s; • The effect of C b and C w on Φ is negligible, and the optimum values are close to zero; • C f has the most effect on Φ, and its optimum value increases with E s , but not with E b ; • M b shows significant variation in its optimum value (2.3-5.0 tons). The average value is higher for h b = 0.6 m than for h b = 0.3 m, in agreement with [3]; • The effect of M b on Φ is noticeable, but relatively small in comparison with C f . Using its average value for all combinations of E b and E s does not significantly affect the solution.
The higher value of Φ for the soft subgrade when v = 100 m/s is because in this case, the velocity is supercritical, as already mentioned. As the load velocity gets closer to the critical velocity, the displacements are significantly amplified, especially in the upward direction. If the load velocity is higher than the critical one, then there is a wider region affected by the displacements and the maximum downward value shifts back gradually from the point of load application. If the two models being compared have different critical velocities, then this limits the validity of the comparison. Obviously, the wave propagation in DSM is somehow limited leading to differences between the two models that cannot be removed by improved optimization. In addition, the omission of horizontal displacements prevents the development of Rayleigh waves in the DSM.
Given the results above, the dampers C b and C w were set to zero for combined optimizations. It should be noted that this would not be the case if material damping was used. Since only radiation/geometric damping is considered, the following relation was adapted from [73]: where f C,rad,s has unit of square meters. As for the mass, M b , although it showed variation with the properties of the foundation, it was preliminarily assumed that the oscillating mass is proportional to the mass density of the ballast, as suggested in [3]:  Table 3, including the resulting value of M b .  Due to the high values of M b obtained by the optimization, it is clear that the dependence on foundation properties is important and it is thus necessary to increase the oscillating mass by contribution from the foundation, as preliminarily indicated in the previous section: where M s is the dynamically activated foundation mass. Then, M will be used in place of M b in the DSM description from Figure 1. This means that the assumptions made in [3] are qualitatively correct, but for the shear spring (K w ) as well as for the dynamically activated mass (M), besides the ballast contribution, foundation contribution must also be added, which is the main new contribution of this paper. In general, it can be concluded that the DSM provides a very good approximation to the 3D FE model, when the velocity of the load is not too close to the critical. Nevertheless, the approximation for velocities close to the critical can still be considered as reasonable, as can be seen in Figure 6.

Ballast Vertical Stiffness
The stress cone model assumes that stresses are transmitted from the effective loadbearing area under the sleeper and then dissipate at a constant angle, α b . This assumption approximates the static stress distribution of the 3D FE model, as shown in Figure 7. w activated mass ( M ), besides the ballast contribution, foundation contribution must also be added, which is the main new contribution of this paper.
In general, it can be concluded that the DSM provides a very good approximation to the 3D FE model, when the velocity of the load is not too close to the critical. Nevertheless, the approximation for velocities close to the critical can still be considered as reasonable, as can be seen in Figure 6.

Ballast Vertical Stiffness
The stress cone model assumes that stresses are transmitted from the effective load-  In [3], the superposition of the stress cones in the longitudinal direction is considered. Here, the superposition of the cones in both directions (i.e., also in the transversal direction between rails, as shown in Figure 8) is taken into account. Figure 9 also illustrates the

Ballast Vertical Stiffness
The stress cone model assumes that stresses are transmitted from the effective load-  In [3], the superposition of the stress cones in the longitudinal direction is considered. Here, the superposition of the cones in both directions (i.e., also in the transversal direction between rails, as shown in Figure 8) is taken into account. Figure 9 also illustrates the In [3], the superposition of the stress cones in the longitudinal direction is considered. Here, the superposition of the cones in both directions (i.e., also in the transversal direction between rails, as shown in Figure 8) is taken into account. Figure 9 also illustrates the effective load-bearing area of the sleeper. According to [45], the effective length of the sleeper loaded under each rail is given by: where l sleep is the length of the sleeper and l g is the track gauge.
Vibration 2020, 3 FOR PEER REVIEW 2 effective load-bearing area of the sleeper. According to [45], the effective length of the sleeper loaded under each rail is given by:    The vertical stiffness of the ballast is the inverse of the integral of the virtual strain over its depth due to a vertical load at the top. Due to superposition, it is necessary to combine three different sections: Figure 9. Geometry of the stress cone with superposition in both directions.
The total depth of the ballast is h b , while h x and h z are the height at which adjacent stress cones overlap in the longitudinal and transversal directions, respectively: where l s is the spacing between sleepers. Thus, the dimensions of the cone-base are: The vertical stiffness of the ballast is the inverse of the integral of the virtual strain over its depth due to a vertical load at the top. Due to superposition, it is necessary to combine three different sections:

Ballast Shear Stiffness
The shear stiffness of the ballast, here designated as K w,b does not come directly from the stress cone assumption, and [3] does not provide a way to estimate it. It is proposed here to use classical expressions from beam theory for shear and bending stiffness, adapted to a representative ballast area. Given that the vertical stiffness of the ballast is calculated based on the stress cone distribution, it is reasonable to assume that the representative ballast area, A b , is the cross-sectional area of the stress cone, as depicted in Figure 10: here to use classical expressions from beam theory for shear and bending stiffness, adapted to a representative ballast area. Given that the vertical stiffness of the ballast is calculated based on the stress cone distribution, it is reasonable to assume that the representative ballast area, b A , is the cross-sectional area of the stress cone, as depicted in Figure 10: Figure 10. Cross-section of the stress cone in the transversal direction.
The shear stiffness , wb K is then equal to:

Subgrade Stiffness
In [3], it is assumed that the reaction modulus of the subgrade is known from experimental measurements, and therefore the vertical stiffness of the subgrade, f K , can be calculated by multiplying this value by the area of the base of stress cone of the ballast. In this paper, a general expression will be proposed based on the Vlasov model [74] that prescribes evolution of the vertical displacement along the active depth of the soil, s h : z h Figure 10. Cross-section of the stress cone in the transversal direction.
The shear stiffness K w,b is then equal to:

Subgrade Stiffness
In [3], it is assumed that the reaction modulus of the subgrade is known from experimental measurements, and therefore the vertical stiffness of the subgrade, K f , can be calculated by multiplying this value by the area of the base of stress cone of the ballast. In this paper, a general expression will be proposed based on the Vlasov model [74] that prescribes evolution of the vertical displacement along the active depth of the soil, h s : where γ is a parameter in m −1 that describes the rate at which the vertical displacement in the foundation decreases with depth. When γ = 0, the displacement decreases linearly with depth and the vertical stress in the soil is constant. Both the subgrade vertical reaction modulus, K s , and the shear reaction modulus, K s,P , are calculated by minimizing the total strain energy due to a uniform surface load ( [74]): These values, when multiplied by the area of the base of stress cone of the ballast (A f = l x l z ), give the subgrade stiffness, K f , and the discrete rotational stiffness due to shear, which can then be used to calculate the shear stiffness of the subgrade, K w,s : Figure 11 shows the variation of K s and K s,P with depth for different values of γ. Except for γ = 0, both values become nearly constant after a certain depth of the model, which is lower when γ is higher. This is in agreement with the results of the 3D FE model, for which the static displacement for h s = 25 m and h s = 50 m was very similar, as were the numerical optimal values of the DSM (Table 2). Applying the principles above to a substructure with two distinct layers shows that, as long as the bottom layer is significantly thicker than the top layer, the total vertical reaction modulus can be approximated as a series of springs, while the shear reaction modulus can be approximated as parallel rotation springs. This is consistent with the fact that the shear stiffness is proportional with the layer's depth, and agrees with the numerical finding that the shear stiffness element of the DSM, K w , can be approximated as a linear combination of G b and G b .  (Table 2). Applying the principles

Dynamically Activated Mass
In [3], the ballast mass is calculated from the volume of the stress cone, which is here adapted to consider superposition in both directions: The other contribution specified in Equation (3) is based on [75] where vibrating mass of the soil under foundation is analyzed. Using these formulas and taking into account that the only half the track is being modelled, the following definition of the subgrade mass is proposed:

Subgrade Damping
For the foundation damping, the formula proposed in [73] specifies: Vibration 2021, 4 166 where A f is the area of the foundation, in this case, the area of the base of the stress cone, and c z is the rate of absorption, when c z = 1, in theory, full absorption of incident waves occurs.

Fitting the Decisive Parameters to the Numerical Results
The expressions proposed in Equations (1), (3), (7), (8), (10), (12)-(17) define all the parameters of the DSM from the known properties of the track, with the exception of the angle of stress distribution in the ballast, α b , the rate of displacement decrease with depth in the subgrade, γ, and the rate of absorption due to radiation damping, c z .
To determine the values for these parameters, the expressions above are fitted to the numerical results obtained in the previous sections. Since there are only a few discrete points to fit, a simple norm of the relative difference between the theoretical and numerical results is used.
The first parameter to optimize is α b , which must fit both ballast depths, h b = 0.3 m and h b = 0.6 m. The resulting optimum value for the first set in Table 2 is α b = 49.8 • , and produces an error of ±9%, while for the second set in Table 2 it is α b = 47.6 • with an error of ±4%. The values are shown in Table 4. Using these values of α b , it is possible to define the area of the base of the stress cone and proceed to the optimization of γ and f K,s . The optimum value for the first set in Table 2 is γ = 0.268 m −1 , with an error of f K,s from −16% to +18%, and for the second set in Table 2 it is γ = 0.324 m −1 , with an error of f K,s from −11% to +15%.
Then, it is possible to use the above determined values of α b and γ to compute the shear stiffness of the ballast and the subgrade. Table 5 shows a summary of these results. It is seen that the error is quite high, nevertheless, this is not very important, because the main shear contribution comes from the foundation. For the subgrade shear stiffness, using the optimized values of α b and γ, the error of f K,w,s ranges from −9% to +34% for the first set in Table 2 and from −29% to −14% for the second set. Both f K,s and f K,w,s are plotted as a function of h s in Figure 12, for a rough average of the optimum values of γ = 0.3 m −1 .
It can be concluded that the approximation is quite reasonable, especially taking into account that the optimum values of γ were obtained only from optimizing of f K,s .
The value of Φ obtained for the static solution using these expressions ranges from 3.3% to 14.8% for the first optimization set in Table 2 and from 3.3% to 13.8% for the second set. Vertical displacements for the representative load of 40 kN are compared in Figure 13. It can be concluded that the agreement between DSM with the proposed expressions and the 3D FE model is excellent.
It is seen that the error is quite high, nevertheless, this is not very important, because the main shear contribution comes from the foundation. For the subgrade shear stiffness, using the optimized values of  b and  , the error of ,, K w s f ranges from -9% to +34% for the first set in Table 2 and from -29% to -14% for the second set. Both The value of  obtained for the static solution using these expressions ranges from 3.3% to 14.8% for the first optimization set in Table 2 and from 3.3% to 13.8% for the second   Table 6. It is clear that the proposed equations provide a good approximation for the lower velocity, but not for the high velocity, for which the optimal mass is 20% to 30% lower, likely because the latter is close to the critical velocity of the track for a soft subgrade. Since the proposed equations do not account for this effect, they work well for low to moderate velocities, but not for high velocities, for which the dynamic amplification is more pronounced.
As for the radiation damping, the numerical optimal values in Table 3 are used to deduce the optimum value of z c . The results are summarized in Table 7. It can be concluded that 0.4 z c = is a good recommendation for railway models with properties inside the range studied here. Table 7. Optimum values of z c , based on Table 3. Further, the same comparison must be carried out for the dynamic case. Using the optimum value of α b the error in the proposed expression can be calculated. The results are summarized in Table 6. It is clear that the proposed equations provide a good approximation for the lower velocity, but not for the high velocity, for which the optimal mass is 20% to 30% lower, likely because the latter is close to the critical velocity of the track for a soft subgrade. Since the proposed equations do not account for this effect, they work well for low to moderate velocities, but not for high velocities, for which the dynamic amplification is more pronounced.
As for the radiation damping, the numerical optimal values in Table 3 are used to deduce the optimum value of c z . The results are summarized in Table 7. It can be concluded that c z = 0.4 is a good recommendation for railway models with properties inside the range studied here. The value of Φ obtained for the steady-state dynamic solution using these expressions ranges from 8.1% to 13.1% for v = 50 m/s and from 8.3% to 43.6% for v = 100 m/s. The higher errors are always observed for the soft subgrade (E s = 50 MPa), particularly for v = 100 m/s, as shown in Figure 14, due to the aforementioned dynamic amplification factor due to the critical velocity. It is of note that the maximum upward and downward displacements for the worst case are relatively close to the ones obtained on the 3D FE model, so the DSM and the proposed expressions provide a reasonable approximation to the rail displacement even when the load is moving with the velocity is close to the velocity of propagation of elastic waves in the soil.
Finally, the results of the DSM are compared with the experimental results published in [46], as was the case for the 3D FE model, using The vertical displacement obtained in the DSM due to the passage of the full train and a single bogie at 220 km/h is presented in Figure 15. It can be concluded that the results of the DSM are very close to the ones of the 3D FE model and provide a good approximation of the experimental results. In fact, the results are in better agreement than the ones of the 2D model presented in [46], where the It is of note that the maximum upward and downward displacements for the worst case are relatively close to the ones obtained on the 3D FE model, so the DSM and the proposed expressions provide a reasonable approximation to the rail displacement even when the load is moving with the velocity is close to the velocity of propagation of elastic waves in the soil.
Finally, the results of the DSM are compared with the experimental results published in [46], as was the case for the 3D FE model, using α b = 50 • , γ = 0.3 m −1 and c z = 0.4. The vertical displacement obtained in the DSM due to the passage of the full train and a single bogie at 220 km/h is presented in Figure 15. It is of note that the maximum upward and downward displacements for the worst case are relatively close to the ones obtained on the 3D FE model, so the DSM and the proposed expressions provide a reasonable approximation to the rail displacement even when the load is moving with the velocity is close to the velocity of propagation of elastic waves in the soil.
Finally, the results of the DSM are compared with the experimental results published in [46], as was the case for the 3D FE model, using The vertical displacement obtained in the DSM due to the passage of the full train and a single bogie at 220 km/h is presented in Figure 15. It can be concluded that the results of the DSM are very close to the ones of the 3D FE model and provide a good approximation of the experimental results. In fact, the results are in better agreement than the ones of the 2D model presented in [46], where the multibody vehicle model was considered, justifying the assumption of moving forces It can be concluded that the results of the DSM are very close to the ones of the 3D FE model and provide a good approximation of the experimental results. In fact, the results are in better agreement than the ones of the 2D model presented in [46], where the multibody vehicle model was considered, justifying the assumption of moving forces adopted here.

Comparison with Other Works
As most of the literature does not present formulas for establishing the model data, as a reference for comparison with other works the model presented in [3] is chosen. In [3], values of K b and M b are calculated using α b = 35 • and subgrade stiffness is calculated using the given subgrade modulus K s . Therefore, the comparison is performed in the following way: the data kept in both models are related to the rail: EI = 6. There is no indication for E s in [3]. Thus, K s = 90 MPa/m given in [3] is used with the optimized γ = 0.3 m −1 and h s = 6 m, to give E s = 295.7 MPa. Table 8 summarizes the differences in both models. It is necessary to highlight that values for K w and M, given in Table 8 for the model developed in this paper, are very high due to the contribution from the foundation, as justified by optimizations. In fact, it would be more adequate to include them in a separate foundation layer.
The representative force of 40 kN is used for results comparison, as in previous sections. Static displacement is depicted in Figure 16.  Table 8 summarizes t differences in both models. It is necessary to highlight that values for w K and M , given in Table 8 for the mo developed in this paper, are very high due to the contribution from the foundation, justified by optimizations. In fact, it would be more adequate to include them in a separ foundation layer. The representative force of 40kN is used for results comparison, as in previous s tions. Static displacement is depicted in Figure 16. It is seen that the high shear stiffness of the model attenuates the upward displa ments which are not present in the 3D FE model. The deflection shape agreement is bet for the model presented in this paper. Further, the stabilized steady-state deflection sha Figure 16. Static vertical displacement of the rail. Comparison of the model proposed in this paper, the model from [3] and the 3D FE model It is seen that the high shear stiffness of the model attenuates the upward displacements which are not present in the 3D FE model. The deflection shape agreement is better for the model presented in this paper. Further, the stabilized steady-state deflection shape is compared to the moving load in Figure 17. . Comparison of the model proposed in this paper, the model from [3] and the 3D FE model. Finally, the maximum displacement values in the absolute value are extracted a function of velocity. In Figure 18, the maximum value over the structure is plotted gether with the value in the load position. The values are normalized by static displa ment. The typical fact can be observed that by getting closer to the critical velocity, maximum downward displacement does not occur at the load position. It is confirm that the increase in oscillating mass as suggested in this paper makes the prediction of critical velocity more reasonable, bearing in mind the value of s E . This concludes comparison. It was shown that the new formulas led to much better agreement with the 3D model, not only in regard to the deflection shape, but also for the prediction of the criti velocity.

Conclusions
The present work establishes the validity and the range of applicability of the DS It proposes general formulas for identifying all DSM properties without the necessity Finally, the maximum displacement values in the absolute value are extracted as a function of velocity. In Figure 18, the maximum value over the structure is plotted together with the value in the load position. The values are normalized by static displacement. The typical fact can be observed that by getting closer to the critical velocity, the maximum downward displacement does not occur at the load position. It is confirmed that the increase in oscillating mass as suggested in this paper makes the prediction of the critical velocity more reasonable, bearing in mind the value of E s . This concludes the comparison. . Comparison of the model proposed in this paper, the model from [3] and the 3D FE model. Finally, the maximum displacement values in the absolute value are extracted a function of velocity. In Figure 18, the maximum value over the structure is plotted gether with the value in the load position. The values are normalized by static displa ment. The typical fact can be observed that by getting closer to the critical velocity, maximum downward displacement does not occur at the load position. It is confirm that the increase in oscillating mass as suggested in this paper makes the prediction of critical velocity more reasonable, bearing in mind the value of s E . This concludes comparison. It was shown that the new formulas led to much better agreement with the 3D model, not only in regard to the deflection shape, but also for the prediction of the criti velocity.

Conclusions
The present work establishes the validity and the range of applicability of the DS It proposes general formulas for identifying all DSM properties without the necessity calling on experimental results. The formulas proposed are simple and straight-forw to use and avoid the necessity of numerical calibration.
It was shown that the DSM shows a good approximation to the 3D FE model for whole range of parameters considered-ballast and subgrade depth and the Young m It was shown that the new formulas led to much better agreement with the 3D FE model, not only in regard to the deflection shape, but also for the prediction of the critical velocity.

Conclusions
The present work establishes the validity and the range of applicability of the DSM. It proposes general formulas for identifying all DSM properties without the necessity of calling on experimental results. The formulas proposed are simple and straight-forward to use and avoid the necessity of numerical calibration.
It was shown that the DSM shows a good approximation to the 3D FE model for the whole range of parameters considered-ballast and subgrade depth and the Young moduli. However, since the DSM does not model the elastic wave propagation in the soil, it is less reliable when the load moves at a speed close to the critical velocity. Good agreement was observed up to a speed of 75% of the Rayleigh wave velocity in the subgrade (8%-13% overall error), while at a speed very close to the Rayleigh waves, the error became more significant (18%).
All DSM components can be determined from basic mechanical and geometrical properties, except for three fundamental values, for which the following values were found to provide a good approximation: Regarding the formulas available in the literature, it was concluded that the oscillating mass of the model cannot represent only the dynamically activated ballast mass, but the contribution from the foundation must be added. The same conclusion was taken regarding the shear stiffness of the model: the value determined for the ballast has to be superposed with the contribution coming from the foundation. The shear stiffness of the model is essential for achieving a good agreement. This means that reduced models that are focusing mainly on the vertical dynamic equilibrium and/or are neglecting the dynamically activated mass in the foundation, cannot represent the rail deflection shape accurately within the range of typical track properties even in the subcritical velocity range.