Next Article in Journal
Comparison of Reduction Methods for Finite Element Geometrically Nonlinear Beam Structures
Next Article in Special Issue
Design and Modeling of Viscoelastic Layers for Locomotive Wheel Damping
Previous Article in Journal
A Generalized Index for the Assessment of Helicopter Pilot Vibration Exposure
Previous Article in Special Issue
Simulation of Torsional Vibration of Driven Railway Wheelsets Respecting the Drive Control Response on the Vibration Excitation in the Wheel-Rail Contact Point
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Track Systems and Development, Infrastructure, Banedanmark, DK-2450 Copenhagen, Denmark
2
Departamento de Engenharia Civil, Faculdade de Ciências e Tecnologia, Universidade Nova de Lisboa, 2829-516 Caparica, Portugal
3
IDMEC, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal
*
Author to whom correspondence should be addressed.
Submission received: 1 January 2021 / Revised: 22 January 2021 / Accepted: 17 February 2021 / Published: 22 February 2021
(This article belongs to the Special Issue Railway Dynamics and Maintenance)

Abstract

:
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 finite 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 fitting 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 deflection shape of the rail in a general way.

1. 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 three-layer 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 advantages over the more detailed 3D FE model, which presents a high computational cost, large number of results to analyze, the need for special boundary conditions, necessity for high level of discretization, etc. Given its computational efficiency, the DSM and other, even simpler models, are often used to study complex phenomena that arise in different railway transportation problems. Examples of published works dealing with specific railway problems include: (i) modelling the vehicle–track dynamic interaction on plain line [4,5,6,7,8,9,10,11,12,13] but also in the case when the line is passing over bridge [14] or crossing a transition zone [15,16]; (ii) assessing and diagnosing the deterioration of track components [17,18,19,20,21,22], including switches and crossings [23,24,25], rail wear/corrugation [26,27,28,29] and hanging sleepers [30,31]; (iii) studying other problems related to wheel–rail interaction, including the effect of wheel flats [32,33,34,35,36,37] and rolling noise [38,39,40,41,42].
Published works with the aim of defining the properties of reduced models with discrete supports consider the ballast pyramid model referred to in [43,44] and the Saller assumption from 1932 which can be consulted in [45]. The pyramid model is further detailed in [4], under the name stress cone model, where it is used not only for the definition of the vertical ballast stiffness, but also for the identification of the oscillating ballast mass. Further improvements include the superposition of the cones between adjacent sleepers, i.e., in the longitudinal direction [3]. Nevertheless, these works have two notable limitations: (i) they do not propose theoretical expressions for all the elements of the model, (ii) they each validate their models and theoretical expressions by comparison with experimental measurements of a single railway track, and therefore cannot prove that they are widely applicable to a range of different track properties. The two limitations are particularly serious when considered together; since some parameters of the models are obtained by fitting its results to experimental measurements for a single set of track parameters, it is not possible to unequivocally determine if the other parameters are intrinsic to the physical phenomena, or if the good approximation is obtained simply by virtue of the fitting process.
Therefore, the present work aims to provide a detailed analysis of the viability and applicability of the DSM, to derive formulas for its components that are still missing and to enhance formulas that are already available. This is achieved by fitting the DSM steady-state vertical displacement of the rail under moving force to reference solutions that cover a wide range of representative values of the track. The reference solutions are obtained by a linear elastic 3D FE model, which is validated by the experimental measurements published by Paixão [46]. The use of a linear elastic model is appropriate because the focus is on the track’s short-term dynamic behavior. Analyzing the deflection shape of the rail is also adequate, because obtaining a good approximation of the steady-state rail displacement is a necessary prerequisite for studying more complex phenomena.
The DSM is fitted to the reference solutions using genetic algorithms, first individually for each combination of track parameters, to evaluate its ability to approximate the results of the 3D FE model and to define proportionality relationships between the properties of the track and those of the DSM. After that, optimization is conducted simultaneously for various combinations of track parameters. Based on these results, theoretical expressions for all parameters of the DSM are derived as functions of the geometry and mechanical properties of the track. It is shown that the currently available formulas have to be adapted in a way to include inertial effects of the dynamically activated part of the foundation and that the contribution of the shear stiffness which must account not only for ballast but also for foundation, is essential. This is in agreement with [47,48], where it is shown that if the shear stiffness of the model is neglected, then the critical velocity of a uniformly moving force tends to zero with the dynamically activated mass of the model, which is not physical. With this conclusion, all similar models that neglect the shear resistance of the model and inertial properties of the foundation are unable to reproduce the deflection shape of the rail in a general way. In addition to these new conclusions, the stress cone formulas are improved by incorporating the superposition of the cones in the lateral direction. All formulas are presented as a simple computational engineering tool, 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.

2. 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 ( E I ) and shear stiffness ( G A * ), and the mass per unit length ( m ) are well-known. Likewise, the fastening system stiffness and damping coefficient ( K p a d and C p a d ) 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 s l e e p ), 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.
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.

3. 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 cross-section 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.
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]). A great variety of such values can be found in the literature: some adopt values of the order of 3 m [58,59], while others consider much higher values, from 10 m [60,61] to 50 m [62,63,64] and as high as 70 m [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.
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 4 G / r and tangential stiffness of 2 G / 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 o e d / 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 o e d 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.
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.3 m, 130 MPa, 1530 kg/m3, 0.2), sub-ballast layer (0.3 m, 200 MPa, 1935 kg/m3, 0.3), capping layer (0.2 m, 2820 MPa, 1935 kg/m3, 0.3), embankment soils (~7 m, 80 MPa, 2040 kg/m3, 0.3) and natural foundation (~3.5 m, 300 MPa, 2040 kg/m3, 0.3). 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/m3, 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/m3, 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.

4. 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 L2-norm of the difference between the rail displacements resulting from the DSM and from the 3D FE model, divided by the L2-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:
  • The depth of the ballast ( h b = { 0.3 , 0.6 }   m );
  • The Young modulus of the ballast ( E b = { 50 , 150 , 300 }   MPa );
  • The depth of the subgrade ( h s = { 3 , 6 , 9 , 25 , 50 }   m );
  • The Young modulus of the subgrade ( E s = { 50 , 100 , 150 }   MPa );
  • The load velocity ( v = { 0 , 50 , 100 }   m / s ).
Since the analysis is linear and velocities are mostly subcritical, the smooth evolution of the response on all parameters can be assumed allowing for interpolation. For the fixed values of the Poisson ratio and density, the lowest velocity of propagation of Rayleigh waves that is considered is approximately 92.3 m/s in the subgrade, which means that the load velocity of v = 100   m / s is supercritical in this case.

4.1. Identification of Dependencies on Mechanical Properties

4.1.1. Stiffness Parameters (Static Solution), Individual Optimization

The parameters to be determined (compare with Figure 1) are K b , K w and K f . First, the individual optimizations are performed, to establish some general trends. Each of the 18 possible combinations of E b , E s and h b are assumed separately, while h s is kept at a constant value of 6 m. The following observations were made:
  • Φ 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 .
  • Above a certain limit, the value of K b does not have a significant impact on Φ , and so the optimum values are spurious, with no clear dependence on E b , E s or h b .
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:
K b = f K , b E b K f = f K , s E s oed K w = f K , w , b G b + f K , w , s G s
where E s o e d 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.

4.1.2. 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.

4.1.3. 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 are: C b , C w , C f and M b . They are optimized for the dynamic case with a load velocity v = { 50 , 100 }   m / s . For simplicity, only a subgrade depth of h s = 6   m is assumed, using the viscoelastic BCs discussed in the previous section. The following observations were made for the individual optimizations:
  • 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]:
C f = f C , r a d , s 3.4 G s ρ s π ( 1 ν s )
where f C , r a d , 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]: M b = f M , b ρ b where f M , b has a unit of cubic meters.

4.1.4. Damping and Mass Parameters (Dynamic Solution), Combined Optimization

The dynamic solutions for v = 50   m / s and v = 100   m / s are optimized for all possible combinations of E b and E s simultaneously, but h b is kept fixed and only h s = 6   m is assumed. The results are summarized in Table 3, including the resulting value of M b .
Some qualitative observations:
  • f C , r a d , s increases with h b in a similar way to f K , s . This suggests that the area of the base of the stress cone influences the rate of radiation damping, which agrees with [73];
  • f M increases by 0.645 and 0.323 m3 with the increase in h b for v = 50   m / s and v = 100   m / s , respectively;
  • f C , r a d , s for v = 100   m / s is slightly lower than for v = 50   m / s for h b = 0.3   m , but slightly higher for h b = 0.6   m . These differences are likely due to the stochasticity of the optimization method;
  • f M for v = 100   m / s is 25%–31% lower than for v = 50   m / s .
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:
M = M b + M s
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.

4.2. Identification of Dependencies on Geometrical and Other Properties

4.2.1. Ballast Vertical Stiffness

The stress cone model assumes that stresses are transmitted from the effective load-bearing 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.
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:
l e = l s l e e p l g
where l s l e e p is the length of the sleeper and l g is the track gauge.
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:
h x = min ( ( l s l b ) / ( 2 tan α b ) , h b ) ,   h ¯ x = h b h x h z = min ( ( l g l e ) / ( 2 tan α b ) , h b ) ,   h ¯ z = h b h z
where l s is the spacing between sleepers. Thus, the dimensions of the cone-base are:
l x = min ( l b + 2 h x tan α b , l s ) l z = l e + tan α b ( h b + h z )
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:
f K , b = 1 / ( f b , 1 + f b , 2 + f b , 3 ) ,   K b = E b / ( f b , 1 + f b , 2 + f b , 3 )
f b , 1 = ln { l b [ l e + 2 tan α b min ( h x , h z ) ] l e [ l b + 2 tan α b min ( h x , h z ) ] } / [ 2 ( l b l e ) tan α b ]
f b , 2 = { ln [ ( 1 2 l b + h x tan α b ) { l e + tan α b [ h z + 1 2 ( l g l e ) ] } ( 1 2 l b + h z tan α b ) { l e + tan α b [ h x + 1 2 ( l g l e ) ] } ] tan α b [ ( 2 l e l b + ( l g l e ) tan α b ) ] , h z < h x ln ( l e + 2 h z tan α b l e + 2 h x tan α b ) / ( 2 l s tan α b ) , h z h x
f b , 3 = ln { l e + tan α b [ h b + 1 2 ( l g l e ) ] l e + tan α b [ max ( h x , h z ) + 1 2 ( l g l e ) ] } / ( l s tan α b )

4.2.2. 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:
A b = ( h b 2 + 2 h b h z h z 2 ) tan α b / 2 + l e h b
The shear stiffness K w , b is then equal to:
K w , b = ( ( G b A b * l s ) 1 + ( 12 E b I b l s 3 ) 1 ) 1 ,   A b * = 5 6 A b ,   I b = h b 3 ( l e 2 + 4 l e l z + l z 2 ) 36 ( l e + l z )

4.2.3. 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 :
f ( y ) = sinh ( γ ( h s y ) ) sinh ( γ 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]):
K ¯ s = 0 h s E s o e d d 2 f d y 2 d y = E s o e d γ sinh ( γ h s ) cosh ( γ h s ) + γ h s 2 sinh 2 ( γ h s ) K ¯ s , P = 0 h s G s f 2 ( y ) d y = G s sinh ( γ h s ) cosh ( γ h s ) γ h s 2 γ sinh 2 ( γ h s )
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 :
K f = K ¯ s A f K w , s = K ¯ s , P A f / l s 2
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 .

4.2.4. 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:
f M , b = f M , b , 1 + f M , b , 2 + f M , b , 3 ,   M b = ( f M , b , 1 + f M , b , 2 + f M , b , 3 ) ρ b
f M , b , 1 = 4 3 tan 2 α b min ( h x , h z ) 3 + ( l b + l e ) tan α b min ( h x , h z ) 2 + l b l e min ( h x , h z )
f M , b , 2 = { ( h x h z ) ( l b l e + 2 3 tan 2 α b ( h x 2 + h x h z + h z 2 + 3 ( l g l e ) ( h x + h z ) ) + 1 2 tan α b ( l b ( l g l e + h x + h z ) + 2 l e ( h x + h z ) ) ) , h z < h x l s ( h z h x ) ( l e + tan α b ( h x + h z ) ) , h z h x
f M , b , 3 = 1 2 l s ( h b max ( h x , h z ) ) ( 2 l e + tan α b ( l g l e + h b + max ( h x , h z ) ) )
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:
M s = 2 ρ s 1 ν s ( 2 l x l z π ) 3 / 2

4.2.5. Subgrade Damping

For the foundation damping, the formula proposed in [73] specifies:
f C , r a d , s = c Z A f
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.

5. 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.
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 α 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 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.

6. 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: E I = 6.62   MNm 2 , m = 60.64   kg / m ; to the sleepers: M s l e e p = 125.5   kg , l s = 0.545   m , l b = 0.273   m , l e = 0.95   m ; to the rail pads: K p a d = 65   MN / m , C p a d = 75   kNs / m and to the ballast: E b = 110   MPa , h b = 0.45   m , ρ b = 1800   kg / m 3 . Ballast vertical and shear damping are neglected as in the final optimization steps.
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.
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.
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.
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.

7. 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:
  • Angle of stress distribution in the ballast: α b = 50 °;
  • Rate of displacement reduction with depth in the subgrade: γ = 0.3   m - 1 ;
  • Rate of absorption for the radiation damping of the subgrade: c z = 0.4 .
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.

Author Contributions

This research article is based on Ph.D. thesis of the first author, A.F.S.R., supervised by the second author, Z.D. Conceptualization, Z.D. and A.F.S.R.; methodology, A.F.S.R. and Z.D.; software, A.F.S.R.; validation, A.F.S.R.; writing—original draft preparation, A.F.S.R.; writing—review and editing, Z.D.; supervision, Z.D. Both authors have read and agreed to the published version of the manuscript.

Funding

This research was funded in the form of a PhD grant, SFRH/BD/75375/2010 and supported by project UIDB/50022/2020.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

This work was supported by the Portuguese Foundation for Science and Technology (FCT), through IDMEC, under LAETA, project UIDB/50022/2020.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Knothe, K.L.; Grassie, S.L. Modelling of railway track and vehicle-track interaction at high frequencies. Veh. Syst. Dyn. 1993, 22, 209–262. [Google Scholar] [CrossRef]
  2. Rodrigues, F.S.; Dimitrovová, Z. Applicability of Simplified Models of Railways Tracks Obtained by Optimization and Fitting Techniques. In Proceedings of the 6th International Conference on Engineering Optimization (EngOpt2018), Lisbon, Portugal, 17–19 September 2018. [Google Scholar] [CrossRef]
  3. Zhai, M.W.; Wang, K.Y.; Lin, J.H. Modelling and experiment of railway ballast vibrations. J. Sound Vib. 2004, 270, 673–683. [Google Scholar] [CrossRef]
  4. Zhai, W.M.; Sun, X. A detailed model for investigating vertical interaction between railway vehicle and track. Veh. Syst. Dyn. 1994, 23, 603–615. [Google Scholar] [CrossRef]
  5. Dahlberg, T. Vertical dynamic train/track interaction—verifying a theoretical model by full-scale experiments. Veh. Syst. Dyn. 1995, 24, 45–57. [Google Scholar] [CrossRef]
  6. Zhai, W.; Cai, Z. Dynamic interaction between a lumped mass vehicle and a discretely supported continuous rail track. Comput. Struct. 1997, 63, 987–997. [Google Scholar] [CrossRef]
  7. Oscarsson, J.; Dahlberg, T. Dynamic train/track/ballast interaction—Computer models and full-scale experiments. Veh. Syst. Dyn. 1998, 29, 73–84. [Google Scholar] [CrossRef]
  8. Oscarsson, J. Dynamic train-track interaction: Variability attributable to scatter in the track properties. Veh. Syst. Dyn. 2002, 37, 59–79. [Google Scholar] [CrossRef]
  9. Oscarsson, J. Simulation of train-track interaction with stochastic track properties. Veh. Syst. Dyn. 2002, 37, 449–469. [Google Scholar] [CrossRef]
  10. Bureika, G.; Subačius, R. Mathematical model of dynamic interaction between wheel-set and rail track. Transport 2002, 17, 46–51. [Google Scholar] [CrossRef] [Green Version]
  11. Lei, X.; Noda, N.A. Analyses of dynamic response of vehicle and track coupling system with random irregularity of track vertical profile. J. Sound Vib. 2002, 258, 147–165. [Google Scholar] [CrossRef]
  12. Kouroussis, G.; Gazetas, G.; Anastasopoulos, I.; Conti, C.; Verlinden, O. Discrete modelling of vertical track-soil coupling for vehicle-track dynamics. Soil Dyn. Earthq. Eng. 2011, 31, 1711–1723. [Google Scholar] [CrossRef]
  13. Nguyen, K.; Goicolea, J.M.; Galbadon, F. Comparison of dynamic effects of high-speed traffic load on ballasted track using a simplified two-dimensional and full three-dimensional model. Proc. Inst. Mech. Eng. Part F J. Rail Rapid Transit 2014, 228, 128–142. [Google Scholar] [CrossRef]
  14. Chen, Z.W.; Zhai, W.M.; Cai, C.B.; Sun, Y. Safety threshold of high-speed railway pier settlement based on train-track-bridge dynamic interaction. Sci. China Technol. Sci. 2015, 58, 202–210. [Google Scholar] [CrossRef]
  15. Varandas, J.N. Long-Term Behavior of Railway Transitions under Dynamic Loading Application to Soft Soil Sites. Ph.D. Thesis, Faculdade de Ciências e Tecnologia, Universidade Nova de Lisboa, Caparica, Portugal, 2013. [Google Scholar]
  16. Varandas, J.N.; Hölscher, P.; Silva, M.A.G. Settlement of ballasted track under traffic loading: Application to transition zones. Proc. Inst. Mech. Eng. Part F J. Rail Rapid Transit 2014, 228, 242–259. [Google Scholar] [CrossRef]
  17. Mauer, L. An interactive track-train dynamic model for calculation of track error growth. Veh. Syst. Dyn. 1995, 24, 209–221. [Google Scholar] [CrossRef]
  18. Esveld, C.; de Man, A. Use of railway track vibration behavior for design and maintenance. In IABSE Symposium Report; IABSE Conference Rotterdam 2013, Assessment, Upgrading and Refurbishment of Infrastructures; International Association for Bridge and Structural Engineering: New York, NY, USA, 2013; Volume 99. [Google Scholar] [CrossRef]
  19. Kaewunruen, S.; Remennikov, A. Monitoring structural degradation of rail pads in laboratory using impact excitation technique. In Proceedings of the 1st International Conference on Structural Condition Assessment, Monitoring, and Improvement, Perth, Australia, 12–14 December 2005; Available online: https://ro.uow.edu.au/engpapers/285 (accessed on 20 February 2021).
  20. Azoh, T.S.; Nzie, W.; Djeumako, B.; Fotsing, B.S. Modeling of train track vibrations for maintenance perspectives: Application. Eur. Sci. J. 2014, 10, 260–275. [Google Scholar]
  21. Sadri, M.; Lu, T.; Steenbergen, M. Railway track degradation: The contribution of a spatially variant support stiffness—Local variation. J. Sound Vib. 2019, 455, 203–220. [Google Scholar] [CrossRef]
  22. Sadri, M.; Lu, T.; Steenbergen, M. Railway track degradation: The contribution of a spatially variant support stiffness—Global variation. J. Sound Vib. 2020, 464, a114992. [Google Scholar] [CrossRef]
  23. Markine, V.; Steenbergen, M.; Shevtsov, I. A dynamic model for analysis of damage of railway switches. In Proceedings of the International Symposium on Dynamic of Vehicles on Roads and Tracks; (IAVSD 2009), KTH, Stockholm, Sweden, 17–21 August 2009. [Google Scholar]
  24. Markine, V.; Steenbergen, M.; Shevtsov, I. Combatting RCF on switch points by tuning elastic track properties. Wear 2011, 271, 158–167. [Google Scholar] [CrossRef]
  25. Tejada, S.D.; Danielsen, H. Track Degradation and Maintenance. In Proceedings of the RailCPH/Banekonference 2017, Copenhagen, Denmark, 15 May 2017. [Google Scholar]
  26. Igeland, A. Railhead corrugation growth explained by dynamic interaction between track and bogie wheelsets. Proc. Inst. Mech. Eng. Part F J. Rail Rapid Transit 1996, 210, 11–20. [Google Scholar] [CrossRef]
  27. Igeland, A.; Ilias, H. Rail head corrugation growth predictions based on non-linear high frequency vehicle/track interaction. Wear 1997, 213, 90–97. [Google Scholar] [CrossRef]
  28. Ilias, H. The influence of railpad stiffness on wheelset/track interaction and corrugation growth. J. Sound Vib. 1999, 227, 935–948. [Google Scholar] [CrossRef]
  29. Jin, X.S.; Wen, Z.F. Effect of discrete track support by sleepers on rail corrugation at a curved track. J. Sound Vib. 2008, 315, 279–300. [Google Scholar] [CrossRef]
  30. Zhang, S.; Xiao, X.; Wen, Z.; Jin, X. Effect of unsupported sleepers on wheel/rail normal load. Soil Dyn. Earthq. Eng. 2008, 28, 662–673. [Google Scholar] [CrossRef]
  31. Zhu, J.Y.; Thompson, D.J.; Jones, C.J.C. On the effect of unsupported sleepers on the dynamic behavior of a railway track. Veh. Syst. Dyn. 2011, 49, 1389–1408. [Google Scholar] [CrossRef]
  32. Dong, R.; Sankar, S.; Dukkipati, R.V. A finite element model of railway track and its application to the wheel flat problem. Proc. Inst. Mech. Eng. Part F J. Rail Rapid Transit 1994, 208, 61–72. [Google Scholar] [CrossRef]
  33. Zhai, W.M.; Wang, Q.C.; Lu, Z.W.; Wu, X.S. Dynamic effects of vehicles on tracks in the case of raising train speeds. Proc. Inst. Mech. Eng. Part F J. Rail Rapid Transit 2001, 215, 125–135. [Google Scholar] [CrossRef]
  34. Wu, T.X.; Thompson, D.J. A hybrid model for the noise generation due to railway wheel flats. J. Sound Vib. 2002, 251, 115–139. [Google Scholar] [CrossRef]
  35. Nielsen, J.C.O.; Oscarsson, J. Simulation of dynamic train-track interaction with state-dependent track properties. J. Sound Vib. 2004, 275, 515–532. [Google Scholar] [CrossRef]
  36. Uzzal, R.U.A.; Ahmed, W.; Rakheja, S. Dynamic analysis of railway vehicle-track interactions due to wheel flat with a pitch-plane vehicle model. J. Mech. Eng. 2008, 39, 86–94. [Google Scholar] [CrossRef] [Green Version]
  37. Liu, X.; Zhai, W. Analysis of vertical dynamic wheel/rail interaction caused by polygonal wheels on high-speed trains. Wear 2014, 314, 282–290. [Google Scholar] [CrossRef]
  38. Thompson, D.J.; Hemsworth, B.; Vincent, N. Experimental validation of the TWINS prediction program for rolling noise, part 1: Description of the model and method. J. Sound Vib. 1996, 193, 123–135. [Google Scholar] [CrossRef]
  39. Thompson, D.J.; Fodiman, P.; Mahé, H. Experimental validation of the TWINS prediction program for rolling noise, part 2: Results. J. Sound Vib. 1996, 193, 137–147. [Google Scholar] [CrossRef]
  40. Wu, T.X.; Thompson, D.J. A double Timoshenko beam model for vertical vibration analysis of railway track at high frequencies. J. Sound Vib. 1999, 224, 329–348. [Google Scholar] [CrossRef]
  41. Wu, T.X.; Thompson, D.J. The effects of local preload on the foundation stiffness and vertical vibration of railway track. J. Sound Vib. 1999, 219, 881–904. [Google Scholar] [CrossRef]
  42. Wu, T.X.; Thompson, D.J. Vibration analysis of railway track with multiple wheels on the rail. J. Sound Vib. 2001, 239, 69–97. [Google Scholar] [CrossRef]
  43. Ahlbeck, D.R.; Meacham, H.C.; Prause, R.H. The development of analytical models for railroad track dynamics. Railr. Track Mech. Technol. 1975, 239–263. [Google Scholar] [CrossRef]
  44. Doyle, N.F. Railway Track Design: A Review of Current Practice; Bureau of Transport Economics (BTE), 1980. Available online: https://www.bitre.gov.au/sites/default/files/op_035.pdf (accessed on 20 February 2021).
  45. Kerr, D. On the determination of the rail support modulus k. Int. J. Solids Struct. 2000, 37, 4335–4351. [Google Scholar] [CrossRef]
  46. Paixão, A.L.M. Transition Zones in Railway Tracks—An Experimental and Numerical Study on the Structural Behavior. 2014. Available online: http://repositorio.lnec.pt:8080/jspui/handle/123456789/1007244 (accessed on 20 February 2021).
  47. Dimitrovová, Z. Critical velocity of a uniformly moving load on a beam supported by a finite depth foundation. J. Sound Vib. 2016, 366, 325–342. [Google Scholar] [CrossRef]
  48. Dimitrovová, Z. Analysis of the critical velocity of a load moving on a beam supported by a finite depth foundation. Int. J. Solids Struct. 2017, 122–123, 128–147. [Google Scholar] [CrossRef]
  49. Birmann, F. Track parameters, static and dynamic. Proc. Inst. Mech. Eng. Conf. Proc. 1965, 180, 73–85. [Google Scholar] [CrossRef]
  50. Newton, S.G.; Clark, R.A. An investigation into the dynamic effects on the track of wheelflats on railway vehicles. J. Mech. Eng. Sci. 1979, 21, 287–297. [Google Scholar] [CrossRef]
  51. Sato, Y. Recent developments in railway track mechanics in JNR. Jpn. Railw. Eng. 1981, 20, 1–5. [Google Scholar]
  52. ANSYS Inc. 12.1 Documentation Release 12.1 by ANSYS; ANSYS Inc.: Canonsburg, PA, USA, 2009. [Google Scholar]
  53. EN 13674-1. Railway Applications—Track—Rail—Part 1: Vignole Railway Rails 46 kg/m and above; Norm. CEN, European Committee for Standardization: Brussels, Belgium, 2011. [Google Scholar]
  54. IMV-019. Fabrico e Fornecimento de Travessas Monobloco UIC54 e 60; Norm 1000-2011035.02; Rede Ferroviária Nacional, REFER EPE: Lisbon, Portugal, 2000. (In Portuguese) [Google Scholar]
  55. Celep, Z.; Bazant, Z.P. Spurious reflection of elastic waves due to gradually changing finite element size. Int. J. Numer. Methods Eng. 1983, 19, 631–646. [Google Scholar] [CrossRef] [Green Version]
  56. Mednikov, A. Investigation of the depth of the dynamically active zone of soil and some of its characteristics. Soil Mech. Found. Eng. 1965, 2, 338–339. [Google Scholar] [CrossRef]
  57. Li, D.; Selig, E.T. Evaluation and Remediation of Potential Railway Subgrade Problems; Tech. Rep.; Association of American Railroads, Chicago Technical Center: Washington, DC, USA, 1995. [Google Scholar]
  58. Teixeira, P.F.; Ferreira, P.A.; López-Pita, A.; Casas, C.; Bachiller, A. The use of bituminous subballast on future high-speed lines in Spain: Structural design and economical impact. Int. J. Railw. 2009, 2, 1–7. Available online: http://www.ijr.or.kr/on_line/admin/files/vol.2_no.1_01.pdf (accessed on 20 February 2021).
  59. Indraratna, B.; Nimbalkar, S.; Christie, D.; Rujikiatkamjorn, C.; Vinod, J. Field assessment of the performance of a ballasted rail track with and without geosynthetics. J. Geotech. Geoenvironmental Eng. 2010, 136, 907–917. [Google Scholar] [CrossRef]
  60. Bronsert, J.; Baeßler, M.; Cuéllar, P.; Rücker, W. Numerical modeling of train-track-interaction at bridge transition zones considering the long-term behavior. In Proceedings of the 11th International Conference on Vibration Problems (ICOVP-2013), Lisbon, Portugal, 9–12 September 2013; Available online: http://www2.dec.fct.unl.pt/~zuzana/icovp/components/com_breezingforms/uploads/246_paper0.pdf (accessed on 20 February 2021).
  61. Bronsert, J.; Baeßler, M.; Cuéllar, P.; Rücker, W. Assessment and optimisation of bridge transition zones on the basis of a numerical model for train-track-bridge interaction. In Proceedings of the 9th International Conference on Structural Dynamics (EURODYN 2014), Porto, Portugal, 30 June–2 July 2014; Available online: https://paginas.fe.up.pt/~eurodyn2014/CD/papers/109_MS03_ABS_1298.pdf (accessed on 20 February 2021).
  62. Hall, L. Simulations and analyses of train-induced ground vibrations in finite element models. Soil Dyn. Earthq. Eng. 2003, 23, 403–413. [Google Scholar] [CrossRef]
  63. Zhai, W.; He, Z.; Song, X. Prediction of high-speed train induced ground vibration based on train-track-ground system model. Earthq. Eng. Eng. Vib. 2010, 9, 545–554. [Google Scholar] [CrossRef]
  64. Celebi, E.; Göktepe, F. Non-linear 2-D FE analysis for the assessment of isolation performance of wave impeding barrier in reduction of railway-induced surface waves. Constr. Build. Mater. 2012, 36, 1–13. [Google Scholar] [CrossRef]
  65. Kaynia, M.; Madshus, C.; Zackrisson, P. Ground vibration from high-speed trains: Prediction and countermeasure. J. Geotech. Geoenvironmental Eng. 2000, 126, 531–537. [Google Scholar] [CrossRef]
  66. Rodrigues, F.S. Viability and Applicability of Simplified Models for the Dynamic Analysis of Ballasted Railway Tracks. Ph.D. Thesis, Faculdade de Ciências e Tecnologia da Universidade Nova de Lisboa, Caparica, Portugal, 2017. [Google Scholar]
  67. Thompson, J.; Verheij, J.W. The dynamic behavior of rail fasteners at high frequencies. Appl. Acoust. 1997, 52, 1–17. [Google Scholar] [CrossRef] [Green Version]
  68. Deeks, J.; Randolph, M.F. Axisymmetric time-domain transmitting boundaries. J. Eng. Mech. 1994, 120, 25–42. [Google Scholar] [CrossRef]
  69. Liu, J.; Du, Y.; Du, X.; Wang, Z.; Wu, J. 3D viscous-spring artificial boundary in time domain. Earthq. Eng. Eng. Vib. 2006, 5, 93–102. [Google Scholar] [CrossRef]
  70. Jesus, H.; Dimitrovová, Z.; Silva, M.A.G. A Statistical Analysis of the Dynamic Response of a Railway Viaduct. Eng. Struct. 2014, 71, 244–259. [Google Scholar] [CrossRef] [Green Version]
  71. Lysmer, J.; Kuhlemeyer, R. Finite dynamic model for infinite media. J. Eng. Mech. Div. 1969, 95, 859–877. [Google Scholar] [CrossRef]
  72. Holland, J.H. Adaptation in Natural and Artificial Systems—An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence; reprint MIT Press, 1992 ed.; University of Michigan Press: Ann Arbor, MI, USA, 1975. [Google Scholar]
  73. Mylonakis, G.; Nikolaou, S.; Gazetas, G. Footings under seismic loading: Analysis and design issues with emphasis on bridge foundations. Soil Dyn. Earthq. Eng. 2006, 26, 824–853. [Google Scholar] [CrossRef]
  74. Jones, R.; Xenophontos, J. The Vlasov foundation model. Int. J. Mech. Sci. 1977, 19, 317–323. [Google Scholar] [CrossRef]
  75. Lysmer, J. Vertical Motion of Rigid Footings. 1965. Available online: https://vulcanhammernet.files.wordpress.com/2017/01/lysmer.pdf (accessed on 20 February 2021).
Figure 1. DSM, based on [4].
Figure 1. DSM, based on [4].
Vibration 04 00013 g001
Figure 2. Stress cone load distribution in the ballast, adapted from [3]: (a) without superposition; (b) with superposition.
Figure 2. Stress cone load distribution in the ballast, adapted from [3]: (a) without superposition; (b) with superposition.
Vibration 04 00013 g002
Figure 3. 3D FE model of the railway track.
Figure 3. 3D FE model of the railway track.
Vibration 04 00013 g003
Figure 4. Summary of ballast properties over the literature: (a) Young’s modulus, (b) Poisson’s ratio, (c) density.
Figure 4. Summary of ballast properties over the literature: (a) Young’s modulus, (b) Poisson’s ratio, (c) density.
Vibration 04 00013 g004
Figure 5. Vertical displacement of the rail: comparison between experimental measurements, fitted numerical results (2D) from [46] and numerical results (3D) obtained here; (a) full vehicle; (b) first bogie.
Figure 5. Vertical displacement of the rail: comparison between experimental measurements, fitted numerical results (2D) from [46] and numerical results (3D) obtained here; (a) full vehicle; (b) first bogie.
Vibration 04 00013 g005
Figure 6. Normalized vertical displacement of the rail for the 3D FE model and the DSM for v = 100   m / s , h b = 0.6   m , h s = 6   m and different combinations of E b and E s ; x F is the point of load application.
Figure 6. Normalized vertical displacement of the rail for the 3D FE model and the DSM for v = 100   m / s , h b = 0.6   m , h s = 6   m and different combinations of E b and E s ; x F is the point of load application.
Vibration 04 00013 g006
Figure 7. Static vertical stress in the 3D FE model of the railway track. Blue: maximum negative stress in the ballast; red: 0% to 10% of the maximum negative stress; grey: values out of range; h b = 0.6   m .
Figure 7. Static vertical stress in the 3D FE model of the railway track. Blue: maximum negative stress in the ballast; red: 0% to 10% of the maximum negative stress; grey: values out of range; h b = 0.6   m .
Vibration 04 00013 g007
Figure 8. Superposition of the stress cones in the transversal direction.
Figure 8. Superposition of the stress cones in the transversal direction.
Vibration 04 00013 g008
Figure 9. Geometry of the stress cone with superposition in both directions.
Figure 9. Geometry of the stress cone with superposition in both directions.
Vibration 04 00013 g009
Figure 10. Cross-section of the stress cone in the transversal direction.
Figure 10. Cross-section of the stress cone in the transversal direction.
Vibration 04 00013 g010
Figure 11. (a) K ¯ s and (b) K ¯ s , P as a function of h s for different values of γ .
Figure 11. (a) K ¯ s and (b) K ¯ s , P as a function of h s for different values of γ .
Vibration 04 00013 g011
Figure 12. Optimum values from Table 2 (crosses) and proposed expression (blue line) when h b = 0.3   m : (a) f K , s ; (b) f K , w , s .
Figure 12. Optimum values from Table 2 (crosses) and proposed expression (blue line) when h b = 0.3   m : (a) f K , s ; (b) f K , w , s .
Vibration 04 00013 g012
Figure 13. Vertical displacement of the rail for the 3D FE model and the DSM for h b = 0.6   m , h s = 6   m and different combinations of E b and E s .
Figure 13. Vertical displacement of the rail for the 3D FE model and the DSM for h b = 0.6   m , h s = 6   m and different combinations of E b and E s .
Vibration 04 00013 g013
Figure 14. Normalized vertical displacement of the rail for the 3D FE model and the DSM, v = 100   m / s ; h b = 0.3   m (a) stiff ballast and subgrade; (b) soft ballast and subgrade ( x F is the point of load application).
Figure 14. Normalized vertical displacement of the rail for the 3D FE model and the DSM, v = 100   m / s ; h b = 0.3   m (a) stiff ballast and subgrade; (b) soft ballast and subgrade ( x F is the point of load application).
Vibration 04 00013 g014
Figure 15. Vertical displacement of the rail: comparison between experimental measurements, the 3D FE model and the DSM using the mechanistic expressions; (a) full vehicle; (b) first bogie.
Figure 15. Vertical displacement of the rail: comparison between experimental measurements, the 3D FE model and the DSM using the mechanistic expressions; (a) full vehicle; (b) first bogie.
Vibration 04 00013 g015
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
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
Vibration 04 00013 g016
Figure 17. Steady-state vertical displacement of the rail: v = 100   m / s . Comparison of the model proposed in this paper, the model from [3] and the 3D FE model.
Figure 17. Steady-state vertical displacement of the rail: v = 100   m / s . Comparison of the model proposed in this paper, the model from [3] and the 3D FE model.
Vibration 04 00013 g017
Figure 18. Maximum downward normalized displacement and normalized displacement under the force as a function of velocity. Comparison of the model proposed in this paper and in [3].
Figure 18. Maximum downward normalized displacement and normalized displacement under the force as a function of velocity. Comparison of the model proposed in this paper and in [3].
Vibration 04 00013 g018
Table 1. Elastic material properties for the three-dimensional railway model.
Table 1. Elastic material properties for the three-dimensional railway model.
PropertyRailSleeperBallastSubgrade
Young’s modulus [MPa]210∙10338∙10350,100,15050,100,150
Poisson’s ratio0.30.20.250.35
Density [kg/m3]78502064 117501900
1 Modified by correction factor.
Table 2. Optimum values for the static case in combined optimizations.
Table 2. Optimum values for the static case in combined optimizations.
h s = { 6 , 25 , 50 }   m h s = { 3 , 6 , 9 }   m
Parameter h b = 0.3   m Parameter h b = 0.3   m Parameter h b = 0.3   m
f K , b [m]1.9111.483 f K , b [m]1.9561.349
f K , s , 6 [m]0.2150.221 f K , s , 3 [m]0.3290.346
f K , s , 25 [m]0.1510.137 f K , s , 6 [m]0.2090.225
f K , s , 50 [m]0.1510.140 f K , s , 9 [m]0.1730.183
f K , w , b [m]0.2500.722 f K , w , b [m]0.1090.499
f K , w , s , 6 [m]3.6743.751 f K , w , s , 3 [m]2.8073.856
f K , w , s , 25 [m]4.9206.398 f K , w , s , 6 [m]4.1755.202
f K , w , s , 50 [m]5.0336.906 f K , w , s , 9 [m]5.0786.872
Φ 2.4%–5.7%3.4%–5.5% Φ 1.7%–4.9%1.7%–6.5%
Table 3. Optimum values for the dynamic case in combined optimizations.
Table 3. Optimum values for the dynamic case in combined optimizations.
v = 50   m / s v = 100   m / s
v = 50   m / s v = 100   m / s v = 50   m / s v = 100   m / s v = 50   m / s v = 100   m / s
f C , r a d , s [m]0.3870.418 f C , r a d , s [m]0.3530.436
f M , b [m3]1.8762.521 f M , b [m3]1.4111.734
M b [t]3.2834.412 M b [t]2.4693.035
Φ 4.6%–8.2%5.2%–11.5% Φ 7.0%–18.5%8.5%–19.0%
Table 4. Optimized value of α b , f K , b [m] and respective error, optimum values from Table 2.
Table 4. Optimized value of α b , f K , b [m] and respective error, optimum values from Table 2.
h b [ m ] Opt. Val. α b = 49.8 °ErrorOpt. Val. α b = 47.6 °Error
f K , b f K , b f K , b f K , b
0.31.9112.086+9%1.9562.029+4%
0.61.4831.344−9%1.3491.302−4%
Table 5. Value of f K , w , b [m] and respective error, optimum values from Table 2.
Table 5. Value of f K , w , b [m] and respective error, optimum values from Table 2.
h b [ m ] Opt. Val. α b = 49.8 °ErrorOpt. Val. α b = 47.6 °Error
f K , w , b f K , w , b f K , w , b f K , w , b
0.30.2500.226−9%0.1090.222+104%
0.60.7220.949+31%0.4990.930+86%
Table 6. Values of the mass and respective error, optimum values from Table 3.
Table 6. Values of the mass and respective error, optimum values from Table 3.
h b [ m ] Opt. Val.
v = 50   m / s
α b = 49.8 °Error Opt. Val.
v = 100   m / s
α b = 49.8 °Error
M M M M
0.33.2833.273−3%2.4693.273+33%
0.64.4124.900+11%3.0354.900+61%
Table 7. Optimum values of c z , based on Table 3.
Table 7. Optimum values of c z , based on Table 3.
h b [ m ] v = 50   m / s v = 100   m / s
0.30.3930.358
0.60.3470.361
Table 8. Model data summary.
Table 8. Model data summary.
ParameterFormulas in This Paper[3]
K w , b [MN/m]27.378.4 *
K w , s [MN/m]500.9---
K w [MN/m]528.2
M b [kg]595.7531.4
M s [kg]3033.6---
M [kg]3629.3
K b [MN/m]168.27137.75
K f [MN/m]88.877.5
C f [kNs/m]30831.15 *
* Means that the value was used without justification.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Rodrigues, A.F.S.; Dimitrovová, Z. Applicability of a Three-Layer Model for the Dynamic Analysis of Ballasted Railway Tracks. Vibration 2021, 4, 151-174. https://0-doi-org.brum.beds.ac.uk/10.3390/vibration4010013

AMA Style

Rodrigues AFS, Dimitrovová Z. Applicability of a Three-Layer Model for the Dynamic Analysis of Ballasted Railway Tracks. Vibration. 2021; 4(1):151-174. https://0-doi-org.brum.beds.ac.uk/10.3390/vibration4010013

Chicago/Turabian Style

Rodrigues, André F. S., and Zuzana Dimitrovová. 2021. "Applicability of a Three-Layer Model for the Dynamic Analysis of Ballasted Railway Tracks" Vibration 4, no. 1: 151-174. https://0-doi-org.brum.beds.ac.uk/10.3390/vibration4010013

Article Metrics

Back to TopTop