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,
and
, point masses modeling the sleeper half-masses,
, and supported further by linear spring–damper elements representing the vertical stiffness and damping of all other layers,
and
, 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
. 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
where
is the spatial coordinate and
is the time. In addition to the parameters defined in the first paragraph of this section,
and
are the bending stiffness and mass per unit length of the beam,
and
are the unknown vertical displacement fields of the beam and discrete masses, respectively, and
means the load. Derivatives are denoted by the respective variables in subscripts preceded by a comma. The load is given by
where
is the Dirac delta function,
and
are the moving forces,
and
are the moving masses,
is the constant velocity and
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
, which will result in continuous supports. In
Figure 1, there is also a shear component
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
with unit N/m is used, then
with the correct unit of N.
In Equation (3),
and
are the masses displacements. Owing to the assumption of the tight contact,
holds with
and
. Therefore, the relevant derivatives using the chain rule are
Subsequently, the loading term reads as follows:
Nevertheless, the solution will be presented in moving coordinates
,
, leading to the following system of equations (compare with [
24,
25,
42]):
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
and
. 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
The other values are related to the reference Winkler beam by
where for simplicity, the previous designations are kept, but all values (except for
) are assumed as longitudinally distributed.
Therefore, the decisive parameters are
and
. Assuming that the sleepers are spaced from
= 0.545 m [
49,
50] to
= 0.711 m [
50], this together with values in
Table 1 identifies a relatively narrow interval for
, but it gives a very large variability for
. A typical spacing in Portugal [
51], in most European and some other countries, e.g., India, is 0.6 m. According to [
50],
should be related to the sleeper type and ballast height. There are studies analyzing how significant savings could be achieved by increasing
[
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
where, besides the parameters specified in Equations (9)–(11),
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
where
and
,
and
represent
,
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
-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
. The determinant can be written as a cubic polynomial with real coefficients for
, and the second equation as a quadratic polynomial with real coefficients for
In order to express directly
, Equation (21), multiplied by
, should be subtracted from Equation (20) multiplied by 2. Then
which is simplified when
to
Thus, finding the real double -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 , , and ) can be solved for valid and consequently and .
The irregular case was defined as the case with only one resonance. This analysis is presented for the simpler case of
. After substituting Equation (23) into Equation (21), a cubic equation is obtained for
. Its discriminant is simplified to
The irregular case occurs when expression (24) is negative, which means that there is only one real root for
. Expression (24) is a fourth-order polynomial for
and third-order polynomial for
. To identify the separation between three and one resonance,
is solved analytically as a function of
.The analytical values of
are given by
where
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
until
and then by
for larger
. It has an asymptotic tendency to
, but this only occurs for
outside the allowable range. It can be observed that for
, there is always an open interval of
for which there is only one resonance. For
, there is only a double root
leading to a collapsed interval, and therefore, one of the three resonances is doubled. The same happens when the selection
and
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,
and for better clarity, the scale does not start at
. As an example for the regular situation,
,
are selected, and for the irregular one:
,
. 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 (
), the false critical velocity (
) and the higher critical velocity (
), as predicted. In the irregular case, as shown in
Figure 3b, it is observed that the only analytically determined value of 0.996 is
and
is replaced by the pseudocritical velocity (
) 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
and
, which corresponds to the case where the interval of irregular cases is collapsed to a single value, proving that
is suppressed by
(both equal to 0.707) and the case is irregular with a unique
of 1.125. The same happens in the case plotted in
Figure 3d for
and
, which corresponds to the lower end of the irregular case interval, where
is again suppressed by
(both equal to 0.553) and the case is irregular with a unique
of 0.658.
This analysis is consistent with correspondence with finite beams according to [
53]. As explained in [
42],
corresponds to local minima in the resonant velocity plot, while
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,
must be determined by parametric analysis. This is important for two reasons. First,
may not provide the relevant information about the lowest velocity where excessive vibrations occur. Second, when it comes to instability,
plays the role of
and therefore must be known in advance to draw correct conclusions regarding the instability of moving inertial objects. In addition,
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
is ambiguous.
A very good estimate of
is provided by extension of the classical formula by considering the involved stiffnesses in series and by summing the masses, i.e., by
which is valid for higher
, approximately for
, which, in fact, means that for such values,
is practically independent on
and thus
can also be used. This identifies the influential parameters on
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
and
.
The results presented in
Figure 4,
Figure 5,
Figure 6 and
Figure 7 allow one to conclude how
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
is by changing the stiffness of the rail pads,
can be altered by the sleepers’ selection; additionally, the reference value
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
where
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
to keep its formalism,
will be used in the following with
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
and also
which shows that whenever there was
in Equations (17)–(19) is now
.
In Equation (36), there is still
, which must be solved for. For this, first, Equation (36) is solved for
.
Now, one can perform the inverse Fourier transform to return to the Laplace image, and by assuming
, one obtains the following:
where
and finally
where without losing generality,
was switched to
.
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
,
,
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
= 0.707,
= 1.899 and
= 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
that fulfills the characteristic equation
where subscript 1 is removed for simplicity and
is proportional to the equivalent flexibility of the model. Because
must be real, this implies finding a
for which
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
or
; 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
(or
) and
, and both ends tend to infinite
. The second branch approaches
from the left where it tends to infinite
, and the other end tends to zero
at infinite
.
is ignored and has no influence on the instability lines.
With decreasing damping, the instability lines become closer to
. From
Figure 8 and
Figure 9, it can be concluded that the left sides of both branches are already quite close to
, 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
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
on or below the line plotted in
Figure 10. The curve in
Figure 10 does not look smooth because the step in
was selected as 0.01 and in
, it was selected as 0.1. The line would be smoother for finer steps.
As an example, one case selected from the irregular region:
,
and
with variable damping levels is discussed. The instability lines are presented in
Figure 11, where intersections with
are shown. This case is also irregular from the point of view of the moving force; therefore,
is plotted in
Figure 11, demonstrating that the role of
is fully secured by
. For the sake of completeness, related
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:
. For
and
, 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
, the closed branch is broken, and there are only two instability branches. The first one is regular, and the second one crosses
. For
, there are also only two branches, but both of them intersect
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
with additional parameters introduced as
,
and
.
The solution steps can follow the same pattern as before, which allows to conclude that the characteristic equation is now
which for
reduces to the characteristic equation of
, as expected.
To simplify the following analysis, it is assumed that
and
. Equation (46) is not as straightforward for tracing the instability lines as Equation (44) was. Nevertheless, Equation (46) is quadratic for
; therefore, it can be solved analytically for
, and then two conditions for finding the real-valued induced frequency can be specified as
Thus, if one of these functions has a real value for real , then the corresponding 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 and are not respected. Any combination of allowable parameters creates situations where or are intersected, and thus, instability occurs in the subcritical velocity range. Additionally, 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
. 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
is 50–497 for 10 t and 4.5–44 for 880 kg. Therefore,
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
-value. In fact, only values as low as 0.1–0.15 are worth studying. Thus,
was stepped by 0.1,
was stepped by 2, real
was stepped by 0.01 and
was stepped by 0.25 starting at 1. The level of damping was selected as
. The tested region was specified by
spanning from 0 to
or
. The obtained results are presented graphically in
Figure 13.
The graphs in
Figure 13 are limited on the top by
, as justified before, and on the left by
or
, which varies with
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
, these values are quite close to
. 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
or higher
than indicated in
Figure 13. For simplicity, higher values than
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
increases by constant steps.
To conclude this analysis, the case with
and
, which was regular for one moving mass and now is irregular, is further examined for
with
. All instability branches are plotted until
to cover also the region beyond
. 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
and the relation to
and
are not preserved; also, the region with higher
values is shown in
Figure 14b,c for different scales of
. These graphs should be accompanied by the associated real-valued
, which are reported separately for each
in
Figure 15.
It is worthwhile to compare
Figure 14 with
Figure 13. If the limit
is imposed, it can be concluded that for
and
, there is an unstable interval in the subcritical velocity range. Crossings with discrete
values reported in
Figure 13 do not indicate the minimum
, but this value can be extrapolated. For
, there is no value lower than
in the subcritical velocity range, and for
, there is only the starting branch with values relatively close to
. 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 cannot be determined analytically, and therefore, parametric analysis must be used to identify . Without it, there would be a serious error in predicting the lowest velocity that should be avoided. Then, at , there is no real resonance, but an increase in displacements is verified and also the role of 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 but only for low values of , 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 , which should be avoided in real scenarios. It merely means that 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 .
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
and distances
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
and
may represent cases to be avoided when
is low, which corresponds to light sleepers such as wooden or similar. In addition,
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
. As for
, this is not true because
is not included in the irregular cases identified in
Figure 13. Regarding
, this is also not true because, as indicated in the legend of
Figure 13, the starting
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 is by changing the stiffness of the rail pads. can be altered by the sleepers’ selection; additionally, the reference value 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.