Next Article in Journal
Exploring the Relationship between Muscular Strength, Flexibility, and mSEBT Test Performance in Saudi Arabian Women
Next Article in Special Issue
Power Dissipation and Wear Modeling in Wheel–Rail Contact
Previous Article in Journal
Meat Analogues from Pea Protein: Effect of Different Oat Protein Concentrates and Post Treatment on Selected Technological Properties of High-Moisture Extrudates
Previous Article in Special Issue
The Effect of Cable Aging on Surge Arresters Designed by Genetic Algorithm
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Instability of Vibrations of Mass(es) Moving Uniformly on a Two-Layer Track Model: Parameters Leading to Irregular Cases and Associated Implications for Railway Design

by
Zuzana Dimitrovová
1,2
1
Department of Civil Engineering, NOVA School of Science and Technology, NOVA University of Lisbon, 2829-516 Caparica, Portugal
2
Institute of Mechanical Engineering (IDMEC), Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal
Appl. Sci. 2023, 13(22), 12356; https://doi.org/10.3390/app132212356
Submission received: 10 September 2023 / Revised: 7 November 2023 / Accepted: 9 November 2023 / Published: 15 November 2023
(This article belongs to the Special Issue Railway Dynamic Simulation: Recent Advances and Perspective)

Abstract

:
Ballasted railway tracks can be modeled using reduced/simplified models composed of several layers of discrete components. This paper deals with the two-layer model, which is very popular due to its computational efficiency. In order to provide some recommendations for track design, it is necessary to identify which set of parameters leads to some irregular/unexpected behavior. In this paper, irregularities are investigated at three levels, namely, (i) the critical velocity of a moving constant force, (ii) the instability of one moving mass, and (iii) the instability of two moving masses. All results are presented in a dimensionless form to cover a wide range of real parameters. Irregular cases are identified by sets of parameters leading to them, which is the main finding of this paper; then, general conclusions are drawn. Regarding the method, all results are obtained analytically or semi-analytically, where “semi” refers to solving the roots of a given polynomial using predefined numerical procedures in symbolic software. No numerical integration is involved in any of the results presented. This means that the results are highly accurate and refer to exact values, so any kind of parametric or sensitivity analyses is readily possible.

1. Introduction

The two-layer railway track model is widely used by other researchers to approximate several phenomena. Here, it will be used to identify the set of parameters under which it exhibits irregular behavior, i.e., behavior that is not expected.
In recent decades, a considerable amount of research has been presented in this field, proving that it is a field of great interest and still very active thanks to the current trend of decarbonization, which requires the transfer of road transport to rail. With this in mind, it is not possible to cover all relevant research in this section. Emphasis will be placed on railway-related research implementing similar solution techniques, avoiding finite and boundary element methods in this review.
Published works can be classified according to several criteria. First, a separation can be made according to the structure type, into finite or infinite, or according to moving objects, into structures subjected to moving force(s) or moving inertial object(s) such as mass(es), oscillator(s), or even simplified vehicle models. After that, a separation can be made according to the supporting structure arrangement, into continuous (2D or 3D) or discrete, which usually consists of several layers. A further separation can distinguish between fully linear and non-linear behavior of the supporting structure components.
When building a literature review, it is customary to start with Frýba’s monograph [1], which includes several cases mentioned in the previous paragraph. The moving force problem is generally much simpler, and fully analytical solutions for infinite structures are given in [2,3]. The critical velocity is also determined, but in such structures with massless foundations, waves cannot propagate in the foundation, so the conclusions are unrealistic. The results are more realistic when the foundation mass is taken into account [4,5,6].
In general, and especially in railway applications, the beam is modeled according to the Euler–Bernoulli theory. Some works are devoted to the analysis of the differences induced in the results by the implementation of linear and non-linear beam theories. This is exemplified on an elastically supported finite beam without an elastic foundation traversed by a moving mass in [7,8], where a similar analysis is made and results obtained using different beam theories are compared. A multi-span beam is analyzed in [9]. The main focus is on the mechanical behavior of the beam, modeled with linear viscoelastic material according to the Kelvin–Voight model, which is suitable for concrete structures but not for applications where the beam is represented by a steel rail. In all these works, besides time discretization, spatial discretization is also used, followed by a numerical solution of the governing equations, which is somewhat beyond the scope of this overview. Although parametric analyses are conducted to extract the maximum deflections, the issue of instability is not addressed, and since there is no elastic foundation, the conclusions drawn are not directly relevant to the analysis presented in this paper.
Pioneering works on the instability of moving inertial objects include [10,11,12]. In [10,11], several masses moving on a finite beam are considered, and in [12], the problem of one mass moving on a viscoelastically supported infinite beam is solved. The problem of one moving mass is further described in detail in [13,14], where the instability is determined using the D-decomposition method. The physical explanation of the phenomenon is related to the anomalous Doppler waves radiated by the moving object. Mathematically, such a situation is detected due to the exponential increase in vibration amplitudes, implying that the steady-state solution is never reached even in the presence of damping. Several moving inertial objects on finite structures have been recently analyzed in [15,16,17]. Other works about infinite structures have also implemented the D-decomposition method together with the dynamic Green’s function [18,19,20,21] or integral transforms [22,23]. It is commonly assumed that the mass is in permanent contact with the beam [22,23]; however, in some works, a contact spring is introduced [18,19,20,21]. The critical velocity of the moving force is then usually not analyzed. A connection with the critical velocity and a different approach to identifying instability by tracing the instability lines is presented in [24,25]. This approach is also suitable for the problem of two moving proximate masses, where a strong dynamic interaction can significantly alter the onset of instability. Using the mentioned approach, the conditions under which the results can be superposed can also be derived, and cases where the dynamic interaction induces instability at a velocity lower than the lowest critical velocity of the moving force can be identified. Additionally, this approach can be readily extended to moving oscillators.
Recent works deal with more complicated moving objects or foundation [26,27,28]. Further extensions to large deflections are considered in [29]. The change in stiffness of the Winkler foundation introduced with a smooth variation is analyzed in [30,31]. A moving force is applied in [30], while a moving oscillator is applied in [31]. A comparison of possible methods for solving similar problems is presented in [32]. When inhomogeneity is considered, the moving element method, successfully applied in [33,34,35,36], or the moving window [37] should also be mentioned, because they are efficient in dealing with such situations. In [36], a computationally efficient and accurate numerical method for the dynamic response of a maglev train passing on an “infinitely” long multi-span bridge is proposed. The results are validated against test cases from the literature. A novel model for the dynamic interaction of a beam with its foundation is introduced in [38]. It is shown that the Winkler and Pasternak moduli are not constant but time-dependent, because they are influenced by the beam–foundation interaction. Therefore, they are determined as a part of the solution. The mass that is dynamically activated in the foundation is determined in a similar way.
The two-layer model of the railway track is quite popular and has already been used to accomplish several objectives, like, e.g., in [39,40]. The three-layer model is more realistic, and its applicability to represent full finite element models is dealt with in [41] where the question of the critical velocity of the moving force is also briefly discussed. Among other recent works on layered models, Dimitrovová [42] can be mentioned. A very long finite three-layer model is considered in [43]. The modal expansion method is applied to a reduced model to increase computational efficiency. A control volume is used to reduce the computational domain, and the structure is traversed by a multi-body model of the vehicle.
Recently, artificial intelligence, machine learning and several other metaheuristic algorithms are gaining popularity and have been successfully implemented in identification of the critical velocity of a railway track [44]. Such an approach is fully justified by the fact that the track corridor includes several track sections, and a separate analysis of each of them would be computationally demanding. However, the present paper is dealing with a much simpler situation where the track section is longitudinally homogeneous, and therefore, it is possible to derive (semi)analytical results where all parameter influences are directly present. Therefore, there is no need to seek help in predictions by metaheuristic approaches.
In this paper, the two-layer model is analyzed from a different point of view. The goal is to establish the range of parameters leading to an irregular behavior, where “irregular” means “different than expected”. The model is introduced in Section 2, where three irregular situations are described in detail. Furthermore, the range of allowable parameters is identified in Section 3. After that, three sections are devoted to the three situations where some irregularity occurs. Finally, some conclusions are drawn in Section 7.
There are no published works that deal so widely with the identified irregularities and summarize the results over a wide range of admissible parameters, which identifies the novel contribution of this paper. Such analysis has a practical importance because unexpected behavior can complicate the railway track design. Some results concerning the two-layer model have already been included in [25,42]. However, [42] is focused on the comparison of layered models, while [25] is specifically dedicated to the identification of full beam vibrations by implementing a new form of results presentation. The present paper has a different focus, and the main objective is to identify sets of parameters leading to unexpected results, which is a problem that has not been analyzed in [25,42] or elsewhere.

2. The Model: Governing Equations and Irregular Situations

The track is modeled as an infinite beam standing for the rail, supported by linear spring–damper elements representing the rail pads, k p and c p , point masses modeling the sleeper half-masses, M s , and supported further by linear spring–damper elements representing the vertical stiffness and damping of all other layers, k f and c f , that is the ballast, sub-ballast, foundation and other layers that might be placed below the sleepers. The rail may be subjected to an axial force N . The model is shown in Figure 1.
The main assumptions for the analysis are listed below:
  • The beam is prismatic and straight, and its material is homogeneous and isotropic.
  • The beam may be subjected to an axial force acting on its axis (considered positive when inducing compression).
  • The beam obeys linear elastic Euler–Bernoulli theory. Other researchers have found this beam theory adequate for modeling the steel rail because the beam height is small with respect to the length affected by the main vibrations. Differences in results obtained by other beam theories are noticeable only at very high frequencies [45]. In addition, the main simplification of the two-layer model affects the foundation model, and therefore, additional details regarding higher-order beam theories do not provide more realistic results.
  • The vertical displacements are measured from the equilibrium position determined by the weight of the model’s components.
  • The initial conditions are assumed to be homogeneous; however, this does not affect either the critical velocity or the conditions for instability.
  • Mass horizontal position is determined by its velocity.
  • No friction is assumed to act at the contact point.
  • The load and vertical displacements are considered positive downward.
  • As usual in several applications, the acting force may or may not represent the moving mass weight.
  • It is assumed that the moving masses are always in contact with the beam, which means that the contact is assumed as tight, and, therefore, the mass displacement and the displacement of the corresponding beam axis position are the same at all times. An extension to a situation with a linear contact spring can be easily implemented, but it is not presented in this paper.
The equations of motion of the model depicted in Figure 1 in fixed coordinates are
E I w , x x x x + N w , x x + m w , t t + c p ( w , t u , t ) + k p ( w u ) = p ( x , t )
M s u , t t c p ( w , t u , t ) k p ( w u ) + k f u k s u , x x + c f u , t = 0
where x is the spatial coordinate and t is the time. In addition to the parameters defined in the first paragraph of this section, E I and m are the bending stiffness and mass per unit length of the beam, w ( x , t ) and u ( x , t ) are the unknown vertical displacement fields of the beam and discrete masses, respectively, and p ( x , t ) means the load. Derivatives are denoted by the respective variables in subscripts preceded by a comma. The load is given by
p ( x , t ) = ( P 1 M 1 w 01 , t t ( t ) ) δ ( x x 1 ) + ( P 2 M 2 w 02 , t t ( t ) ) δ ( x x 2 )
where δ is the Dirac delta function, P 1 and P 2 are the moving forces, M 1 and M 2 are the moving masses, v is the constant velocity and d = x 2 x 1 is the distance between the moving masses. Without affecting generality too much, all parameters are assumed to be distributed to simplify the solution. Thus, discrete values are divided by the distance between the sleepers l s , which will result in continuous supports. In Figure 1, there is also a shear component k s because the model should have some shear stiffness. Placing it between the sleepers is not very realistic, but this is immaterial for the analysis. A different treatment must be implemented to make this value distributed, because the condition of uniform distribution must be applied on the inverse value. Therefore, if for a particular track a value of k ¯ s with unit N/m is used, then k s = k ¯ s l s with the correct unit of N.
In Equation (3), w 01 ( t ) and w 02 ( t ) are the masses displacements. Owing to the assumption of the tight contact, w 0 j ( t ) = w ( x j , t ) , j = 1 , 2 holds with x 1 = v t and x 2 = v t + d . Therefore, the relevant derivatives using the chain rule are
w 01 , t t ( t ) = v 2 w , x x ( x , t ) + 2 v w , x t ( x , t ) + w , t t ( x , t )   with   x = x 1
w 02 , t t ( t ) = v 2 w , x x ( x , t ) + 2 v w , x t ( x , t ) + w , t t ( x , t )   with   x = x 2
Subsequently, the loading term reads as follows:
( P 1 M 1 ( w , t t ( x , t ) + 2 v w , x t ( x , t ) + v 2 w , x x ( x , t ) ) ) δ ( x v t ) + ( P 2 M 2 ( w , t t ( x , t ) + 2 v w , x t ( x , t ) + v 2 w , x x ( x , t ) ) ) δ ( x v t d )
Nevertheless, the solution will be presented in moving coordinates r = v x t , t = t , leading to the following system of equations (compare with [24,25,42]):
E I w , r r r r + N w , r r + m ( w , t t 2 v w , r t + v 2 w , r r ) + c p ( w , t v w , r u , t + v u , r ) + k p ( w u ) = ( P 1 M 1 w , t t ) δ ( r ) + ( P 2 M 2 w , t t ) δ ( r d )
M s ( u , t t 2 v u , r t + v 2 u , r r ) c p ( w , t v w , r u , t + v u , r ) k p ( w u ) k s u , r r + k f u + c f ( u , t v u , r ) = 0
In the following text, the term “active point” will be used to spatially locate a moving object on a beam. The moving object can be either a force or a mass upon which the force acts. The designation “critical velocity” will be used for the critical velocity of a moving constant force. In this sense, it corresponds to a resonance effect because if the force is moving with its critical velocity in the absence of damping, then the steady-state solution is never achieved because the deflections tend to infinity and never stabilize. The critical velocity also indicates the separation between two distinct types of the steady-state beam shape. Up to the lowest critical velocity, approaching the lowest value from below, the maximum deflection of the steady-state solution over the entire beam is achieved at the active point. The minimum deflection (in absolute value) is significantly lower than the maximum. As already mentioned, there is no solution for the critical velocity. For higher velocities, zero deflection is achieved at the active point, and maximum and minimum deflections across the entire beam are the same. Similar properties can be found at higher critical velocity; in particular, the jump to zero deflection at the active point when crossing the critical velocity is preserved.
Since the model has two layers, it is expected to have two critical velocities. However, by analyzing the equations, it can be shown that there will always be one or three resonances. In the case of three resonances, after sorting them by numerical value, the value in the middle will identify the so-called “false critical velocity” [42]. This means that resonance will occur in the sense that a steady-state solution will never be reached, but the expected properties listed above will not be met. Additionally, it will be shown that the false critical velocity has no effect on stability. The regular case is therefore defined as the case with three resonances and the irregular case is therefore defined as the case with only one resonance. In the case of one resonance, the missing critical velocity has to be supplemented with the so-called “pseudocritical velocity”, which will always be lower than the critical velocity.
To access the instability conditions, either for one or two moving masses, instability lines will be traced in the velocity-moving mass plane. The approach presented here is simpler than the D-decomposition method used by other authors [13,14,18,19,20,21,22,23]. It simply means searching for the real-valued mass induced frequency that fulfills the characteristic equation in a damped case. This is because such real-valued mass induced frequency always marks the separation between the regions with a stable and unstable unsteady harmonic part of the full solution.
In the case of one moving mass, instability never occurs in the subcritical velocity range; however, it is necessary to highlight that in the case of only one resonance, the pseudocritical velocity must play the role of the lower critical velocity. In the regular case, there are only two instability line branches, which tend to infinite moving mass near (pseudo)critical velocities and to zero moving mass when the velocity is approaching infinity. The instability branches for lower damping are below those for higher damping and do not intersect. Moreover, the branches asymptotically approach (pseudo)critical velocities from the left in a better way than from the right. The irregular cases are thus characterized by a situation where the instability branches intersect the higher critical velocity. In such a case, there are usually more than two branches, and the additional ones can also be closed.
Unfortunately, in the case of two moving masses, the (pseudo)critical velocities do not provide any reasonable indication. There is an infinite number of instability branches, whose asymptotes are not linked to any of the (pseudo)critical velocities, and these branches generally intersect both (pseudo)critical velocities. The irregular situation is thus identified as the one where the dynamic interaction shifts the onset of instability into the subcritical range of velocities. However, since mathematically, the lower (pseudo)critical velocity is always intersected by an instability branch, this condition will be considered as an indication for the irregular case only when this intersection occurs for a realistic value of the moving mass.
It should be noted that the instability branches cannot terminate suddenly. A detailed analysis has shown that any branch always has a continuation until it is closed or tends asymptotically to infinite or zero moving mass. The asymptotic tendency has its origin in three factors: (i) the moving mass tends to infinity because the real-valued mass-induced frequency tends to zero; (ii) the moving mass tends to infinity because the equivalent flexibility of the supporting structure tends to zero; and (iii) the moving mass tends to zero because the velocity tends to infinity.

3. Range of Allowable Values

In the following discussion, the range of allowable values is specified. As for the rail, the range of possible values is quite narrow; basically, there are two guide sets of values for the 54E1 and 60E1 rails, determining the limits on E I and m . More variability can be attributed to other data. They are listed in Table 1 and correspond to the values given in [46,47,48].
The range of allowable parameters is defined with respect to their dimensionless counterparts, for which a reference Winkler beam is selected. First, the dimensionless moving spatial coordinate is introduced as
ξ = χ r   where   χ = k f 4 E I 4
The other values are related to the reference Winkler beam by
μ s = M s m ,   κ p = k p k f
η p = c p 2 m k f ,   η f = c f 2 m k f ,   η N = N 2 k f E I ,   η s = k s 2 k f E I
where for simplicity, the previous designations are kept, but all values (except for N ) are assumed as longitudinally distributed.
Therefore, the decisive parameters are μ s and κ p . Assuming that the sleepers are spaced from l s = 0.545 m [49,50] to l s = 0.711 m [50], this together with values in Table 1 identifies a relatively narrow interval for μ s ( 1 ; 5.5 ) , but it gives a very large variability for κ p ( 0.03 ; 40000 ) . A typical spacing in Portugal [51], in most European and some other countries, e.g., India, is 0.6 m. According to [50], l s should be related to the sleeper type and ballast height. There are studies analyzing how significant savings could be achieved by increasing l s [52]; nevertheless, these are mainly theoretical studies, and larger spacings than specified above are not yet common in practice. For practical reasons, a shorter distance is inadequate, having in mind the typical sleeper base.
There is little data on shear stiffness, so it will be set to zero in most examples. Non-zero values would make the model generally stiffer. Similarly, a compressive normal force would make the model softer and a tensile one stiffer. However, since this is not the main concern, it will also be set to zero in most examples. Regarding the viscous damping coefficients, very different values can be found in the literature, so it was decided to directly vary the dimensionless parameters.

4. Critical Velocity of a Constant Moving Force

To determine the critical velocity of a moving force, Equations (1)–(3) are first switched to moving coordinates and then changed to a dimensionless form. After removing the time-dependent terms not required for the steady-state regime, they read as
w ˜ , ξ ξ ξ ξ + 4 α 2 w ˜ , ξ ξ + 4 η N w ˜ , ξ ξ + 8 η p α ( w ˜ , ξ + u ˜ , ξ ) + 4 κ p ( w ˜ u ˜ ) = 8 η P 1 δ ( ξ )
μ s α 2 u ˜ , ξ ξ η s u ˜ , ξ ξ 2 η p α ( w ˜ , ξ + u ˜ , ξ ) κ p ( w ˜ u ˜ ) 2 η f α u ˜ , ξ + u ˜ = 0
where, besides the parameters specified in Equations (9)–(11),
w ˜ = w w s t ,   u ˜ = u w s t   with   w s t = P 1 χ 2 k f ,   and   η P 1 = P 1 P 1
v r e f = 4 k f E I m 2 4 = 1 χ k f m ,   α = v v r e f
Next, the Fourier transform is applied. In the transferred space, it is more convenient to write the equations in the matrix form to facilitate the analytical solution
[ D 1 + D 2 D 2 D 2 D 2 + D 3 ] { W U } = { 2 η P 1 0 }
where
D 1 ( p ) = p 4 4 α 2 p 2 η N p 2
D 2 ( p ) = 2 i η p α p + κ p
D 3 ( p ) = μ s α 2 p 2 2 i η f α p + η s p 2 + 1
and W , U and p represent w ˜ , u ˜ and ξ in the Fourier space. Since the critical velocities should be determined without damping, expressions (18) and (19) can be further simplified.
To analytically determine the critical velocity, it is necessary to find at least one real double p -root of the determinant of system (16). This means finding a real-valued solution to a system of two equations, where the first one is the determinant and the second one is its derivative with respect to p . The determinant can be written as a cubic polynomial with real coefficients for p ¯ = p 2 , and the second equation as a quadratic polynomial with real coefficients for p ¯
p ¯ 3 ( μ s α 2 η s ) p ¯ 2 ( 4 ( μ s α 2 η s ) ( η N + α 2 ) + κ p + 1 ) + 4 p ¯ ( α 2 ( κ p μ s + κ p + 1 ) + η N ( κ p + 1 ) η s κ p ) 4 κ p = 0
3 p ¯ 2 ( μ s α 2 η s ) 2 p ¯ ( 4 ( μ s α 2 η s ) ( η N + α 2 ) + κ p + 1 ) + 4 ( α 2 ( κ p μ s + κ p + 1 ) + η N ( κ p + 1 ) η s κ p ) = 0
In order to express directly α 2 , Equation (21), multiplied by p ¯ , should be subtracted from Equation (20) multiplied by 2. Then
α 2 = 8 κ p p ¯ 3 η s 4 p ¯ ( η N ( κ p + 1 ) η s κ p ) 4 p ¯ ( κ p μ s + κ p + 1 ) p ¯ 3 μ s
which is simplified when η N = η s = 0 to
α 2 = 8 κ p 4 p ¯ ( κ p μ s + κ p + 1 ) p ¯ 3 μ s
Thus, finding the real double p -root is an easy task using the predefined root-finding procedures in any symbolic software. Equation (22) or Equation (23) is substituted into Equation (21), which (for the given data κ p , μ s , η s and η N ) can be solved for valid p ¯ and consequently p and α .
The irregular case was defined as the case with only one resonance. This analysis is presented for the simpler case of η N = η s = 0 . After substituting Equation (23) into Equation (21), a cubic equation is obtained for p ¯ 2 = p 4 . Its discriminant is simplified to
μ s 3 κ p 3 ( κ p 8 ) + 3 μ s 2 κ p 2 ( κ p 2 ) 2 + 3 μ s κ p ( κ p 3 3 κ p 2 ) + ( κ p + 1 ) 4
The irregular case occurs when expression (24) is negative, which means that there is only one real root for p ¯ 2 = p 4 . Expression (24) is a fourth-order polynomial for κ p and third-order polynomial for μ s . To identify the separation between three and one resonance, μ s is solved analytically as a function of κ p .The analytical values of μ s are given by
θ = 1 3 arccos ( 2 B 1 3 9 B 1 B 2 + 27 B 3 2 ( B 1 2 3 B 2 ) 3 / 2 )
μ s j = B 1 3 + 2 3 B 1 2 3 B 2 cos ( θ + ( j 1 ) 2 π 3 ) , j = 1 , 2 , 3
where
B 1 = 3 ( κ p 2 ) 2 κ p ( κ p 8 )
B 2 = 3 ( κ p 3 3 κ p 2 ) κ p 2 ( κ p 8 )
B 3 = ( κ p + 1 ) 4 κ p 3 ( κ p 8 )
The separation of regions with three and only one resonance is shown in Figure 2, where irregular stands for one resonance, regular stands for three resonances and the cases falling on the boundary have, in fact, two resonances because one of the three resonances is doubled.
The curve in Figure 2 is obtained by plotting μ s 3 until κ p = 2 and then by μ s 1 for larger κ p . It has an asymptotic tendency to κ p = 8 , but this only occurs for μ s outside the allowable range. It can be observed that for μ s > 1 , there is always an open interval of κ p for which there is only one resonance. For μ s = 1 , there is only a double root κ p = 1 leading to a collapsed interval, and therefore, one of the three resonances is doubled. The same happens when the selection κ p and μ s lies on the curve indicated in Figure 2. In these cases, the critical velocity properties are suppressed, because the false critical velocity is decisive and the case is irregular.
The results of the parametric analyses of four selected cases are presented in Figure 3; in all of them, η N = η s = 0 and for better clarity, the scale does not start at α = 0 . As an example for the regular situation, μ s = 1 , κ p = 0.36 are selected, and for the irregular one: μ s = 1.1 , κ p = 0.6 . In the regular case, as shown in Figure 3a, it can be seen that the analytically determined values for the resonances 0.681, 0.707 and 0.913 indicate the lower critical velocity ( α c r 1 ), the false critical velocity ( α f c r ) and the higher critical velocity ( α c r 2 ), as predicted. In the irregular case, as shown in Figure 3b, it is observed that the only analytically determined value of 0.996 is α c r 2 and α c r 1 is replaced by the pseudocritical velocity ( α p c r ) with a value of 0.692. Other cases concern the situations at the boundary of the domain shown in Figure 2, i.e., when one of the roots is double. Figure 3c shows when μ s = 1 and κ p = 1 , which corresponds to the case where the interval of irregular cases is collapsed to a single value, proving that α c r 1 is suppressed by α f c r (both equal to 0.707) and the case is irregular with a unique α c r 2 of 1.125. The same happens in the case plotted in Figure 3d for μ s = 3 and κ p = 0.125 , which corresponds to the lower end of the irregular case interval, where α c r 1 is again suppressed by α f c r (both equal to 0.553) and the case is irregular with a unique α c r 2 of 0.658.
This analysis is consistent with correspondence with finite beams according to [53]. As explained in [42], α c r j corresponds to local minima in the resonant velocity plot, while α f c r corresponds to a local maximum, and the double root corresponds to the inflection point.
In conclusion, the irregular cases only mean that the first peak in displacements cannot be found by solving for the resonances, and therefore, α p c r must be determined by parametric analysis. This is important for two reasons. First, α c r may not provide the relevant information about the lowest velocity where excessive vibrations occur. Second, when it comes to instability, α p c r plays the role of α c r 1 and therefore must be known in advance to draw correct conclusions regarding the instability of moving inertial objects. In addition, α p c r can be dominant like in Figure 3b or not. When it is barely visible, extreme values are not always reached at the same α , and the determination of α p c r is ambiguous.
A very good estimate of α c r 1 is provided by extension of the classical formula by considering the involved stiffnesses in series and by summing the masses, i.e., by
α c r 1 = κ p ( 1 + κ p ) ( 1 + μ s ) 2 4 1 + η s η N
which is valid for higher κ p , approximately for κ p > 100 , which, in fact, means that for such values, α c r 1 is practically independent on κ p and thus
α c r 1 1 1 + μ s 1 + η s η N
can also be used. This identifies the influential parameters on α c r j and thus suggests how this value can be adapted to suit the needs. To complete the analysis, graphs of resonances and an indication of the irregular cases is further summarized in Figure 4, Figure 5, Figure 6 and Figure 7, also with the influence on η N and η s .
The results presented in Figure 4, Figure 5, Figure 6 and Figure 7 allow one to conclude how α c r 1 is affected by model parameters as well as how to change the model parameters to increase or decrease this value respecting the actual needs. The easiest way to alter parameter κ p is by changing the stiffness of the rail pads, μ s can be altered by the sleepers’ selection; additionally, the reference value k f can be adjusted by ballast layer height or by other reinforcing layers below it. The changes resulting from modifying the sleepers’ spacing are not very important, and introducing a different rail profile than that intended for such a railway is usually not possible.

5. Instability of One Moving Mass

To solve the problem of instability of one moving mass, the previously omitted time-dependent terms must be included back; thus, Equations (12) and (13) are now given by
w ˜ , ξ ξ ξ ξ + 4 η N w ˜ , ξ ξ + 4 α 2 w ˜ , ξ ξ + 4 w ˜ , τ τ 8 α w ˜ , ξ τ + 8 η p ( w ˜ , τ α w ˜ , ξ u ˜ , τ + α u ˜ , ξ ) + 4 κ p ( w ˜ u ˜ ) = ( 8 η P 1 4 η M 1 w ˜ , τ τ ) δ ( ξ )
μ s ( u ˜ , τ τ 2 α u ˜ , ξ τ + α 2 u ˜ , ξ ξ ) 2 η p ( w ˜ , τ α w ˜ , ξ u ˜ , τ + α u ˜ , ξ ) κ p ( w ˜ u ˜ ) η s u ˜ , ξ ξ + u ˜ + 2 η f ( u ˜ , τ α u ˜ , ξ ) = 0
where
τ = χ v r e f t ,   η M 1 = M 1 χ m
In the solution method, the Laplace transform must be applied first and then the Fourier one second. Initial conditions are assumed to be homogeneous. The following definition is used for the Laplace transform
F ˜ ( ξ , q ¯ ) = 0 f ( ξ , τ ) e q ¯ τ d τ
to keep its formalism, q ¯ = i q will be used in the following with q designating the frequency. In the transformed space, the equations can be again written in the matrix—formally the same as Equation (16) except for the right-hand side
[ D 1 + D 2 D 2 D 2 D 2 + D 3 ] { W U } = { 2 η P 1 i q + η M 1 q 2 W ˜ ( 0 , i q ) 0 }
and also
D 1 ( p , q ) = p 4 4 ( q α p ) 2 η N p 2
D 2 ( p , q ) = 2 i η p ( q α p ) + κ p
D 3 ( p , q ) = μ s ( q α p ) 2 + 2 i η f ( q α p ) + η s p 2 + 1
which shows that whenever there was α p in Equations (17)–(19) is now q α p .
In Equation (36), there is still W ˜ ( 0 , i q ) , which must be solved for. For this, first, Equation (36) is solved for W ( p , i q ) .
W ( p , i q ) = 1 D 1 ( p , q ) + D 2 ( p , q ) D 3 ( p , q ) D 2 ( p , q ) + D 3 ( p , q ) ( 2 η P 1 i q + η M 1 q 2 W ˜ ( 0 , i q ) )
Now, one can perform the inverse Fourier transform to return to the Laplace image, and by assuming ξ = 0 , one obtains the following:
W ˜ ( 0 , i q ) = 2 η P 1 i q K ( 0 , q ) 2 π η M 1 q 2 K ( 0 , q )
where
K ( ξ , q ) = e i p ξ d p D 1 ( p , q ) + D 2 ( p , q ) D 3 ( p , q ) D 2 ( p , q ) + D 3 ( p , q )
and finally
W ˜ ( ξ , q ) = 2 η P 1 i q K ( ξ , q ) 2 π η M 1 q 2 K ( 0 , q )
where without losing generality, W ˜ ( ξ , i q ) was switched to W ˜ ( ξ , q ) .
The final response in the time domain is obtained by the inverse Laplace transform, which in most cases can be well approximated by the sum of residues as a result of contour integration.
However, this paper is concerned with the identification of intervals of velocities where the mass movement is unstable, and more specifically, to the identification of the irregular cases. The case with μ s = 1 , κ p = 300 , η N = η s = 0 and variable damping levels is selected as an example of the regular, that is, the expected behavior. According to the previous section, this case has three resonances, which are calculated as α c r 1 = 0.707, α f c r = 1.899 and α c r 2 = 4.573. Without involving the D-decomposition method, instability lines can be directly determined by tracing the real-valued-induced frequency. This simply means to find a real q that fulfills the characteristic equation
2 π η M q 2 K ( 0 , q ) = 0
where subscript 1 is removed for simplicity and K ( 0 , q ) is proportional to the equivalent flexibility of the model. Because η M must be real, this implies finding a q for which K ( 0 , q ) is real, and thus, the full analysis can be kept within the real domain.
It can be easily verified that the lowest branch of the instability lines never intersects α c r 1 or α p c r ; therefore, instability always occurs in the supercritical velocity range. In the regular cases, instability lines are formed by two branches, as demonstrated in Figure 8. The first branch is delimited by α c r 1 (or α p c r ) and α c r 2 , and both ends tend to infinite η M . The second branch approaches α c r 2 from the left where it tends to infinite η M , and the other end tends to zero η M at infinite α . α f c r is ignored and has no influence on the instability lines.
With decreasing damping, the instability lines become closer to α c r j . From Figure 8 and Figure 9, it can be concluded that the left sides of both branches are already quite close to α c r j , and therefore, the decrease in damping does not have much influence. However, the right part of the first branch is highly affected by the damping level.
The irregular cases are identified by the fact that α c r 2 is intersected by one of the branches. Additionally, there can be a third branch, which in some cases can be closed. It was proven by parametric analysis within the ranges of allowable parameters that the irregular cases only occur for κ p on or below the line plotted in Figure 10. The curve in Figure 10 does not look smooth because the step in κ p was selected as 0.01 and in μ s , it was selected as 0.1. The line would be smoother for finer steps.
As an example, one case selected from the irregular region: μ s = 1.1 , κ p = 0.6 and η N = η s = 0 with variable damping levels is discussed. The instability lines are presented in Figure 11, where intersections with α c r 2 are shown. This case is also irregular from the point of view of the moving force; therefore, α p c r is plotted in Figure 11, demonstrating that the role of α c r 1 is fully secured by α p c r . For the sake of completeness, related q values are plotted in Figure 12.
The level of damping is indicated in the legend of Figure 11 and Figure 12. It is considered that both values are equal: η p = η f . For η p = η f = 1 10 6 and η p = η f = 1 10 4 , there are three instability branches from which one of them is closed. This is not very clear from Figure 11, because the closed branch is overlaid for these two levels of damping. The other two branches have similar behavior as in the regular cases. For η p = η f = 0.01 , the closed branch is broken, and there are only two instability branches. The first one is regular, and the second one crosses α c r 2 . For η p = η f = 0.05 , there are also only two branches, but both of them intersect α c r 2 and have no similarity with the regular behavior.

6. Instability of Two Moving Masses

When two masses are traversing the beam, then the only difference with respect to the previous analysis is the right-hand side of Equation (32), which reads as
4 ( 2 η P 1 η M 1 w ˜ , τ τ ) δ ( ξ ) + 4 ( 2 η P 2 η M 2 w ˜ , τ τ ) δ ( ξ d ˜ )
with additional parameters introduced as η P 2 = P 2 P 1 , η M 2 = M 2 χ m and d ˜ = d χ .
The solution steps can follow the same pattern as before, which allows to conclude that the characteristic equation is now
( π 2 η M 1 q 2 K ( 0 , q ) ) ( π 2 η M 2 q 2 K ( 0 , q ) ) 4 η M 1 η M 2 q 4 K ( d ˜ , q ) K ( d ˜ , q ) = 0
which for d ˜ = 0 reduces to the characteristic equation of η M = η M 1 + η M 2 , as expected.
To simplify the following analysis, it is assumed that η M 1 = η M 2 = η M and η P 1 = η P 2 = η P . Equation (46) is not as straightforward for tracing the instability lines as Equation (44) was. Nevertheless, Equation (46) is quadratic for η M ; therefore, it can be solved analytically for η M , and then two conditions for finding the real-valued induced frequency can be specified as
K ( d ˜ , q ) K ( d ˜ , q ) K ( 0 , q )   or   K ( d ˜ , q ) K ( d ˜ , q ) + K ( 0 , q )
Thus, if one of these functions has a real value for real q , then the corresponding η M can be computed from the associated root, and after that, the negative values can be discarded.
After detailed analysis, it can be concluded that when two moving masses are considered, then α c r j and α p c r are not respected. Any combination of allowable parameters creates situations where α c r 1 or α p c r are intersected, and thus, instability occurs in the subcritical velocity range. Additionally, α c r 2 does not serve to indicate some position of any of the asymptotes.
To identify the irregular cases, it is necessary to state when the instability occurs in the subcritical velocity range for some realistic η M . In railway applications, the mass and moving force are not linked by a mass-to-weight relation. Owing to the primary and the secondary suspension, it is a bit ambiguous what mass should be used in such simplified cases. It could be 10 t as a typical force transmitted by one wheel or just 880 kg, which is the wheel mass. Since according to Table 1, the value of χ can vary between 0.3 and 2.7 m−1, then the range for η M is 50–497 for 10 t and 4.5–44 for 880 kg. Therefore, η M = 100 was selected to be a limit value.
To avoid tracing the instability lines for all possible cases, the analysis is performed by identifying crossings with a specific q -value. In fact, only values as low as 0.1–0.15 are worth studying. Thus, μ s was stepped by 0.1, κ p was stepped by 2, real q was stepped by 0.01 and d ˜ was stepped by 0.25 starting at 1. The level of damping was selected as η p = η f = 0.05 . The tested region was specified by α spanning from 0 to α c r 1 or α p c r . The obtained results are presented graphically in Figure 13.
The graphs in Figure 13 are limited on the top by η M = 100 , as justified before, and on the left by α c r 1 or α p c r , which varies with κ p and cannot therefore be directly indicated in these graphs. In Figure 13a, there are only starting crossings, implying the beginning of a larger unstable region with a higher degree of instability, but for lower μ s , these values are quite close to α c r 1 . In Figure 13b,c, there are always two crossings, meaning that they delimitate an unstable interval, and after that, stability is recovered. This also means that in such cases, the rate of instability is generally low: lower than in the region with the starting branch. No irregular situation was found for d ˜ = 1 or higher μ s than indicated in Figure 13. For simplicity, higher values than κ p = 300 were not tested because the tendency for such high values is predictable. This is obvious from the fact that the points in graphs of Figure 13 are becoming closer as κ p increases by constant steps.
To conclude this analysis, the case with μ s = 1 and κ p = 300 , which was regular for one moving mass and now is irregular, is further examined for d ˜ = 1 ; 1.25 ; 1.5 ; 1.75 with η p = η f = 0.05 . All instability branches are plotted until α = 10 to cover also the region beyond α c r 2 . The results of this analysis are presented in Figure 14 and Figure 15. In Figure 14a, α is limited by 1.1, which is more realistic. However, to show how complicated instability branches can be, how the number of asymptotes is affected by d ˜ and the relation to α c r j and α p c r are not preserved; also, the region with higher α values is shown in Figure 14b,c for different scales of η M . These graphs should be accompanied by the associated real-valued q , which are reported separately for each d ˜ in Figure 15.
It is worthwhile to compare Figure 14 with Figure 13. If the limit η M = 100 is imposed, it can be concluded that for d ˜ = 1.5 and d ˜ = 1.75 , there is an unstable interval in the subcritical velocity range. Crossings with discrete q values reported in Figure 13 do not indicate the minimum η M , but this value can be extrapolated. For d ˜ = 1 , there is no value lower than η M = 100 in the subcritical velocity range, and for d ˜ = 1.25 , there is only the starting branch with values relatively close to α c r 1 . It should also be mentioned that by increasing the level of damping, the situation becomes generally worse [24,25].

7. Conclusions

The two-layer model of the ballasted railway track is a generally reliable model used by other researchers to analyze various aspects of its dynamic response induced by a moving railway vehicle. The model was experimentally validated 40 years ago. These results are published in [41]. There may be an issue in fitting the stiffness, damping and other parameters of the model, but once this is accomplished, the results obtained are realistic and trustworthy. Since wide variations in all parameters are considered in the manuscript, it is not possible to find a representative track for each combination to run the experiment and validate the results.
This paper presents a detailed analysis of the two-layer model; however, this is from a different point of view than in previous works. The main objective is to identify the range of parameters leading to some irregular behavior, where “irregular” means “different than expected”. This has a practical importance because unexpected behavior can complicate the railway track design. There are no published works that deal so widely with the identified irregularities and summarize the results over the wide range of admissible parameters, which identifies the novel contribution of this paper.
The irregular situations are defined differently for three distinct situations: one constant moving force and the associated critical velocities, and one or two moving masses and the associated velocity intervals of instability, identified by the instability lines. The sets of parameters corresponding to each of these three irregularities overlap but are not the same. To obtain a track with expected behavior, all of them must be avoided. All results are presented in a dimensionless form, mostly analytically or semi-analytically, but without any numerical integration or discretization.
The irregular case related to one moving force merely alerts that α c r 1 cannot be determined analytically, and therefore, parametric analysis must be used to identify α p c r . Without it, there would be a serious error in predicting the lowest velocity that should be avoided. Then, at α p c r , there is no real resonance, but an increase in displacements is verified and also the role of α p c r in instability issues is preserved. Therefore, this irregularity only requires a different approach in the determination of the relevant values without suggestions for the railway design.
Regarding the instability of one moving mass, it was concluded that the irregular cases occur throughout the allowable range of μ s but only for low values of κ p , corresponding to a strong foundation and soft rail pads. However, even this irregularity does not imply any recommendation for railway design because the irregularity affects the dynamic response at velocities around α c r 2 , which should be avoided in real scenarios. It merely means that α c r 2 is not playing the expected role. An important conclusion that was drawn is that the instability of one moving mass never occurs in the subcritical velocity range, considering the role of α p c r .
Therefore, the irregular cases that must be avoided in railway design are the ones that involve two moving proximate masses. However, it is not only important to know what values of η M and distances d ˜ lead to the instability in the subcritical velocity range; it is also essential to know how deep these values are located there, meaning what is the corresponding α . Another important fact is the rate of instability, which is generally lower in short, closed intervals of instability and higher at the beginning of a larger, possibly semi-infinite interval. The distance between the masses also plays a role. It can certainly be concluded that d ˜ = 1.5 and d ˜ = 1.75 may represent cases to be avoided when μ s is low, which corresponds to light sleepers such as wooden or similar. In addition, d ˜ = 1.25 is also to be avoided, even if α -values are not so low, but the rate of instability is expected to be higher. There is no indication that a lower value implies a worse case except perhaps for μ s . As for d ˜ , this is not true because d ˜ = 1 is not included in the irregular cases identified in Figure 13. Regarding κ p , this is also not true because, as indicated in the legend of Figure 13, the starting κ p value of several lines is quite high.
Definition of the dimensionless parameters and the presented graphs allow one to conclude how a certain result is affected by the model parameters as well as how to change the model parameters to increase or decrease this result value, respecting the actual needs. The easiest way to alter parameter κ p is by changing the stiffness of the rail pads. μ s can be altered by the sleepers’ selection; additionally, the reference value k f can be adjusted by ballast layer height or by adding or improving other reinforcing layers below it. Changes resulting from modifying the sleepers’ spacing are not very important, and introducing a different rail profile than that intended for such a railway is usually not possible.

Funding

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

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Frýba, L. Vibration of Solids and Structures under Moving Loads, Research Institute of Transport, Prague (1972), 3rd ed.; Thomas Telford: London, UK, 1999. [Google Scholar]
  2. Basu, D.; Kameswara Rao, N.S.V. Analytical solutions for Euler-Bernoulli beam on visco-elastic foundation subjected to moving load. Int. J. Numer. Anal. Methods Geomech. 2013, 37, 945–960. [Google Scholar] [CrossRef]
  3. Froio, D.; Rizzi, E.; Simões, F.M.F.; Da Costa, A.P. Universal analytical solution of the steady-state response of an infinite beam on a Pasternak elastic foundation under moving load. Int. J. Solids Struct. 2018, 132–133, 245–263. [Google Scholar] [CrossRef]
  4. Dieterman, H.A.; Metrikine, A.V. Steady state displacements of a beam on an elastic half-space due to uniformly moving constant load. Eur. J. Mech. A/Solids 1997, 16, 295–306. [Google Scholar]
  5. 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]
  6. 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]
  7. Kiani, K.; Nikkhoo, A. On the limitations of linear beams for the problems of moving mass-beam interaction using a meshfree method. Acta Mech. Sin. 2011, 28, 164–179. [Google Scholar] [CrossRef]
  8. Kiani, K.; Nikkhoo, A.; Mehri, B. Prediction capabilities of classical and shear deformable beam models excited by a moving mass. J. Sound Vib. 2009, 320, 632–648. [Google Scholar] [CrossRef]
  9. Kiani, K.; Nikkhoo, A.; Mehri, B. Assessing dynamic response of multispan viscoelastic thin beams under a moving mass via generalized moving least square method. Acta Mech. Sin. 2010, 26, 721–733. [Google Scholar] [CrossRef]
  10. Nelson, H.D.; Conover, R.A. Dynamic stability of a beam carrying moving masses. J. Appl. Mech. Trans. ASME 1971, 38, 1003–1006. [Google Scholar] [CrossRef]
  11. Benedetti, G.A. Dynamic stability of a beam loaded by a sequence of moving mass particles. J. Appl. Mech. Trans. ASME 1974, 41, 1069–1071. [Google Scholar] [CrossRef]
  12. Duffy, D.G. The response of an infinite railroad track to a moving, vibrating mass. J. Appl. Mech. 1990, 57, 66–73. [Google Scholar] [CrossRef]
  13. Metrikine, A.V.; Dieterman, H.A. Instability of vibrations of a mass moving uniformly along an axially compressed beam on a visco-elastic foundation. J. Sound Vib. 1997, 201, 567–576. [Google Scholar] [CrossRef]
  14. Metrikine, A.V.; Verichev, S.N. Instability of vibrations of a moving two-mass oscillator on a flexibly supported Timoshenko beam. Arch. Appl. Mech. 2001, 71, 613–624. [Google Scholar] [CrossRef]
  15. Yang, B.; Gao, H.; Liu, S. Vibrations of a Multi-Span Beam Structure Carrying Many Moving Oscillators. Int. J. Struct. Stab. Dyn. 2018, 18, 1850125. [Google Scholar] [CrossRef]
  16. Roy, S.; Chakraborty, G.; DasGupta, A. Coupled dynamics of a viscoelastically supported infinite string and a number of discrete mechanical systems moving with uniform speed. J. Sound Vib. 2018, 415, 184–209. [Google Scholar] [CrossRef]
  17. Nassef, A.S.E.; Nassar, M.M.; EL-Refaee, M.M. Dynamic response of Timoshenko beam resting on nonlinear Pasternak foundation carrying sprung masses. Iran. J. Sci. Technol. Trans. Mech. Eng. 2019, 43, 419–426. [Google Scholar] [CrossRef]
  18. Mazilu, T.; Dumitriu, M.; Tudorache, C. On the dynamics of interaction between a moving mass and an infinite one-dimensional elastic structure at the stability limit. J. Sound Vib. 2011, 330, 3729–3743. [CrossRef]
  19. Mazilu, T. Interaction between a moving two-mass oscillator and an infinite homogeneous structure: Green’s functions method. Arch. Appl. Mech. 2010, 80, 909–927. [Google Scholar] [CrossRef]
  20. Mazilu, T.; Dumitriu, M.; Tudorache, C. Instability of an oscillator moving along a Timoshenko beam on viscoelastic foundation. Nonlinear Dyn. 2012, 67, 1273–1293. [Google Scholar] [CrossRef]
  21. Mazilu, T. Instability of a train of oscillators moving along a beam on a viscoelastic foundation. J. Sound Vib. 2013, 332, 4597–4619. [Google Scholar] [CrossRef]
  22. Stojanović, V.; Kozić, P.; Petković, M.D. Dynamic instability and critical velocity of a mass moving uniformly along a stabilized infinity beam. Int. J. Solids Struct. 2017, 108, 164–174. [Google Scholar] [CrossRef]
  23. Stojanović, V.; Petković, M.D.; Deng, J. Stability and vibrations of an overcritical speed moving multiple discrete oscillators along an infinite continuous structure. Eur. J. Mech. A/Solids 2019, 75, 367–380. [Google Scholar] [CrossRef]
  24. Dimitrovová, Z. Dynamic interaction and instability of two moving proximate masses on a beam on a Pasternak viscoelastic foundation. Appl. Math. Model. 2021, 100, 192–217. [Google Scholar] [CrossRef]
  25. Dimitrovová, Z. Two-layer model of the railway track: Analysis of the critical velocity and instability of two moving proximate masses. Int. J. Mech. Sci. 2022, 217, 107042. [Google Scholar] [CrossRef]
  26. Stojanović, V.; Deng, J.; Petković, M.D.; Milić, D. Non-stability of a bogie moving along a specific infinite complex flexibly beam-layer structure. Eng. Struct. 2023, 295, 116788. [Google Scholar] [CrossRef]
  27. Stojanović, V.; Deng, J.; Milić, D.; Petković, M.D. Dynamics of moving coupled objects with stabilizers and unconventional couplings. J. Sound Vib. 2024, 570, 118020. [Google Scholar] [CrossRef]
  28. Dimitrovová, Z. Semi-analytical analysis of vibrations induced by a mass traversing a beam supported by a finite depth foundation with simplified shear resistance. Meccanica 2020, 55, 2353–2389. [Google Scholar] [CrossRef]
  29. Jahangiri, A.; Attari, N.K.A.; Nikkhoo, A.; Waezi, Z. Nonlinear dynamic response of an Euler–Bernoulli beam under a moving mass–spring with large oscillations. Arch. Appl. Mech. 2020, 90, 1135–1156. [Google Scholar] [CrossRef]
  30. Fărăgău, A.B.; Metrikine, A.V.; van Dalen, K.N. Transition radiation in a piecewise-linear and infinite one-dimensional structure–a Laplace transform method. Nonlinear Dyn. 2019, 98, 2435–2461. [Google Scholar] [CrossRef]
  31. Fărăgău, A.B.; Keijdener, C.; de Oliveira Barbosa, J.M.; Metrikine, A.V.; van Dalen, K.N. Transition radiation in a nonlinear and infinite one-dimensional structure: A comparison of solution methods. Nonlinear Dyn. 2021, 103, 1365–1391. [Google Scholar] [CrossRef]
  32. Fărăgău, A.B.; Mazilu, T.; Metrikine, A.V.; Lu, T.; van Dalen, K.N. Transition radiation in an infinite one-dimensional structure interacting with a moving oscillator—The Green’s function method. J. Sound Vib. 2021, 492, 115804. [Google Scholar] [CrossRef]
  33. Koh, C.G.; Ong, J.S.Y.; Chua, D.K.H.; Feng, J. Moving element method for train-track dynamics. Int. J. Numer. Methods Eng. 2003, 56, 1549–1567. [Google Scholar] [CrossRef]
  34. Ang, K.K.; Dai, J. Response analysis of high-speed rail system accounting for abrupt change of foundation stiffness. J. Sound Vib. 2013, 332, 2954–2970. [Google Scholar] [CrossRef]
  35. Tran, M.T.; Ang, K.K.; Luong, V.H. Vertical dynamic response of non-uniform motion of high-speed rails. J. Sound Vib. 2014, 333, 5427–5442. [Google Scholar] [CrossRef]
  36. Dai, J.; Lim, J.G.Y.; Ang, K.K. Dynamic response analysis of high-speed maglev-guideway system. J. Vib. Eng. Technol. 2023. [Google Scholar] [CrossRef]
  37. Xu, Y.; Zhu, W.; Fan, W.; Yang, C.; Zhang, W. A new three-dimensional moving Timoshenko beam element for moving load problem analysis. J. Vib. Acoust. Trans. ASME 2020, 142, 031001. [Google Scholar] [CrossRef]
  38. Elhuni, H.; Basu, D. Novel Nonlinear Dynamic Beam-Foundation Interaction Model. ASCE J. Eng. Mech. 2021, 147, 04021012. [Google Scholar] [CrossRef]
  39. 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]
  40. Sadri, M.; Lu, T.; Steenbergen, M. Railway track degradation: The contribution of a spatially variant support stiffness-Global variation. J. Sound Vib. 2020, 464, 114992. [Google Scholar] [CrossRef]
  41. 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. [Google Scholar] [CrossRef]
  42. Dimitrovová, Z. On the Critical Velocity of Moving Force and Instability of Moving Mass in Layered Railway Track Models by Semianalytical Approaches. Vibration 2023, 6, 113–146. [Google Scholar] [CrossRef]
  43. Yang, C.J.; Xu, Y.; Zhu, W.D.; Fan, W.; Zhang, W.H.; Mei, G.M. A three-dimensional modal theory-based Timoshenko finite length beam model for train-track dynamic analysis. J. Sound Vib. 2020, 479, 115363. [Google Scholar] [CrossRef]
  44. Ramos, A.; Castanheira-Pinto, A.; Colaço, A.; Fernández-Ruiz, J.; Costa, P.A. Predicting Critical Speed of Railway Tracks Using Artificial Intelligence Algorithms. Vibration 2023, 6, 895–916. [Google Scholar] [CrossRef]
  45. Grassie, S.L.; Gregory, R.W.; Harrison, D.; Johnson, K.L. The dynamic response of railway track to high frequency vertical excitation. J. Mech. Eng. Sci. 1982, 24, 77–90. [Google Scholar] [CrossRef]
  46. Rodrigues, A.S.F. 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]
  47. Van Dalen, K.N. Ground Vibration Induced by a High-Speed Train Running over Inhomogeneous Subsoil, Transition Radiation in Two-Dimensional Inhomogeneous Elastic Systems. Master’s Thesis, Department of Structural Engineering, TUDelft, Delft, The Netherlands, 2006. [Google Scholar]
  48. Chen, Y.-H.; Huang, Y.-H.; Shih, C.-T. Response of an Infinite Timoshenko Beam on a Viscoelastic Foundation to a Harmonic Moving Load. J. Sound Vib. 2001, 241, 809–824. [Google Scholar] [CrossRef]
  49. 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]
  50. Sañudo, R.; Miranda, M.; Alonso, B.; Markine, V. Sleepers Spacing Analysis in Railway Track Infrastructure. Infrastructures 2022, 7, 83. [Google Scholar] [CrossRef]
  51. 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]
  52. Ortega, R.S.; Pombo, J.; Ricci, S.; Miranda, M. The importance of sleepers spacing in railways. Constr. Build. Mater. 2021, 300, 124326. [Google Scholar] [CrossRef]
  53. Dimitrovová, Z.; Rodrigues, A.S.F. Critical Velocity of a Uniformly Moving Load. Adv. Eng. Softw. 2012, 50, 44–56. [Google Scholar] [CrossRef]
Figure 1. The two-layer model of a ballasted railway track subjected to an axial force and traversing by two proximate masses acted upon by constant vertical forces.
Figure 1. The two-layer model of a ballasted railway track subjected to an axial force and traversing by two proximate masses acted upon by constant vertical forces.
Applsci 13 12356 g001
Figure 2. Separation of regions with one and three resonances in κ p - μ s plane for the simplified case with η N = η s = 0 .
Figure 2. Separation of regions with one and three resonances in κ p - μ s plane for the simplified case with η N = η s = 0 .
Applsci 13 12356 g002
Figure 3. Results of parametric analyses: (a) three resonances; (b) one resonance; (c,d) three resonances but one of them is double. In the legend: “max” maximum displacement over the entire beam; “min” minimum displacement over the entire beam; “AP” deflection at the active point.
Figure 3. Results of parametric analyses: (a) three resonances; (b) one resonance; (c,d) three resonances but one of them is double. In the legend: “max” maximum displacement over the entire beam; “min” minimum displacement over the entire beam; “AP” deflection at the active point.
Applsci 13 12356 g003aApplsci 13 12356 g003b
Figure 4. Graphs of α c r 1 as a function of κ p in logarithmic scale for discrete values of μ s from 1 to 5.5 by 0.5 (ten curves) and η N = η s = 0 : (a) κ p > 100 (Formula (30) is not plotted because it perfectly matches the analytical value; no irregular cases are detected); (b) κ p 0.03 ; 100 (Formula (30) is plotted in green but cannot be used as already mentioned; the region of irregular cases in agreement with Figure 2 is clearly identified).
Figure 4. Graphs of α c r 1 as a function of κ p in logarithmic scale for discrete values of μ s from 1 to 5.5 by 0.5 (ten curves) and η N = η s = 0 : (a) κ p > 100 (Formula (30) is not plotted because it perfectly matches the analytical value; no irregular cases are detected); (b) κ p 0.03 ; 100 (Formula (30) is plotted in green but cannot be used as already mentioned; the region of irregular cases in agreement with Figure 2 is clearly identified).
Applsci 13 12356 g004
Figure 5. Graphs of all resonances and the associated Fourier variable as a function of κ p in logarithmic scale, for discrete values of μ s from 1 to 5.5 by 0.5 (ten curves) and η N = η s = 0 : (a) α c r 1 (orange), α f c r (blue) and α c r 2 (yellow) for κ p 0.03 ; 100 (region of irregular cases is clearly identified); (b) Fourier variables associated to the previous values in the same color scheme indicating the reason for missing resonances.
Figure 5. Graphs of all resonances and the associated Fourier variable as a function of κ p in logarithmic scale, for discrete values of μ s from 1 to 5.5 by 0.5 (ten curves) and η N = η s = 0 : (a) α c r 1 (orange), α f c r (blue) and α c r 2 (yellow) for κ p 0.03 ; 100 (region of irregular cases is clearly identified); (b) Fourier variables associated to the previous values in the same color scheme indicating the reason for missing resonances.
Applsci 13 12356 g005
Figure 6. Graphs of resonances as a function of κ p in logarithmic scale for discrete values of μ s from 1 to 5.5 by 0.5 (ten curves) and η N = 0 , η s = 0.2 : (a) α c r 1 for κ p > 100 (Formula (30) is not plotted because it perfectly matches the semi-analytical value; no irregular cases are detected); (b) α c r 1 (orange), α f c r (blue) and α c r 2 (yellow) for κ p 0.03 ; 100 (region of irregular cases is clearly identified).
Figure 6. Graphs of resonances as a function of κ p in logarithmic scale for discrete values of μ s from 1 to 5.5 by 0.5 (ten curves) and η N = 0 , η s = 0.2 : (a) α c r 1 for κ p > 100 (Formula (30) is not plotted because it perfectly matches the semi-analytical value; no irregular cases are detected); (b) α c r 1 (orange), α f c r (blue) and α c r 2 (yellow) for κ p 0.03 ; 100 (region of irregular cases is clearly identified).
Applsci 13 12356 g006
Figure 7. Graphs of resonances as a function of κ p in logarithmic scale for discrete values of μ s from 1 to 5.5 by 0.5 (ten curves) and η N = 0.2 , η s = 0 : (a) α c r 1 for κ p > 100 (Formula (30) is not plotted because it perfectly matches the semi-analytical value; no irregular cases are detected); (b) α c r 1 (orange), α f c r (blue) and α c r 2 (yellow) for κ p 0.03 ; 100 (region of irregular cases is clearly identified).
Figure 7. Graphs of resonances as a function of κ p in logarithmic scale for discrete values of μ s from 1 to 5.5 by 0.5 (ten curves) and η N = 0.2 , η s = 0 : (a) α c r 1 for κ p > 100 (Formula (30) is not plotted because it perfectly matches the semi-analytical value; no irregular cases are detected); (b) α c r 1 (orange), α f c r (blue) and α c r 2 (yellow) for κ p 0.03 ; 100 (region of irregular cases is clearly identified).
Applsci 13 12356 g007
Figure 8. Instability lines for a regular case with μ s = 1 , κ p = 300 , η N = η s = 0 and two levels of damping as indicated in the legend. (CV1 = α c r 1 , FCV = α f c r , CV2 = α c r 2 . Both CV1 and CV2 are represented by the same dashed line because they have the primary effect on the instability lines. This introduces no ambiguity because, by definition, CV1 < CV2.).
Figure 8. Instability lines for a regular case with μ s = 1 , κ p = 300 , η N = η s = 0 and two levels of damping as indicated in the legend. (CV1 = α c r 1 , FCV = α f c r , CV2 = α c r 2 . Both CV1 and CV2 are represented by the same dashed line because they have the primary effect on the instability lines. This introduces no ambiguity because, by definition, CV1 < CV2.).
Applsci 13 12356 g008
Figure 9. Real-valued mass-induced frequency linked to the instability lines from Figure 8. (CV1 = α c r 1 , FCV = α f c r , CV2 = α c r 2 . Both CV1 and CV2 are represented by the same dashed line because they have the primary effect on the instability lines. This introduces no ambiguity because, by definition, CV1 < CV2).
Figure 9. Real-valued mass-induced frequency linked to the instability lines from Figure 8. (CV1 = α c r 1 , FCV = α f c r , CV2 = α c r 2 . Both CV1 and CV2 are represented by the same dashed line because they have the primary effect on the instability lines. This introduces no ambiguity because, by definition, CV1 < CV2).
Applsci 13 12356 g009
Figure 10. Separation of regular and irregular regions in μ s - κ p plane for the simplified case with η N = η s = 0 .
Figure 10. Separation of regular and irregular regions in μ s - κ p plane for the simplified case with η N = η s = 0 .
Applsci 13 12356 g010
Figure 11. Instability lines for μ s = 1.1 , κ p = 0.6 , η N = η s = 0 and four levels of damping as indicated in the legend (PCV = α p c r , CV = α c r 2 ).
Figure 11. Instability lines for μ s = 1.1 , κ p = 0.6 , η N = η s = 0 and four levels of damping as indicated in the legend (PCV = α p c r , CV = α c r 2 ).
Applsci 13 12356 g011
Figure 12. Real-valued mass induced frequency linked to the instability lines from Figure 11 (PCV = α p c r , CV = α c r 2 ).
Figure 12. Real-valued mass induced frequency linked to the instability lines from Figure 11 (PCV = α p c r , CV = α c r 2 ).
Applsci 13 12356 g012
Figure 13. Identification of irregular cases for two moving masses: (a) d ˜ = 1.25 ; (b) d ˜ = 1.5 ; (c) d ˜ = 1.75 . The numbers in the legend indicate μ s , starting κ p and q , respectively. The arrow indicates the direction of κ p increase. For better clarity, single values are indicated by markers. An additional letter “S” indicates that this branch is starting without terminating the instability region.
Figure 13. Identification of irregular cases for two moving masses: (a) d ˜ = 1.25 ; (b) d ˜ = 1.5 ; (c) d ˜ = 1.75 . The numbers in the legend indicate μ s , starting κ p and q , respectively. The arrow indicates the direction of κ p increase. For better clarity, single values are indicated by markers. An additional letter “S” indicates that this branch is starting without terminating the instability region.
Applsci 13 12356 g013aApplsci 13 12356 g013b
Figure 14. Instability lines for case characterized by μ s = 1 , κ p = 300 , η N = η s = 0 and η p = η f = 0.05 : (ac) parts correspond to different scales of α and η M , the dotted line 0.05 identifies the behavior of a single moving mass (CV2 = α c r 2 ). (In (b), CV1 = α c r 1 , FCV = α f c r , CV2 = α c r 2 . Both CV1 and CV2 are represented by the same dashed line because they have the primary effect on the instability lines. This introduces no ambiguity because, by definition, CV1 < CV2.).
Figure 14. Instability lines for case characterized by μ s = 1 , κ p = 300 , η N = η s = 0 and η p = η f = 0.05 : (ac) parts correspond to different scales of α and η M , the dotted line 0.05 identifies the behavior of a single moving mass (CV2 = α c r 2 ). (In (b), CV1 = α c r 1 , FCV = α f c r , CV2 = α c r 2 . Both CV1 and CV2 are represented by the same dashed line because they have the primary effect on the instability lines. This introduces no ambiguity because, by definition, CV1 < CV2.).
Applsci 13 12356 g014aApplsci 13 12356 g014b
Figure 15. Real-valued mass-induced frequency linked to the instability lines from Figure 14: (a) d ˜ = 1 ; (b) d ˜ = 1.25 ; (c) d ˜ = 1.5 ; (d) d ˜ = 1.75 . (CV1 = α c r 1 , FCV = α f c r , CV2 = α c r 2 . Both CV1 and CV2 are represented by the same dashed line because they have the primary effect on the instability lines. This introduces no ambiguity because, by definition, CV1 < CV2.).
Figure 15. Real-valued mass-induced frequency linked to the instability lines from Figure 14: (a) d ˜ = 1 ; (b) d ˜ = 1.25 ; (c) d ˜ = 1.5 ; (d) d ˜ = 1.75 . (CV1 = α c r 1 , FCV = α f c r , CV2 = α c r 2 . Both CV1 and CV2 are represented by the same dashed line because they have the primary effect on the instability lines. This introduces no ambiguity because, by definition, CV1 < CV2.).
Applsci 13 12356 g015aApplsci 13 12356 g015bApplsci 13 12356 g015c
Table 1. Range of typical values of the main components of the two-layer model.
Table 1. Range of typical values of the main components of the two-layer model.
ParameterApproximate Range
E I (MNm2)4.9–6.4
m (kg/m) distributed54–60
M s (kg) concentrated40–160 [46]
k p (MN/m) concentrated20–5000 [46]
k f (MN/m2) distributed0.22–1000 [47,48]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Dimitrovová, Z. Instability of Vibrations of Mass(es) Moving Uniformly on a Two-Layer Track Model: Parameters Leading to Irregular Cases and Associated Implications for Railway Design. Appl. Sci. 2023, 13, 12356. https://doi.org/10.3390/app132212356

AMA Style

Dimitrovová Z. Instability of Vibrations of Mass(es) Moving Uniformly on a Two-Layer Track Model: Parameters Leading to Irregular Cases and Associated Implications for Railway Design. Applied Sciences. 2023; 13(22):12356. https://doi.org/10.3390/app132212356

Chicago/Turabian Style

Dimitrovová, Zuzana. 2023. "Instability of Vibrations of Mass(es) Moving Uniformly on a Two-Layer Track Model: Parameters Leading to Irregular Cases and Associated Implications for Railway Design" Applied Sciences 13, no. 22: 12356. https://doi.org/10.3390/app132212356

APA Style

Dimitrovová, Z. (2023). Instability of Vibrations of Mass(es) Moving Uniformly on a Two-Layer Track Model: Parameters Leading to Irregular Cases and Associated Implications for Railway Design. Applied Sciences, 13(22), 12356. https://doi.org/10.3390/app132212356

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop