Semi-Analytical Approach and Green’s Function Method: A Comparison in the Analysis of the Interaction of a Moving Mass on an Infinite Beam on a Three-Layer Viscoelastic Foundation at the Stability Limit—The Effect of Damping of Foundation Materials

In this paper, the interaction of a mass moving uniformly on an infinite beam on a three-layer viscoelastic foundation is analyzed with the objective of determining the lowest velocity at the stability limit, called, in this context, the critical velocity. This issue is important for rail transport and, in particular, for the high-speed train, because the moving mass is the basic model of a vehicle, and the infinite beam on a three-layer viscoelastic foundation is the usual mechanical representation of the railway track. In addition to this, the advantages and disadvantages of the two implemented methods, namely, the semi-analytical approach and the Green’s function method, are summarized in terms of computational time, the precision of the obtained results, limitations, and the feasibility of implementation. All results are presented in a dimensionless form to cover a wide range of possible scenarios. Some results may be considered academic, however, results related to a particular railway track are also included. Particular attention is paid to the influence of the damping of materials in the foundation upon the critical velocity of the moving mass. Regarding the semi-analytical approach, it is demonstrated that the critical velocities can be obtained in an exact manner by tracing the branches of the so-called instability lines in the velocity–moving-mass plane. This analysis can be maintained within the real domain. As for the time series, they can be determined by a numerical inverse Laplace transform. Moreover, thanks to the analytical form of the final result in the Fourier domain, each value corresponding to a specific time instant can be obtained directly, that is, without the previous time history. Regarding the Green’s function method, this is used to verify a few points delimiting the stable and unstable regions of the moving mass with the help of the D-decomposition approach. Additionally, a numerical algorithm based on the Green’s function and convolution integral written for dimensionless quantities is used to calculate the time series of the moving mass. In addition to identifying the critical velocity of the moving mass, its connection with the critical velocity of the moving force is emphasized, and the possibility of validating the results on long finite beams using modal expansion is presented and described.


Introduction
The railway has the function of providing for the circulation of railway vehicles in a safe, economical, and comfortable way.For all of these aspects to be guaranteed, it is necessary that there is an adequate transfer of efforts to the ground, taking into account the high demands that the rolling stock imposes on the rail; therefore, one of the most challenging problems in the development of railway passenger transport is the identification of the most suitable technical solutions capable of ensuring an increase in the traffic speed under standardized conditions with regard to the comfort of passengers and the ride quality of the vehicles.From this point of view, the quality of the track plays a determining role in terms of the intensity of mechanical stresses to which the entire vehicle/track system is exposed during operation.Indeed, many reports have indicated a sudden increase in track deformation and acceleration at the level of the rolling gear of the vehicle associated with exceeding a certain speed threshold when trains run on excessively flexible tracks [1,2].The starting point for overcoming such obstacles is a better knowledge of the nature of the interaction phenomena between the vehicle and the track based on theoretical and experimental approaches.
Mechanical representations of the track must consider its structure and its static and dynamic properties associated with the frequency domain of interest.In the following discussion, the structure of the classic track with concrete sleepers is considered.Figure 1 shows the main components of a classic track.The two rails of the track are placed on concrete sleepers, separated by rail pads, and clamped to the sleepers by the rail fastening system, which includes the sleeper bolt and a tension clamp.The frame, consisting of the two rails and the concrete sleepers, is placed on the crushed stone prism (ballast).The elasticity and damping of the track are mainly ensured by the materials from which the rail pads, ballast, and subgrade (the layer under the ballast) are made.

Introduction
The railway has the function of providing for the circulation of railway vehicles in a safe, economical, and comfortable way.For all of these aspects to be guaranteed, it is necessary that there is an adequate transfer of efforts to the ground, taking into account the high demands that the rolling stock imposes on the rail; therefore, one of the most challenging problems in the development of railway passenger transport is the identification of the most suitable technical solutions capable of ensuring an increase in the traffic speed under standardized conditions with regard to the comfort of passengers and the ride quality of the vehicles.From this point of view, the quality of the track plays a determining role in terms of the intensity of mechanical stresses to which the entire vehicle/track system is exposed during operation.Indeed, many reports have indicated a sudden increase in track deformation and acceleration at the level of the rolling gear of the vehicle associated with exceeding a certain speed threshold when trains run on excessively flexible tracks [1,2].The starting point for overcoming such obstacles is a better knowledge of the nature of the interaction phenomena between the vehicle and the track based on theoretical and experimental approaches.
Mechanical representations of the track must consider its structure and its static and dynamic properties associated with the frequency domain of interest.In the following discussion, the structure of the classic track with concrete sleepers is considered.Figure 1 shows the main components of a classic track.The two rails of the track are placed on concrete sleepers, separated by rail pads, and clamped to the sleepers by the rail fastening system, which includes the sleeper bolt and a tension clamp.The frame, consisting of the two rails and the concrete sleepers, is placed on the crushed stone prism (ballast).The elasticity and damping of the track are mainly ensured by the materials from which the rail pads, ballast, and subgrade (the layer under the ballast) are made.In the range of very low frequencies (0.5-20 Hz), that is, in the range in which the natural frequencies of the vehicle carbody and its suspensions are located, the track model can be composed of bodies with concentrated parameters [3][4][5] or even be rigid [6,7].This is because, on the one hand, the receptance of the vehicle train at the wheels is much higher than the track receptance in this frequency range, and due to this, the wheels adopt the displacement caused by the track's irregularities, while the rail remains practically fixed.On the other hand, the bending waves of the rail are evanescent, so the coupling between vehicle wheels due to this type of wave can be neglected.
At low and medium frequencies (6-700 Hz), the track vibrates like a two-DOF oscillator, exhibiting two natural frequencies and an antiresonance frequency between them.In addition to evanescent waves, there are also propagating waves that arise between the low natural frequency and the antiresonance frequency and at frequencies higher than the high natural frequency [8].The track model is based on continuous medium theory, in which the rail is modeled as an Euler-Bernoulli or Timoshenko beam, the sleepers are assimilated by a purely inertial layer, and the viscoelastic properties of the rail pads and ballast prism are modeled using the Winkler foundation with damping, which comprises the so-called two-layer model.Additional shear effects can be introduced considering a In the range of very low frequencies (0.5-20 Hz), that is, in the range in which the natural frequencies of the vehicle carbody and its suspensions are located, the track model can be composed of bodies with concentrated parameters [3][4][5] or even be rigid [6,7].This is because, on the one hand, the receptance of the vehicle train at the wheels is much higher than the track receptance in this frequency range, and due to this, the wheels adopt the displacement caused by the track's irregularities, while the rail remains practically fixed.On the other hand, the bending waves of the rail are evanescent, so the coupling between vehicle wheels due to this type of wave can be neglected.
At low and medium frequencies (6-700 Hz), the track vibrates like a two-DOF oscillator, exhibiting two natural frequencies and an antiresonance frequency between them.In addition to evanescent waves, there are also propagating waves that arise between the low natural frequency and the antiresonance frequency and at frequencies higher than the high natural frequency [8].The track model is based on continuous medium theory, in which the rail is modeled as an Euler-Bernoulli or Timoshenko beam, the sleepers are assimilated by a purely inertial layer, and the viscoelastic properties of the rail pads and ballast prism are modeled using the Winkler foundation with damping, which comprises the so-called two-layer model.Additional shear effects can be introduced considering a two-parameter (i.e., Pasternak) foundation [9,10].However, it would be more realistic to include dynamically activated mass in the foundation as well [11].A characteristic feature of this frequency range is the fact that the rail response does not change significantly along the sleeper spacing, which justifies neglecting the influence of equidistantly spaced sleepers on the rail receptance.
When the frequency increases, the response of the rail changes dramatically due to the influence of the bending modes of the rail on the sleepers.Among them, the most influential is the first eigenmode, whose wavelength is equal to twice the distance between the sleepers.The rail receptance shows a maximum at the mid-span between the sleepers at the pinned-pinned resonance frequency and a minimum (antiresonance) above the sleepers but at a slightly higher frequency.To simulate this rail behavior, the track model should consist of a beam resting on discrete supports [12][13][14].
Several discretization methods for the continuous elastic medium can be used to improve the accuracy of the prediction of the dynamic regime of the track.Thus, for problems specific to high frequencies, the finite element method is used because it allows a more accurate mechanical description of the track components [15][16][17].Particularly complex behavior of the components in a ballast bed can be studied using the discrete element method [18,19].To solve problems related to running noise, the boundary element method is used [20,21].In this context, the moving element method [22,23] and the new moving beam element [24] or the novel model of nonlinear dynamic interaction [25] should also be mentioned.
The track representation is useful for investigating many of the issues that the vehicle/track interaction raises for railway engineers and researchers.An interesting topic from a technical and academic point of view is the phenomenon of the instability of the vehicle/track system, which can be observed in certain cases when the velocity of the vehicle exceeds the velocity of propagation of the elastic waves generated in the track.This occurs when the track infrastructure is built on poorly consistent soils [1,2,26] and/or the rails are subject to excessive expansion [27,28].Indeed, when the velocity of a moving object on a continuous elastic structure exceeds the velocity of propagation of the elastic waves induced in the structure, the moving object radiates anomalous Doppler waves that produce a negative damping effect at the object/structure interface, resulting in the destabilization of the whole system under certain circumstances [29,30].And this phenomenon could occur in the case of a railway vehicle if the above conditions are met.
The problem of the stability of a moving object on an elastic structure with applications to the vehicle/track system was treated in the past by considering different mechanical representations of the vehicle (moving object) and the track (elastic structure).A railway vehicle can be represented by simple models, such as a moving mass [28, [31][32][33], two moving masses [34,35], a two-mass oscillator [36][37][38] or a train of two-mass oscillators [39][40][41], a three-mass oscillator [42], or a three-body oscillator [43,44].The set of moving inertial objects was also dealt with in pioneering works [45,46] and, more recently, in [47].A more complex vehicle model can also be used to study the stability issue [48,49].The track is modeled using an Euler-Bernoulli infinite beam on a continuous foundation with one elastic layer [28,29,31,32,34] or two [35,50,51] or three elastic layers [52].Instead of the Euler-Bernoulli beam, the Timoshenko beam can also be used [42,43,53].In order to take into account the influence of the sleeper bay, a model of an infinite beam on a foundation with discrete supports should be considered for the track [54].
Other works concern the limitations of classical beam theories [55][56][57].Large deformations are investigated in [58], initial conditions are accounted for in [59], and the non-uniform motion of the inertial object is analyzed in [60].Related problems include the variation in foundation stiffness [61,62] or analyses concerned with tracks other than classical ballasted tracks [63].Related problems are also of interest, for instance, complex issues in which a moving object interacts with an elastic structure, as happens in a pantograph-catenary system [64,65].
The main objective of the analyses presented in this paper is to identify the critical velocity of the moving object, which, in this context, will indicate the lowest velocity at which the moving object/elastic structure system becomes unstable.In this regard, it is necessary to emphasize that the critical velocity of the moving force, i.e., of an object without inertia, is different and not related to the anomalous Doppler effect.Nevertheless, the critical velocity of the moving force is closely related to the aim of this paper.This topic has been deeply analyzed over the course of several years: pioneering works are summarized in the monograph in [66], more complicated models are approached analytically in [67], and more recently, artificial intelligence has also been used for this complicated problem in real situations [68].
When investigating the quality of a system to be stable according to the velocity of a moving object, in certain cases, depending on the nature of the moving object and the elastic structure, zones of stability alternating with zones of instability have been identified [42,51,52], but the first zone on the velocity axis is always a zone of stability, and therefore, the critical velocity of a moving object refers to the transition from this zone to the first zone of instability that follows it.
The determination of stability/instability zones and critical velocity is related to the nature of the roots of the characteristic equation associated with the linear equations of motion of the moving object/elastic structure system.All roots of the characteristic equation have a negative real part in the stability zone, and at least one solution of the characteristic equation has a positive real part in the instability zone.The difficulty of the problem lies in the fact that the characteristic equation does not have a polynomial form, and therefore, the solution to the problem can be obtained by a graph-analytical method, the so-called D-decomposition method [69].This method maps the zones where the characteristic equation has stable or unstable roots for a given value of the moving object velocity in the complex plane associated with the parameter of interest of the moving object/elastic structure system by sweeping the angular frequency axis [26,28,29,36,[42][43][44]48,50,70].
Another objective of the problem of the stability of a moving object on an elastic structure refers to the description of the dynamic behavior when the system loses its stability to assess the potential consequences.Achieving this goal requires taking into account the loss of contact between the moving object and the structure.In this case, the equations of motion become nonlinear, and solving them in the time domain can be performed using a method based on the Green's function [32,40].In this way, a nonlinear stability analysis based on the Hopf bifurcation is possible [42].
Relatively recently, a new semi-analytical method for analyzing the response of an elastic structure to a moving object has been proposed [10,11,34,35,38,51,52,59]. The main objective of this method is to determine the full dynamic response of the structure directly as a sum of residues.In this way, the vibrations of the structure can be monitored from the initial instant with respect to inhomogeneous initial conditions [59] and without loss of accuracy.Very good results are achieved with one-layer and two-layer models [10,34,35,38].An indication of the stability of the moving object is then obtained as a secondary result with respect to the main objective mentioned above.For a more complex supporting structure, the sum of residues must be supplemented by additional numerical integration [11] due to discontinuities in the equivalent flexibility of the supporting structure.However, these discontinuities never occur in unstable regions; therefore, stability or instability can also be identified directly from the characteristic equation, as in the D-decomposition method.In the solution method, the variable related to the frequency is modified so that unstable cases are identified by the existence of at least one root with a negative imaginary part.The identification of instability starts with the selection of the velocity of the moving object, and for such a value, a real-valued frequency providing a real-valued equivalent flexibility is found using the classic intersection method, and then the corresponding moving mass is calculated as a real-valued result.Each time the equivalent flexibility is required, it is precisely calculated by contour integration, which ensures the high accuracy of the obtained results.A smooth change in velocity will cause a smooth change in the value of the moving mass.As a result, the so-called instability lines, which can have several branches, are obtained.After scanning the whole velocity-moving-mass plane, by choosing a particular value for the moving mass, velocity intervals of stability and instability can be obtained from the intersection with instability branches.These intervals also indicate the number of unstable components contained in the full solution [52].
With regard to the previous discussion, it can be concluded that the numerical evaluation of the dynamic behavior of railways or, in general, structures subjected to moving loads is undergoing great development, both with the use of numerical methods and with the help of other approaches, such as semi-analytical methods or methods implementing the dynamic Green's function.
This paper is focused on an infinite beam supported by three viscoelastic layers, which, due to its computational efficiency and relatively good approximation to reality, is commonly used by railway engineers for the evaluation of several relevant aspects.The new findings that are presented concern the instability of a moving mass.The critical velocity in this context will be used for the lowest velocity that separates stable and unstable behavior.However, the new contribution of this paper lies not only in the presentation of relevant results but also in the detailed comparison of the semi-analytical and dynamic Green's function approach from the viewpoint of several aspects, such as computational efficiency, the accuracy of the obtained results, ease of implementation, limitations, etc.
It is demonstrated that there are several situations with rather unpredictable behavior, especially at low damping levels, as already shown for simpler, one-layer [34] and two-layer models [35,51].In the semi-analytical approach, the so-called instability lines are traced, indicating the moving mass as a function of velocity in a limiting situation where one of the induced frequencies is real-valued.It is shown that these lines have several branches.They can have quite a strange development, as they can form a closed curve.Open curves must end with asymptotes tending to infinite mass at a fixed velocity, except for at least one end of an instability branch, which tends to zero mass at infinite velocity.
It is also shown that the expected relation to the critical velocity of the moving force is not confirmed in all cases.First, three values of such a critical velocity do not always exist.Sometimes, there is only one, and then the others have to be compensated by pseudocritical velocities that may or may not be dominant, which, in turn, affects their influence on the instability lines.Second, especially for low damping levels, there are more asymptotes than these velocities predict.Third, there are cases where the instability lines intersect the lines of critical velocities of the moving force, which contradicts the regular behavior of the Winkler-Pasternak beam.
This paper is organized as follows: First, the model, along with its basic assumptions and simplifications, is presented in Section 2, where the governing equations are also given.The solutions obtained by the two analyzed methods, that is, the semi-analytical approach and the Green's function method, are derived in Sections 3 and 4, respectively.The set of dimensionless parameters necessary to describe the model are identified in Section 3.They are also used in Section 4, where some additional parameters are also introduced.Section 5 summarizes the range of allowable real railway track model parameters alongside the range of admissible dimensionless parameters.Parameters identifying a specific track section from the literature are also given.Several cases are selected to demonstrate the results.They refer to the specific track section as well as to other track sections identified by dimensionless parameters within the allowable intervals, including the damping of materials in the foundation.Finally, the two solution methods are compared.The main conclusions are drawn in Section 6.

Mechanical Model and Governing Equations
As already mentioned in the Introduction, the three-viscoelastic-layer railway track model is widely used for its simplicity, numerical efficiency, and ability to effectively estimate the real behavior, as demonstrated in [71,72].Its description was already given in previous works on layered models [51,52], but it is repeated here for completeness.In this model, the rail is modeled by a beam that is supported by linear spring-damper elements standing for rail pads, which are connected to point masses representing sleepers' masses, which are, in turn, supported by linear spring-damper elements modeling the ballast, including its shear resistance.Further point masses represent the dynamically activated ballast mass, and these are supported by other linear spring-damper components modeling the foundation.This model is adapted for the analysis of transversal vibrations; therefore, the rail pad and foundation springs act only in the vertical direction, and the rotational degree of freedom is neglected.The shear resistance of the model is therefore introduced exclusively by the ballast shear stiffness elements.
It would be more correct to introduce discrete supports according to the sleeper spacing; nevertheless, as indicated in the Introduction, discrete supports are only important at high frequencies, more precisely around the pinned-pinned value.Therefore, continuous supports are introduced in this paper because they are more suitable for semi-analytical derivations, as well as for derivations exploiting the dynamic Green's function.Discrete supports are thus replaced by continuous ones.Such a simplification still provides a model that serves as a very good approximation of reality.It is also worthwhile to remark that small deflection variations between individual supports can equally be modeled by a harmonic component added to the constant moving force.
A constant mass is assumed to move over the model at a constant velocity, as depicted in Figure 2. The aim of the analysis is to determine the lowest velocity that marks the separation between stable and unstable behavior-the so-called critical velocity.The assumptions for such an analysis are listed below [51]:  The equations of motion can be written as follows [52]: where p(x, t) is the load, and w(x, t), us(x, t), and ub(x, t) are the unknown deflections of the beam, the sleepers' masses, and the ballast's masses, respectively.EI and m stand for the bending stiffness and mass per unit length of the beam, and N is the axial force.kp, cp, kb, cb, and kf, cf are stiffness and viscous damping parameters related to the rail pads, ballast, and foundation, respectively.ks models the ballast shear stiffness, and ms and mb are the sleepers' and ballast's point masses.x is the fixed spatial coordinate, and t is the time.In the semi-analytical approach, fixed (tight) contact between the mass and the beam is assumed, which means that the mass displacement and the displacement of the respective beam axis are always the same.An extension to the situation with a linear contact spring can be easily implemented in conformity with previous works [38], but it is not presented in this paper for simplicity.The Green's function method, on the other hand, has better numerical stability when a contact spring is implemented, as will be demonstrated in Section 5.
The equations of motion can be written as follows [52]: where p(x, t) is the load, and w(x, t), u s (x, t), and u b (x, t) are the unknown deflections of the beam, the sleepers' masses, and the ballast's masses, respectively.EI and m stand for the bending stiffness and mass per unit length of the beam, and N is the axial force.k p , c p , k b , c b , and k f , c f are stiffness and viscous damping parameters related to the rail pads, ballast, and foundation, respectively.k s models the ballast shear stiffness, and m s and m b are the sleepers' and ballast's point masses.x is the fixed spatial coordinate, and t is the time.Derivatives are designated by the variable in the subscript position, preceded by a comma.In this formulation, all parameters are assumed to be distributed.To achieve this, the discrete values are divided by the sleepers' spacing l s , except for the shear stiffness, which is modeled by a vertical spring, and therefore, a uniform distribution has to be applied to the inverse value.Therefore, if a discrete value of k s with the unit N/m is used, then the distributed value is k s = k s l s with the unit N.
As already mentioned, a constant mass M and a constant force P move with the same uniform velocity v.Then, the load p(x, t) is given by where δ is the Dirac delta function, and w 0 (t) is the displacement of the mass.The initial conditions are considered homogeneous, as already specified in assumption (v), and the boundary conditions state that the deflection and its slope vanish for x → ±∞ .

Solution of the Governing Equations
In order to obtain the deflection fields, firstly, it is convenient to switch the fixed spatial coordinate to a moving one, defined by r = x − vt; then [52], Further, dimensionless parameters are introduced because all results will be presented in a dimensionless form.Actual track-specific data are also used, but, in general, all parameters are considered within the range of allowable values that will be defined in Section 5.
A Winkler beam with bending stiffness EI, mass per unit length m, and foundation stiffness k f is chosen as a reference.Then: where v ref is the critical velocity of the force moving over the reference beam, and α stands for the velocity ratio.The dimensionless spatial coordinate and time follow as Displacement fields are related to the maximum static deflection induced by force P on the reference beam; thus, The mass and stiffness ratios are the damping ratios read and the axial force and shear ratios are Finally, the moving mass and force ratios are Defining the moving force is certainly abundant, but it is mathematically more correct and lends itself to further generalizations.
Using the definitions above, one obtains the following: ) For the stability analysis, the Laplace transform must be applied first to capture the initial instant because the transient part of the vibrations is essential: In fact, q = iq, but Equation ( 19) is used to keep the formalism of the Laplace transform.For homogeneous initial conditions, one obtains Then, the Fourier transform is applied as follows: Without losing generality, functional dependence can be adapted to q instead of q, leading to The expressions introduced by Equations ( 24)-( 26) can be simplified by introducing Then, Equations ( 24)-( 26) can be written as (the dependence on the variables is omitted for simplicity, except for the term to be solved) In order to remove W(0, q), firstly, Equation ( 33) is solved for W(p, q): where d 3 is the determinant of the system (33), and is the subdeterminant of order 2, obtained as a determinant of the matrix of Equation ( 33) with the first line and column cut.Now, one can perform the inverse Fourier transform to recover the Laplace image: Then, ξ = 0 can be substituted, and by introducing one obtains Substituting Equation (38) into Equation ( 36), one finally obtains It is worth noting that these derivations are conceptually the same as in previous works [10,52]; the only difference compared to other layered models is the definition of the function K(ξ, q), which is proportional to the equivalent flexibility of the supporting structure.
The final step lies in the inverse Laplace transform, which is defined by Likewise, as in previous works, it is found that it is more convenient to switch the real and imaginary axes; thus, the final deflection shape of the beam can be numerically evaluated as follows: where a > 0 must be chosen so that all discontinuities and poles of the function to be integrated lie above the line (−ia − ∞, ia + ∞) in the complex q-plane.
With the help of contour integration methods, the following can be derived: where ⌢ w tr (ξ, τ) is the truly transient part, and res designates the residue.Then, where the first term stands for the steady-state solution, and q M j in the second term represents mass-induced frequencies as solutions of the characteristic equation.Methods for the determination of mass-induced frequencies and their properties are detailed in previous works [10,38,59].K(ξ, 0), K ξ, q M j , and K ,q 0, q M j can be evaluated accurately by contour integration and analytical terms for residues of simple and double poles.Therefore, contour integration is used in two stages.One instance is for K and K ,q calculations in the p-complex plane for a given q.In this case, the p-poles can be determined numerically, and the residues can be obtained analytically.The other usage is more complicated, because it relates to evaluation in the q-complex plane to derive Equation (43).To achieve this, the q-poles of the function to be integrated in Equation ( 41) must be determined.It is clear from Equation ( 41) that one pole is obvious, that is, 0, and defines the steady-state part of the solution.The other poles, q M j , are solutions of the characteristic equation π − 2η M q 2 K(0, q) = 0.They are referred to as mass-induced frequencies, and the associated residues identify the harmonic part of the transient vibration.The last part of Equation ( 43), w tr (ξ, τ), also belongs to the transient part of the solution and is obtained numerically.The reason for its existence is that the K-function has step discontinuities of bounded values.For contour integration, branch cuts have to be introduced to remove such locations from the complex q-plane.Integration along branch cuts leads to w tr (ξ, τ).While w tr (ξ, τ) is generally less important for one-or two-layer models, unfortunately, this is not the case for three-layer models.There are many regions where q M j does not exist, and then the sum of residues is insufficient to approximate the full solution.Therefore, a numerical calculation according to Equation (41) must usually be performed.

Critical Velocity of a Moving Mass
The critical velocity of the moving mass is defined here as the lowest α that marks the separation between stable and unstable regions.An important finding in the semi-analytical approach is that the discontinuity lines as a function of q are always located at positions with the q of positive imaginary parts, and Equation ( 43) implies that instability occurs when the imaginary part of q M j is negative.Hence, for the stability analysis, which is the goal of this paper, it is not important to analyze w tr (ξ, τ) or to determine the discontinuity lines.Instability lines are identified by real-valued poles q M j , since, in the damped case, assuming all problem parameters are fixed and taking q M j as a function of α, this necessarily implies switching between q M j with positive and negative imaginary parts.These values are not affected by the discontinuity lines, and thus, their tracing as a function of α is not interrupted.They can be determined in a very simple way.Taking the characteristic equation π − 2η M q 2 K(0, q) = 0, it is clear that it suffices to test K(0, q) as a function of α and find a real q for which K(0, q) is also real.This can be realized by simple iteration on the real q-axis.All such positions can be determined without any doubt that some value was missed.Calculations can be performed with very high accuracy using the exact value of K(0, q) obtained semi-analytically by contour integration.The corresponding η M is then calculated.Further, all branches are plotted on the α-η M plane.
Instability lines therefore mark the separation between q M j with positive and negative imaginary parts, meaning between stable and unstable behavior, or between two unstable regions with different degrees of instability.These cases can be easily distinguished from the graphs of branches of instability lines.Then, for any selected η M , it is possible to identify α-intervals where the system is unstable.This is a simple task compared to searching for q M j in general, which implies solving the complex roots of a complex equation.

Critical Velocity of a Moving Force
In order to analyze the results indicated in the previous section, it is important to know the critical velocities of the structure, i.e., the critical velocities that can be assigned to one moving constant force, because there is a tight relationship between the critical velocities of a moving force and a moving mass.
The critical velocity of the moving force can be determined by analyzing the steadystate deflection.When the force is passing over the structure at a critical velocity, then resonance occurs in the undamped case, and the deflection tends to infinity.
For the steady-state analysis, the time dependence can be removed.In addition, the loading term is reduced to one moving constant force, and thus, the previous Equations ( 16)-( 18) simplify to Then, in the Fourier space, by using Equation ( 23), one obtains the following: With the introduction of which is formally the same as in Equation (33), only all terms q − αp are simplified to −αp, and the loading term includes only the moving force: Thus, where d st 3 and d st 2 have the same meaning as for the moving mass, only the terms that are involved are adapted to the steady-state situation, as indicated by the superscript.Contour integration can be used to accurately calculate deflection shapes, avoiding the numerical problems inherent to the FFT, including the necessity of introducing damping for numerical stability.
To identify the critical velocity, it is necessary to find a double real p-root of d st 3 in the undamped case and the corresponding α st cr , that is, to find the real pair p, α st cr > 0 that solves d st 3 = 0 and d st 3,p = 0 when η p = η b = η f = 0.It is easy to verify that for η p = η b = η f = 0, d st 3 is a cubic polynomial for α = α 2 and a fourth-order polynomial for p = p 2 with real coefficients.α is real and positive, and so is α.p is real and p is therefore real and positive because zero values are excluded.d st 3,p can be divided by 2p, and then the resulting expression is still a cubic polynomial for α = α 2 and a cubic polynomial for p = p 2 .Simplifications like in [35] are not possible; however, it is possible to numerically solve the exact value of α st cr .It is also possible to show that there can only be 1, 3, or 5 valid pairs p, α st cr .To find the roots in the simpler case where η N = η s = 0, it is possible to consider pd st 3,p /(2p) − 3d st 3 , which is a quadratic polynomial for α.Its two solutions can be introduced into d st 3 or d st 3,p /(2p).In both cases, all p-roots can be obtained directly with Maple's "solve" function.If necessary, the number of digits can be increased, but no additional numerical methods are needed.Real roots are always paired with their opposite values, and complex roots exist with all four combinations of signs.Only positive real p stays for a valid solution, and only this value yields a positive α.Since complex roots are always in groups of four, there can be 2, 6, or 10 real roots, except in cases with multiple roots.It is then found that there are 1, 3, or 5 resonant velocities.

Long Finite Beam-Eigenmode Expansion
It is also advantageous to perform the analysis on long finite beams because all results presented for infinite beams can be validated on long finite beams, as demonstrated in previous works [10,34,35,52].It was proven that simple supports can be introduced, but the load should start to act a little farther from the support.In the context of this paper, it is important to analyze the critical velocities of the moving force, because this analysis allows identifying the so-called false critical velocity.
To carry out such an analysis, the vibration modes and orthogonality conditions must be determined.In order to keep the analysis in the real domain, undamped modes are considered.Free undamped vibrations in fixed coordinates fulfill the following: Admitting harmonic vibrations w(ξ, τ) = w(ξ)e iϖτ , u s (ξ, τ) = u s (ξ)e iϖτ , and u b (ξ, τ) = u b (ξ)e iϖτ , one obtains the following: Then, assuming w(ξ) = We ipξ , u s (ξ) = Ve ipξ , and u b (ξ) = Ue ipξ yields the following: which can be written in a matrix form: The nullity of the determinant can be used to solve for the natural frequencies.For the given data and admitting that L is the beam length, p = jπ/(Lχ) can be used owing to the simple supports.For the analysis presented here, it is better to keep p.The determinant is a cubic equation for ϖ 2 : Therefore, Equation ( 68) has an analytical solution that gives three simple real values for ϖ 2 when the discriminant is positive, which is always true due to the physical nature of such an equation.Owing to its length, the analytical solution is not included here.
The resonance condition expresses the equality of excitation and natural frequencies: Then, the critical velocity should correspond to α at a local minimum.The stationary condition d dp can be solved for p, but only real and positive solutions are valid solutions.It can be verified that each frequency indicates several stationary values, but at most, five of them identify valid p and α.Specifically, there are only one, three, or five valid pairs (α, p) at resonance.When ordered according to α, then the odd values identify local minima and thus the critical velocity, but the even values stand for local maxima, identifying the false critical velocity.The false critical velocity still induces a resonance in the sense that the steady-state solution tends to infinity, but the usual properties of the critical velocity are not preserved, and furthermore, these values have no influence on the instability lines.

Green's Function Method
In this section, the above problem is solved using the D-decomposition method by applying the Green's function in the frequency domain.The verification of the obtained results is carried out by evaluating the time response of the system using a numerical method for calculating the convolution integral [70].To avoid numerical instability caused by the fact that the mass/beam contact is perfectly rigid, a linear elastic element with sufficiently high stiffness is introduced so that it does not significantly affect the result regarding the critical velocity.

Stability Issue
Equations of motion ( 16)-( 18) describe both the steady-state behavior caused by a static load and the perturbed one.Stability depends on a small free vibration around the equilibrium point, and therefore, the homogeneous equations of motion must be extracted from above in matrix form: where L ξ,τ stands for the matrix differential operator with where and The quantities in Equation ( 76) are the system perturbations around the steady-state position.
Im Z(α, ±iq)) becomes the border between the stable and unstable parts of the complex plane, and the curve (Re Z(α, ζ ± iq), Im Z(α, ζ ± iq)) with ζ > 0 passes through an unstable part.When the point (1, 0) belongs to the stable part (see Figure 3a or Figure 4a), the system is stable.Figures 3b and 4b show the position of the point (1, 0) when the system is unstable.Unstable part of the complex plane can be to the left of the curve (Re Z(α, ±iq), Im Z(α, ±iq)), as in Figure 3, or to the right, as in Figure 4.

Time-Domain Response
In this section, an algorithm for calculating the time series of a moving mass on an infinite beam on a three-layer viscoelastic foundation based on Green's function theory and the convolution theorem is presented in terms of dimensionless quantities.
To avoid numerical instability when the contact between the moving mass and the beam is rigid, an elastic element of stiffness kc is inserted between the moving mass and the beam.
The nonhomogeneous term in Equation (1) becomes where Pc(t) is the contact force between the moving mass and the beam.In addition, the equation of motion of the moving mass, where z(t) is the vertical mass displacement, and the equation of the moving mass/beam contact where H(.) is the Heaviside step function, should be further considered.By introducing the Heaviside step function, the loss of contact between the moving mass and the beam is simulated.
Changing the reference frame, the nonhomogeneous term in Equation ( 6) becomes Pc(t)δ(r), and the equation of the contact reads For rewriting in dimensionless quantities, three new dimensionless parameters must be introduced, in conformity with Equations ( 11), (12), and (15):

Time-Domain Response
In this section, an algorithm for calculating the time series of a moving mass on an infinite beam on a three-layer viscoelastic foundation based on Green's function theory and the convolution theorem is presented in terms of dimensionless quantities.
To avoid numerical instability when the contact between the moving mass and the beam is rigid, an elastic element of stiffness k c is inserted between the moving mass and the beam.
The nonhomogeneous term in Equation (1) becomes where P c (t) is the contact force between the moving mass and the beam.In addition, the equation of motion of the moving mass, where z(t) is the vertical mass displacement, and the equation of the moving mass/ beam contact where H(.) is the Heaviside step function, should be further considered.By introducing the Heaviside step function, the loss of contact between the moving mass and the beam is simulated.
Changing the reference frame, the nonhomogeneous term in Equation ( 6) becomes P c (t)δ(r), and the equation of the contact reads For rewriting in dimensionless quantities, three new dimensionless parameters must be introduced, in conformity with Equations ( 11), (12), and (15): The interaction of a moving mass on an infinite beam on a three-layer viscoelastic foundation is then described by Equations ( 16)- (18), where the right-hand side of Equation ( 16) is adjusted to and by Equations ( 101) and ( 103), which, in the dimensionless form, are given by The solutions of the equations of motion of the beam ( 16)-( 18) with the modification (105) and of the moving mass (106) with (107) can be specified using the time-domain Green functions associated with the beam and the moving mass by where g w (ξ, ε, τ) is the time-domain Green function of the beam with respect to the moving reference frame, and g z (τ) is the time-domain Green function of the moving mass.The time-domain Green function of the beam in the moving reference frame is obtained by applying the inverse Fourier transform to the frequency-domain Green function of the beam: where the frequency-domain Green function of the beam takes the following form: However, for ξ = ε = 0, Equation ( 110) is reduced to where g w (τ) = g w (0, 0, τ).The above integral can be evaluated using a numerical method at any time instant [32,42].
The time-domain Green function of the moving mass is Finally, the interaction of the moving mass on the infinite beam on a three-layer viscoelastic foundation is described by Equations ( 108)-( 109) and ( 106), which can be solved numerically.To this end, the time-domain Green function of the beam, g w (τ nς), and the dimensionless contact force, η P c (ς), are assumed to be linear within the time-step intervals τ i-1 < ς < τ i : where Dτ = τ i − τ i−1 .For τ = τ n , it reads where By introducing Equations ( 115) and (116) into Equation (117), a linear equation for η P c ,n is obtained, and the dimensionless contact force at time τ n is given as its solution.The dimensionless displacement of the beam and of the moving mass can then be calculated using Equations ( 115) and (116).

Allowable Intervals of Dimensionless Parameters
In this section, the allowable ranges of dimensionless parameters are identified.To accomplish this, one must first know the reasonable ranges of possible real values.As for the rail, the range of possible values is quite narrow; basically, there are only two guide data sets for 54 E1 and 60 E1, specifying values of EI and m.The sleepers can be wooden or concrete, and therefore, their mass varies from approximately 80 kg to 320 kg [71].Since only one rail is considered, only half the sleeper mass is introduced.Admitting sleeper spacings from 0.545m [73] to 0.711 m [74], the distributed masses for half sleepers can be determined.There are studies on larger sleeper spacings [75] that could provide significant cost savings; nevertheless, these values are not yet common in practice and will therefore not be considered.As for the limit values for k f , they are taken from [76,77], respectively.
As for the ballast contribution, the stress cone (prism) theory can be used to estimate some realistic values [72,73].Admitting the effective sleeper length according to [71] and then varying the sleepers' spacing from 0.545 to 0.711 m, the sleepers' base from 0.22 to 0.35 m, the ballast height from 0.2 to 0.6 m, and the angle of distribution from 25 • to 50 • , the formulas from [72] indicate the stiffness coefficient range 0.847-3.261m and the volume of dynamically activated mass coefficient range 0.069-0.630m 3 .Additionally, the shear contribution coefficient range is 0.03-0.5 m.Admitting the ballast Young's modulus from 50 MPa to 400 MPa, the ballast density within 1200-2600 kg/m 3 , and the Poisson ratio within 0.1-0.45,real values can be determined and then distributed according to the sleeper spacing.All of these values are summarized in Table 1.The admissible ranges of the dimensionless parameters are summarized in Table 2.These ranges are extended by a certain margin, since all considered real values are not exact limits.η s is sometimes omitted in similar models so that the lower bound is allowed to be zero and only the upper bound is estimated.Regarding the damping coefficients, very different values can be found in the literature, so it was decided to directly vary the dimensionless parameters.
Table 2. Range of some dimensionless values of the main parts of the three-layer model.

Dimensionless Parameter
Approximate Range (with Margins) As for η M , academic values will be used here to explain in detail all the tendencies affecting the onset of instability.In general, there are studies that associate mass with force as its weight, and then the mass could be as high as 10 t, but if only the wheel mass is considered, then it can be as low as 880 kg.Additionally, there are also intermediate approaches that add a corresponding axle mass or part of the boogie mass to the wheel mass.According to Table 1, the value of χ varies between 0.3 and 2.7 m −1 , giving the range of η M as 50-500 for 10 t and 4.4-44 for 880 kg.
In addition to studies with parameters selected from the allowable ranges, one specific case of a real railway is also considered.This case is taken from [73].Using the real values from [73], the dimensionless parameters are calculated as

Test Case from [73]
In this section, a real case of the railway track from [73] is selected for analysis.First, the critical velocity of the moving force is determined because it provides important data for further study.According to Section 3.3, it can be shown that this case has only one resonance.A valid pair for this is p, α st cr = (1.456, 0.983), which can also be confirmed by the analysis described in Section 3.4.Parametric analysis reveals that there are two non-dominant pseudocritical velocities (PCV1 and PCV2) lower than the critical one (CV).Nevertheless, by gradually increasing κ p from 0.8 up to 8, which can be easily achieved by increasing the stiffness of the rail pads, still in a range related to soft rubber, it can be seen that while CV is increasing, the PCVs are gaining its importance.The results of these parametric analyses are summarized in Figure 5.
increasing the stiffness of the rail pads, still in a range related to soft rubber, it can be seen that while CV is increasing, the PCVs are gaining its importance.The results of these parametric analyses are summarized in Figure 5.Further investigation shows how these values affect the instability lines.Even if the CVs and PCVs are different for different j, there are basically three distinct regions where the instability branches are contained.For the original value κp = 0.8, the PCVs are not dominant, and therefore, there is only one instability branch, as shown in Figure 6.For the highest value tested, κp = 8, the branches of the instability lines are contained in all three regions and do not cross the CV or PCVs, as demonstrated in Figure 7.For completeness, the corresponding real-valued qM values are also plotted in Figure 8.
A detailed analysis shows that there is a branch in the first region already from j = 2 and in the second region from j = 3.This is important from an analytical point of view; from a practical point of view, it is also necessary to evaluate whether the ηM values in these branches are feasible.Further investigation shows how these values affect the instability lines.Even if the CVs and PCVs are different for different j, there are basically three distinct regions where the instability branches are contained.For the original value κ p = 0.8, the PCVs are not dominant, and therefore, there is only one instability branch, as shown in Figure 6.increasing the stiffness of the rail pads, still in a range related to soft rubber, it can be seen that while CV is increasing, the PCVs are gaining its importance.The results of these parametric analyses are summarized in Figure 5.Further investigation shows how these values affect the instability lines.Even if the CVs and PCVs are different for different j, there are basically three distinct regions where the instability branches are contained.For the original value κp = 0.8, the PCVs are not dominant, and therefore, there is only one instability branch, as shown in Figure 6.For the highest value tested, κp = 8, the branches of the instability lines are contained in all three regions and do not cross the CV or PCVs, as demonstrated in Figure 7.For completeness, the corresponding real-valued qM values are also plotted in Figure 8.
A detailed analysis shows that there is a branch in the first region already from j = 2 and in the second region from j = 3.This is important from an analytical point of view; from a practical point of view, it is also necessary to evaluate whether the ηM values in these branches are feasible.For the highest value tested, κ p = 8, the branches of the instability lines are contained in all three regions and do not cross the CV or PCVs, as demonstrated in Figure 7.For completeness, the corresponding real-valued q M values are also plotted in Figure 8.

Other Test Cases
In this section, results related to several other cases are presented.It is concluded that they can be classified as regular, i.e., with expected behavior, and irregular.The regular cases behave as in Section 5.2, meaning that the following apply: There is at most one instability branch in each of the three regions delimited by CVs and PCVs; (ii) No branch intersects CV or PCV; (iii) Branches that correspond to lower damping are below the ones with higher damping, and they do not cross; (iv) In the first two regions, the branches asymptotically tend to infinite ηM, and in the last region, they asymptotically tend to infinite ηM from the left and zero ηM from the right.
This regularity is not associated with the number of resonances.Irregular cases then violate one of these properties, especially at low damping levels.

Other Test Cases
In this section, results related to several other cases are presented.It is concluded that they can be classified as regular, i.e., with expected behavior, and irregular.The regular cases behave as in Section 5.2, meaning that the following apply: There is at most one instability branch in each of the three regions delimited by CVs and PCVs; (ii) No branch intersects CV or PCV; (iii) Branches that correspond to lower damping are below the ones with higher damping, and they do not cross; (iv) In the first two regions, the branches asymptotically tend to infinite ηM, and in the last region, they asymptotically tend to infinite ηM from the left and zero ηM from the right.
This regularity is not associated with the number of resonances.Irregular cases then violate one of these properties, especially at low damping levels.A detailed analysis shows that there is a branch in the first region already from j = 2 and in the second region from j = 3.This is important from an analytical point of from a practical point of view, it is also necessary to evaluate whether the η M values in these branches are feasible.

Other Test Cases
In this section, results related to several other cases are presented.It is concluded that they can be classified as regular, i.e., with expected behavior, and irregular.The regular cases behave as in Section 5.2, meaning that the following apply: (i) There is at most one instability branch in each of the three regions delimited by CVs and PCVs; (ii) No branch intersects CV or PCV; (iii) Branches that correspond to lower damping are below the ones with higher damping, and they do not cross; (iv) In the first two regions, the branches asymptotically tend to infinite η M , and in the last region, they asymptotically tend to infinite η M from the left and zero η M from the right.
This regularity is not associated with the number of resonances.Irregular cases then violate one of these properties, especially at low damping levels.
Cases selected for further analysis are listed in Table 3, in which "nd" and "d" designate dominant and non-dominant, respectively.Cases selected for further analysis are listed in Table 3, in which "nd" and "d" designate dominant and non-dominant, respectively.In Case 1 (Figure 9), all five resonances are clearly marked.The first two are quite close; however, the sudden jump to zero deflection at the active point for velocities between CV1 and PCV is undoubtedly noticeable.In Case 2 (Figure 10), PCV1 is dominant, but PCV2 is not.Nevertheless, both are ambiguous because, for an α-step of 0.001, the extreme values on the entire beam are achieved at the same α reported in Table 3, but the extreme values at the active point occur at 0.148 and 0.581, respectively.In Case 1 (Figure 9), all five resonances are clearly marked.The first two are quite close; however, the sudden jump to zero deflection at the active point for velocities between CV1 and PCV is undoubtedly noticeable.
In Case 2 (Figure 10), PCV1 is dominant, but PCV2 is not.Nevertheless, both are ambiguous because, for an α-step of 0.001, the extreme values on the entire beam are achieved at the same α reported in Table 3, but the extreme values at the active point occur at 0.148 and 0.581, respectively.In Case 3 (Figure 11), PCV1 is dominant, and the extreme values are reached at the same α with an accuracy of 0.001, while in Case 4 (Figure 12), PCV1 is non-dominant and highly ambiguous because the extreme values are achieved for the maximum at 0.3037 and 0.3291 and for the minimum at 0.3288 and at the active point at 0.2932 with an accuracy of 0.0001.In particular, Case 2 shows that CV cannot be determined only analytically, as this could lead to the incorrect conclusion that instability can only occur for α > 2.388.This means that when the number of resonances is less than five, an additional parametric analysis must be performed.However, once the PCVs are determined, it is difficult to discern their dominance, i.e., their direct influence on the instability lines.In Case 3 (Figure 11), PCV1 is dominant, and the extreme values are reached at the same α with an accuracy of 0.001, while in Case 4 (Figure 12), PCV1 is non-dominant and highly ambiguous because the extreme values are achieved for the maximum at 0.3037 and 0.3291 and for the minimum at 0.3288 and at the active point at 0.2932 with an accuracy of 0.0001.In particular, Case 2 shows that CV cannot be determined only analytically, as this could lead to the incorrect conclusion that instability can only occur for α > 2.388.This means that when the number of resonances is less than five, an additional parametric analysis must be performed.However, once the PCVs are determined, it is difficult to discern their dominance, i.e., their direct influence on the instability lines.Further analysis will show the instability lines and their connection to the CVs in both regular and irregular cases.Cases 1 and 2 are visualized in Figures 13 and 14.In Figures 13 and 14, features of the regular behavior that is consistent with the one-and two-layer models can be observed: (i) the instability branches are contained in the three regions delimited by the vertical lines of CV or PCV; (ii) the instability lines for lower damping levels are below the lines for higher damping, and they do not cross; (iii) there is only one instability branch in each region; (iv) with decreasing damping, the instability lines approach the CV or PCV faster from the left than from the right; and (v) instability branches tend In Case 3 (Figure 11), PCV1 is dominant, and the extreme values are reached at the same α with an accuracy of 0.001, while in Case 4 (Figure 12), PCV1 is non-dominant and highly ambiguous because the extreme values are achieved for the maximum at 0.3037 and 0.3291 and for the minimum at 0.3288 and at the active point at 0.2932 with an accuracy of 0.0001.In particular, Case 2 shows that CV cannot be determined only analytically, as this could lead to the incorrect conclusion that instability can only occur for α > 2.388.This means that when the number of resonances is less than five, an additional parametric analysis must be performed.However, once the PCVs are determined, it is difficult to discern their dominance, i.e., their direct influence on the instability lines.
Further analysis will show the instability lines and their connection to the CVs in both regular and irregular cases.Cases 1 and 2 are visualized in Figures 13 and 14.In Figures 13 and 14, features of the regular behavior that is consistent with the one-and two-layer models can be observed: (i) the instability branches are contained in the three regions delimited by the vertical lines of CV or PCV; (ii) the instability lines for lower damping levels are below the lines for higher damping, and they do not cross; (iii) there is only one instability branch in each region; (iv) with decreasing damping, the instability lines approach the CV or PCV faster from the left than from the right; and (v) instability branches tend asymptotically to infinite η M because q M is tending to zero, except for the last region, where the right-hand end tends to zero η M for infinite α.In addition, the first two regions can be left without any branches, as in Figure 6.Furthermore, FCV has no effect on the instability lines, but this property is preserved even in irregular cases.In Case 1, there is only one exception for very low damping.The instability line for η p = η b = η f = 1 × 10 −10 adopts a strange shape in the first region, which will be discussed in the next section.Further analysis will show the instability lines and their connection to the CVs in both regular and irregular cases.Cases 1 and 2 are visualized in Figures 13 and 14.In Figures 13 and 14, features of the regular behavior that is consistent with the one-and two-layer models can be observed: (i) the instability branches are contained in the three regions delimited by the vertical lines of CV or PCV; (ii) the instability lines for lower damping levels are below the lines for higher damping, and they do not cross; (iii) there is only one instability branch in each region; (iv) with decreasing damping, the instability lines approach the CV or PCV faster from the left than from the right; and (v) instability branches tend asymptotically to infinite ηM because qM is tending to zero, except for the last region, where the right-hand end tends to zero ηM for infinite α.In addition, the first two regions can be left without any branches, as in Figure 6.Furthermore, FCV has no effect on the instability lines, but this property is preserved even in irregular cases.In Case 1, there is only one exception for very low damping.The instability line for ηp = ηb = ηf = 1 × 10 −10 adopts a strange shape in the first region, which will be discussed in the next section.Further, Cases 3 and 4 are irregular, which means that properties (i)-(v) are not preserved.This is demonstrated in Figures 15 and 16.
In Figure 15, the α-scale is not extended up to CV3 for better clarity.On the other hand, for the same reason, the α-scale does not cover the non-dominant PCV in Figure 16.In both cases, the instability branches can be seen to intersect CV and cross each other, and at low damping levels, strange shapes and additional branches appear.Here, in regard to (v), the instability lines can tend to infinite η M due to K(0, q) tending to zero, which is a new feature.
In summary, by fixing a particular moving mass ratio, meaning by tracing a horizontal line in the graphs of Figures 6, 7a, 13a, 14a, 15 and 16, the dependence of the critical velocity on the damping level is demonstrated.It is seen that in regular cases, there is a smooth increase in the critical velocity with the damping level; however, this dependence can be suddenly interrupted for a certain damping value or changed in irregular cases.Some of the irregularities will be discussed further in the next section.Further, Cases 3 and 4 are irregular, which means that properties (i)-(v) are not preserved.This is demonstrated in Figures 15 and 16  Further, Cases 3 and 4 are irregular, which means that properties (i)-(v) are not preserved.This is demonstrated in Figures 15 and 16  Further, Cases 3 and 4 are irregular, which means that properties (i)-(v) are not preserved.This is demonstrated in Figures 15 and 16

Influence of the Damping of the Materials in the Foundation
In this section, the effect of damping in foundation materials is discussed.This issue is very sensitive to the critical velocity of a moving load [78,79].The critical velocity of one moving mass in the undamped case is expected to match the critical velocities of the moving force.Therefore, the instability lines for very low damping values should be close to the vertical lines indicating CVs or PCVs.However, it has been observed that very low damping values can lead to rather peculiar forms of the instability branches.Additionally, there are more branches than predicted by the number of regions, and some of the new branches can even be closed.
Figure 17 shows the initial part of the first branch for two damping levels in Case 1.The η M -scale in Figure 17 is extended to academic values to demonstrate that even though the properties of the regular behavior are preserved, there is no consistency in the sense that the critical velocity of the moving mass varies smoothly as the damping decreases.

Influence of the Damping of the Materials in the Foundation
In this section, the effect of damping in foundation materials is discussed.This issue is very sensitive to the critical velocity of a moving load [78,79].The critical velocity of one moving mass in the undamped case is expected to match the critical velocities of the moving force.Therefore, the instability lines for very low damping values should be close to the vertical lines indicating CVs or PCVs.However, it has been observed that very low damping values can lead to rather peculiar forms of the instability branches.Additionally, there are more branches than predicted by the number of regions, and some of the new branches can even be closed.
Figure 17 shows the initial part of the first branch for two damping levels in Case 1.The ηM-scale in Figure 17 is extended to academic values to demonstrate that even though the properties of the regular behavior are preserved, there is no consistency in the sense that the critical velocity of the moving mass varies smoothly as the damping decreases.
Figure 17 shows that for realistic ηM, there is a sudden jump, and two more closed intervals of instability appear at a specific damping value, covering α-values significantly lower than the previous ones.Although the exponential growth of unstable behavior is very slow in such cases, it will lead to instability.The presented results are calculated with very high precision achieved by a high number of significant digits, which is possible in symbolic software.Figure 17 shows that for realistic η M , there is a sudden jump, and two more closed intervals of instability appear at a specific damping value, covering α-values significantly lower than the previous ones.Although the exponential growth of unstable behavior is very slow in such cases, it will lead to instability.The presented results are calculated with very high precision achieved by a high number of significant digits, which is possible in symbolic software.
Figures 18 and 19 show the instability branches for low damping levels for Cases 3 and 4, again with the η M -scale extended to academic values.
It can be observed in Figure 18 that, especially at low damping levels, the instability lines have rather complicated shapes.There are six, six, and seven branches for In Figure 19, it is possible to observe how the number of instability branches increases with decreasing damping.Although the η M -scale is academic, there are other branches not covered by the scale.In Figure 19, PCV1, which was identified as non-dominant and highly ambiguous (between 0.2932 and 0.3291), is not plotted, but it can be seen that it serves to delimit the first region where instability branches appear.However, instability is unrealistic in this case due to very high η M .3).
In Figure 19, it is possible to observe how the number of instability branches increases with decreasing damping.Although the ηM-scale is academic, there are other branches not covered by the scale.In Figure 19, PCV1, which was identified as non-dominant and highly ambiguous (between 0.2932 and 0.3291), is not plotted, but it can be seen that it serves to delimit the first region where instability branches appear.However, instability is unrealistic in this case due to very high ηM.3).
There are five, six, and six branches for   3).

Comparison between the Results Obtained Using the Semi-Analytical Method and Green's Function Method
In this section, some results regarding the critical velocity of the moving mass and its time evolution obtained using the semi-analytical method and the method based on the Green's function are presented.In fact, several of the results presented in the previous sections have been verified using the D-decomposition method based on the Green's function, and to avoid redundant presentation, only one case is presented here.3).

Comparison between the Results Obtained Using the Semi-Analytical Method and Green's Function Method
In this section, some results regarding the critical velocity of the moving mass and its time evolution obtained using the semi-analytical method and the method based on the Green's function are presented.In fact, several of the results presented in the previous sections have been verified using the D-decomposition method based on the Green's function, and to avoid redundant presentation, only one case is presented here.
Figure 20 shows the D-decomposition diagrams determining the critical velocity in Case 3 for η p = η b = η f = 0.01 and η M = 60.When the moving mass velocity is 0.46369, the system is stable, and when the velocity is 0.46370, the system is unstable, meaning that the critical velocity is between these two values.This result corresponds to the one displayed in Figure 15 (yellow curve).
Figure 20 shows the D-decomposition diagrams determining the critical veloc Case 3 for ηp = ηb = ηf = 0.01 and ηM = 60.When the moving mass velocity is 0.4636 system is stable, and when the velocity is 0.46370, the system is unstable, meaning th critical velocity is between these two values.This result corresponds to the one disp in Figure 15 (yellow curve).When calculating the time series of the moving mass and the beam using the rithm based on the Green's function presented in Section 4.2, the result depends o stiffness of the elastic element inserted between the moving mass and the beam.
Figure 21 shows the time series of the beam at the moving contact point at a d sionless velocity of 0.46369, considering two values for the dimensionless stiffness elastic element, namely, κc = 94.25049 and κc = 9425.049.The calculation was perform the MATLAB environment using a laptop with an Intel Core i5-8265U CPU @ 1.60 1.80 GHz microprocessor and installed RAM of 8.00 GB.The calculation time is ap mately 7850 s, of which 6610 s is for the calculation of the Green's function and t maining 1240 s is for the calculation of the convolution integral with the determinat the quantities of interest: the displacement of the moving mass and the beam and th tact force.The time step used was Δτ = 0.0258.When calculating the time series of the moving mass and the beam using the algorithm based on the Green's function presented in Section 4.2, the result depends on the stiffness of the elastic element inserted between the moving mass and the beam.
Figure 21 shows the time series of the beam at the moving contact point at a dimensionless velocity of 0.46369, considering two values for the dimensionless stiffness of the elastic element, namely, κ c = 94.25049 and κ c = 9425.049.The calculation was performed in the MATLAB environment using a laptop with an Intel Core i5-8265U CPU @ 1.60 GHz 1.80 GHz microprocessor and installed RAM of 8.00 GB.The calculation time is approximately 7850 s, of which 6610 s is for the calculation of the Green's function and the remaining 1240 s is for the calculation of the convolution integral with the determination of the quantities of interest: the displacement of the moving mass and the beam and the contact force.The time step used was ∆τ = 0.0258.
Figure 20 shows the D-decomposition diagrams determining the critical veloc Case 3 for ηp = ηb = ηf = 0.01 and ηM = 60.When the moving mass velocity is 0.4636 system is stable, and when the velocity is 0.46370, the system is unstable, meaning th critical velocity is between these two values.This result corresponds to the one disp in Figure 15 (yellow curve).When calculating the time series of the moving mass and the beam using the rithm based on the Green's function presented in Section 4.2, the result depends o stiffness of the elastic element inserted between the moving mass and the beam.
Figure 21 shows the time series of the beam at the moving contact point at a d sionless velocity of 0.46369, considering two values for the dimensionless stiffness elastic element, namely, κc = 94.25049 and κc = 9425.049.The calculation was perform the MATLAB environment using a laptop with an Intel Core i5-8265U CPU @ 1.60 1.80 GHz microprocessor and installed RAM of 8.00 GB.The calculation time is app mately 7850 s, of which 6610 s is for the calculation of the Green's function and t maining 1240 s is for the calculation of the convolution integral with the determinat the quantities of interest: the displacement of the moving mass and the beam and th tact force.The time step used was Δτ = 0.0258.When the dimensionless stiffness of the moving mass/beam contact is low (Figure 21a), the system is unstable, but for sufficiently high dimensionless stiffness, the time series shows that the system is stable (Figure 21b).In the semi-analytical approach, fixed contact is assumed; therefore, in this case, the vibrations induced by the moving mass are stable.These results are included in Figure 21b, which shows that the two curves overlap.For a numerical evaluation using the semi-analytical approach, it is necessary to make several choices.First, given that K(0, q), like the full function to be integrated into Equation (41), has a symmetric real part and an antisymmetric imaginary part when q ∈ (−ia − ∞, −ia + ∞), and by exploiting Euler's formula with q = −ia + q r , the following holds: w(0, τ) = 1 π −ia+∞ −ia−∞ −2iη C K(0,q) q(π−2η M q 2 K(0,q)) e iqτ dq = ∞ 0 2e aτ Re −2iη C K(0,q) qπ(π−2η M q 2 K(0,q)) cos(q r τ) − Im −2iη C K(0,q) qπ(π−2η M q 2 K(0,q)) sin(q r τ) dq r .
(120) It is clear from Equation (120) that a, which is positive, must be chosen not only with respect to the definition of the inverse Laplace transform but also with respect to the time interval that is considered for the evaluation.When performing calculations in MATLAB, the numerical precision is compromised at very high τ due to the inherent double precision of this software.Unfortunately, the numerical evaluation of two embedded integrals with exact infinite limits is extremely time-consuming and does not allow any simplification.Therefore, it is best to use a simple trapezoidal rule for the evaluation of the integral in Equation (120).To this end, all values not involving τ can be prepared and then used with any chosen τ.This means that no previous time results are necessary, and thus, the only criterion for the time step is the smoothness of the resulting curve.In this case, ∆τ = 5 was found to be sufficient.In this approach, the calculation time is dependent on the number of values used in the trapezoidal rule, because after that, the evaluation is practically instantaneous.
When performing calculations in MATLAB, it is more convenient to calculate the K-function values numerically using a predefined MATLAB integration subroutine.This is because, for contour integration, double precision could compromise the pole and residue evaluation.The calculation time is practically the same if the p-scale is extended to infinity or some reasonable value.Therefore, the most important choices concern the scale and the step of q r .Generally, the smaller the a-value, the larger the concentration of the function to be integrated around q r = 0, i.e., around q = −ia.This allows for a shorter q r -interval but requires a lower q r -step.After calculating some values, their graph can indicate whether the choices are adequate or not.Numerically, each choice should be confirmed by convergence analysis.One of the criteria to be met, in addition to being independent of refined choices, is a zero value at zero time.For the graph in Figure 21b, a = 0.0001, q r -step= 5 × 10 −5 , and q r limited by 2.5 were chosen, which involved the evaluation of 50,000 values using the trapezoidal rule and took 947 s on a laptop with the same characteristics as specified before.
To evaluate the accuracy of the results given in Figure 21b, the RMS value of the difference between the two graphs was evaluated using discrete values with a time step of ∆τ = 5, that is, using 1001 values, and the result was 0.44%, which is very good.
To conclude, Case 4, with a low level of damping η p = η b = η f = 1 × 10 −6 , a moving mass ratio η M = 60, and two values of the velocity ratio, α = 0.4517 and α = 0.4518, is selected to demonstrate how the results obtained by the Green's function method approach those obtained by the semi-analytical method by increasing the stiffness of the contact spring.The results are presented in Figure 22. α = 0.4517 belongs to the stable interval, and α = 0.4518 belongs to the unstable one, with the exact separation being α = 0.451710491628607, but since the degree of instability is very low, the difference would only be seen after a very long time.As expected, by increasing the contact stiffness, the contact becomes closer to rigid, and the results of both methods are identical.α = 0.4517 belongs to the stable interval, and α = 0.4518 belongs to the unstable one, with the exact separation being α = 0.451710491628607, but since the degree of instability is very low, the difference would only be seen after a very long time.As expected, by increasing the contact stiffness, the contact becomes closer to rigid, and the results of both methods are identical.

Conclusions
In this paper, a detailed analysis of an infinite beam on a three-layer viscoelastic foundation subjected to a moving mass is presented.Three main objectives are listed: (i) the identification of the lowest velocity at the stability limit; (ii) the analysis of the damping influence of the foundation materials; and (iii) a comparison of the methods implemented in this paper, namely, the semi-analytical approach and the approach based on the Green's function method.
Using the Green's function method to calculate the time series of a moving mass on an infinite beam on a three-layer foundation requires elastic contact between the moving mass and the beam and the presence of damping in the foundation materials to ensure numerical stability.The advantage of this method is the possibility of taking into account the nonlinearity of the contact between the moving mass and the beam, namely, the nonlinear character of the elastic contact and the possibility of contact loss [32,40,42].On the other hand, the D-decomposition method developed in terms of the Green's function method can be applied to calculate the critical velocity of the moving mass on a beam on a three-layer foundation with high accuracy when the foundation is damped.
The main advantages of the semi-analytical approach are listed as follows: (i) there are no numerical problems in cases without any damping; (ii) there are no numerical issues in considering a rigid contact (however, a linear contact spring can also be included); and (iii) many related results, such as the critical velocity of the moving force, the analysis

Conclusions
In this paper, a detailed analysis of an infinite beam on a three-layer viscoelastic foundation subjected to a moving mass is presented.Three main objectives are listed: (i) the identification of the lowest velocity at the stability limit; (ii) the analysis of the damping influence of the foundation materials; and (iii) a comparison of the methods implemented in this paper, namely, the semi-analytical approach and the approach based on the Green's function method.
Using the Green's function method to calculate the time series of a moving mass on an infinite beam on a three-layer foundation requires elastic contact between the moving mass and the beam and the presence of damping in the foundation materials to ensure numerical stability.The advantage of this method is the possibility of taking into account the nonlinearity of the contact between the moving mass and the beam, namely, the nonlinear character of the elastic contact and the possibility of contact loss [32,40,42].On the other hand, the D-decomposition method developed in terms of the Green's function method can be applied to calculate the critical velocity of the moving mass on a beam on a three-layer foundation with high accuracy when the foundation is damped.
The main advantages of the semi-analytical approach are listed as follows: (i) there are no numerical problems in cases without any damping; (ii) there are no numerical issues in considering a rigid contact (however, a linear contact spring can also be included); and (iii) many related results, such as the critical velocity of the moving force, the analysis on finite beams, etc., are presented in analytical form.There are just no limitations; the results are highly accurate, and the formulations developed in this paper can be easily used by other researchers.Unfortunately, the main advantage of the new approach ( [10,34,35,38,59]), which is so effective in one-and two-layer models, is lost in the three-layer model, because the mass-induced frequencies are not well defined in most cases, and therefore, the analytical expression for the unsteady harmonic part cannot be used in such cases.A regularization approach to overcome this difficulty is currently in preparation.This is why the time series in this paper are presented with the numerical evaluation of the inverse Laplace transform.However, the accuracy and calculation times have been shown to be very similar to the Green's function method.In summary, time series in simpler models, like one-or two-layer models, were much simpler and more advantageous because it was possible to identify the so-called unsteady harmonic part, which was determined analytically, and thus, the numerical evaluation was practically instantaneous and could be extended to infinite time without loss of accuracy.However, other advantages of the semi-analytical approach related to the determination of instability lines are preserved for all layered models, and these lines can be determined very quickly and accurately while keeping the main part of the calculations in the real domain.
(i) The beam is straight and prismatic, and it is made of isotropic homogeneous material.(ii) The beam can withstand an axial force, in accordance with Figure 2. (iii) The beam obeys the linear elastic Euler-Bernoulli theory.(iv) Vertical displacements are measured from the equilibrium position corresponding to the deflection induced by the weight of the model components.(v) The initial conditions are homogeneous; nevertheless, this has no effect on the critical velocity.(vi) The velocity of the moving mass determines its horizontal position.(vii) No friction acts at the contact point.(viii) Loads and vertical displacements are assumed to be positive when acting downward.(ix) As is usual in several applications, the acting force may or may not represent the moving mass weight.Materials 2024, 17, x FOR PEER REVIEW 7 of 41

Figure 2 .
Figure 2. Railway track model with three viscoelastic layers subjected to an axial force and a uniformly moving mass subjected to a vertical force.

Figure 2 .
Figure 2. Railway track model with three viscoelastic layers subjected to an axial force and a uniformly moving mass subjected to a vertical force.

Figure 9 .
Figure 9. Parametric analysis related to Case 1 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a-c) correspond to different α-scales).

Figure 9 .
Figure 9. Parametric analysis related to Case 1 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a-c) correspond to different α-scales).

Figure 10 .
Figure 10.Parametric analysis related to Case 2 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a-c) correspond to different α-scales).

Figure 11 .
Figure 11.Parametric analysis related to Case 3 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a,b) correspond to different α-scales).

Figure 10 .
Figure 10.Parametric analysis related to Case 2 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a-c) correspond to different α-scales).

Materials 2024 , 41 Figure 10 .
Figure 10.Parametric analysis related to Case 2 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a-c) correspond to different α-scales).

Figure 11 .
Figure 11.Parametric analysis related to Case 3 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a,b) correspond to different α-scales).

Figure 11 . 41 Figure 12 .
Figure 11.Parametric analysis related to Case 3 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a,b) correspond to different α-scales).Materials 2024, 17, x FOR PEER REVIEW 29 of 41

Figure 12 .
Figure 12.Parametric analysis related to Case 4 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a,b) correspond to different α-scales).

Figure 12 .
Figure 12.Parametric analysis related to Case 4 (max/min-global maximum/minimum deflection over the whole beam; AP-deflection at the active point; (a,b) correspond to different α-scales).
Figure17shows that for realistic η M , there is a sudden jump, and two more closed intervals of instability appear at a specific damping value, covering α-values significantly lower than the previous ones.Although the exponential growth of unstable behavior is very slow in such cases, it will lead to instability.The presented results are calculated with very high precision achieved by a high number of significant digits, which is possible in symbolic software.Figures18 and 19show the instability branches for low damping levels for Cases 3 and 4, again with the η M -scale extended to academic values.It can be observed in Figure18that, especially at low damping levels, the instability lines have rather complicated shapes.There are six, six, and seven branches forη p = η b = η f = 1 • 10 −4 , η p = η b = η f = 1 • 10 −5, and η p = η b = η f = 1 • 10 −10 , respectively.Each mentioned case has one of the branches closed.For lower damping, these branches attain realistic values of η M .In Figure19, it is possible to observe how the number of instability branches increases with decreasing damping.Although the η M -scale is academic, there are other branches not covered by the scale.In Figure19, PCV1, which was identified as non-dominant and highly ambiguous (between 0.2932 and 0.3291), is not plotted, but it can be seen that it serves to delimit the first region where instability branches appear.However, instability is unrealistic in this case due to very high η M .

Figures 18 and 19
Figures 18 and 19 show the instability branches for low damping levels for Cases 3 and 4, again with the ηM-scale extended to academic values.

Figure 21 .
Figure 21.Time series of the beam at the moving point for (a) Green's function method wi 94.25049;(b) Green's function method with κc = 9425.049(GFM) and the semi-analytical app (SAM).

Figure 21 .
Figure 21.Time series of the beam at the moving point for (a) Green's function method wi 94.25049;(b) Green's function method with κc = 9425.049(GFM) and the semi-analytical app (SAM).

Figure 21 .
Figure 21.Time series of the beam at the moving point for (a) Green's function method with κ c = 94.25049;(b) Green's function method with κ c = 9425.049(GFM) and the semi-analytical approach (SAM).

Table 1 .
Ranges of typical values of the main parts of the three-layer model.

Table 3 .
Other selected cases to exemplify the influence of the resonances and the classification.

Table 3 .
Other selected cases to exemplify the influence of the resonances and the classification.First, the effect of the number of resonances is demonstrated in Figures 9-12