Instability of a Moving Bogie: Analysis of Vibrations and Possibility of Instability in Subcritical Velocity Range

1. Introduction

One of the effective ways to combat climate change and reduce the carbon footprint is to prioritize rail over road transport. This naturally implies the need to increase the railway network capacity by increasing travel speed, number of trains and axle loads. Such increased demands emphasize the need of ensuring the passenger safety and comfort. The upper limit for the travel speed is usually derived from the value of the critical velocity, which in this context corresponds to the critical velocity of the moving force, which is equal to the lowest velocity of wave propagation in the supporting structure, resulting from the dynamic interaction of all the parts involved and can be called the critical velocity for resonance.

Another commonly used indicator for ensuring the safety of train passage is the receptance curve, which, however, only indicates the natural frequencies of the supporting structure. The frequencies that cause unstable behavior cannot be identified using the receptance curves because they depend on the mass and velocity of the moving object. Therefore, it is necessary to study the transient vibrations caused by the moving vehicle and the associated so-called induced frequencies, which can identify whether the motion will be stable or unstable. Nevertheless, the problem of instability was somehow overlooked because it was shown that a single moving mass or oscillator can achieve unstable behavior only in the supercritical velocity range, and since this region is a priori excluded from the design of the track, the issue of instability was not found to be crucial.

However, the dynamic interaction of two moving proximate masses can lead to a different situation, which surprisingly worsens with increased foundation damping. Such cases have already been identified in the author’s previous works.

Moving load problems are still very active scientific topics related to transportation engineering applied to rail or road transport. Several simplifications have been introduced to solve these problems more efficiently and to present the results in a more compact form suitable for optimization and sensitivity analyses. The simplifications are applied to the supporting structure or to the moving vehicle, or both, depending on the problem being solved. In addition to the analysis of the influence of several parameters on dynamic performance, topics of particular interest are critical velocity for resonance, critical frequency for resonance, instability and dynamic interaction of proximate inertial objects.

Typical simplifications at the level of supporting structure led to layered models consisting of beams, point masses and spring-damper elements. These models are quite popular for their simplicity and computational efficiency. Such simplified models, together with simplified vehicle models, are used in more practical studies as well as in some more theoretical studies. The applicability of the three-layer model in railway applications is analyzed in [1] and an overview of some features of the one-, two- and three-layer models is given in [2].

Many research works have been developed during the last decades; thus, this review will essentially focus on the most recent publications.

In [3] a critical discussion of developments in the field of dynamic response of rigid road pavements to a moving vehicle is presented. Rigid pavements consist of an elastic concrete plate or jointed plates resting on a supporting foundation. Such models are also used in railway engineering, only with different parameters. The supporting medium can be modeled as a system of spring-damper elements or as a homogeneous or layered half-plane/space with infinite or finite depth. The paper summarizes the works where the material behavior of the subgrade layers is isotropic as well as anisotropic, viscoelastic, poroelastic or inelastic. In [4] another more detailed review by the same authors on both types of pavements, i.e. rigid and flexible, is given. In these works, vehicles are simulated by constant concentrated or distributed over a finite line/area moving forces, while their inertial effect and thus the possibility of instability is omitted. In [5], the plane strain model is used in the moving mass problem to analyze instability issues.

In [6], although also focused on road pavements, the guiding structure model corresponds to the one-layer model, which is often used in railway applications. The model is crossed by a two-mass oscillator with a contact spring-damper element, also very common in railway applications. The harmonic roughness of the roadway is taken into account, but no instability or critical velocity are analyzed. Since there is only one contact point, the dynamic interaction of proximate objects cannot be observed. The work focuses on the analysis of the effect of the nonlinearity of the foundation, which is addressed by the Adomian decomposition method. Moving coordinates are implemented and the results are limited to the steady-state part of the response, since only a single Fourier transform is applied. Such an approach in fact removes the inertial effects of the moving object. The research findings are that there is a critical value for the nonlinear stiffness, which increases with the vehicle velocity and the foundation damping.

As for the railway applications, the nonlinearity of the supporting structure is analyzed in [7]. A two-layer model is implemented, and a simplified vehicle model is used. A new nonlinear stiffness model that comprehensively characterizes the ballast properties by a displacement-dependent stiffness, frequency-dependent stiffness, hysteresis, and time/space-varying features is developed. The model is validated by experimental measurements processing the ground penetrating radar signal using the adaptive optimal kernel time–frequency representation method.

Moving force problems are simpler than moving mass problems, especially for finite length supporting structures, for which analytical solutions exist and can be formalized by modal expansion. However, a relatively recent study [8] still analyzes this problem for a double-beam structure focusing on the effect of various boundary conditions. Modal expansion is efficient on finite structures, and although the motion of mass on finite Euler-Bernoulli beams without elastic foundation and damping is a classic problem that has already been addressed by several researchers, there are still new results and physical insights worth publishing for such a simple structure as shown in [9]. Moreover, in [10] a new equivalent lumped parameter model is proposed to describe the vibrations of simply supported beams without an elastic foundation under the action of a moving force, which is a classic problem with a fully analytical solution. In [11], vibrations of a similar structure subjected to a moving mass are solved using Green’s function method, with the aim to analyze the effect of boundary conditions, mass velocity and its ratio with respect to the beam mass.

More practical studies usually deal with typical degradation on railway lines, that is transition zones, rail corrugation, maintenance optimization or possible failures. For example, in [12] the effect of temperature on concrete cracking in unballasted lines is analyzed; in [13] failure of suspension damper in railway vehicles based on the cross-correlation analysis of bogie accelerations is presented.

Regarding the transition zones, it is known that the sudden change in the vertical stiffness of railway tracks increases the dynamic loads and causes numerous defects on ballasted tracks. The main geometric parameters of the approach slab for short-span bridges are optimized in [14]. A finite element model is built and subjected to a series of moving forces. The model was validated with field measurements acquired using a laser/camera-based measuring technique. The obtained geometric sensitivity analysis results show that the rail and ballast displacements can be significantly reduced by using a geometrically optimized approach slab along the transition zone of a railway bridge.

In [15], a revolutionary new design recently proposed to mitigate the negative effects of transition zones is evaluated for extreme conditions. This new design of a transition structure is called SHIELD (Safe Hull Inspired Energy Limiting Design), and was designed to reverse the propagating waves by inverting the shape of the embankment zone, thus directing the wave energy towards the foundation instead of upwards to damage the railroad. It was found that the novel design effectively mitigates dynamic amplifications and leads to a smooth distribution of deformation energy across subcritical, critical and supercritical velocity regimes in both directions of travel, meaning that the expected traffic-induced degradation will be as uniform as possible in the longitudinal direction.

Transition zones are directly linked to the floating sleeper phenomenon, for which a novel approach in modelling is presented in [16], providing an accurate mathematical model with a minimum number of degrees of freedom supported by extensive field measurements of a railway track. Although the observed phenomenon is very complex, a simplified mathematical model of one degree of freedom has proven to be sufficiently accurate for its characterization.

Regarding rail corrugation, a new hypothesis for the formation process of short pitch rail corrugation is presented in [17]. A dynamic finite element (FE) wheel-track model is used to verify the hypothesis by reproducing corrugation initiation and its consistent growth. The FE model is composed of discrete elements as in other works with semianalytical or analytical solutions. The supporting structure is in fact equivalent to a two-layer model extended to 3D, and both the railpads and the foundation, which in fact only corresponds to the ballast layer, are modeled by a set of spring-damper elements acting in three directions. The wheelset and rails are modeled in detail using 3D elements, and then the remaining mass of the vehicle is lumped and placed above the primary suspension. It is found that longitudinal compression modes are responsible for corrugation initiation. The proposed mechanism can explain field observations. In [18] the finite element model is validated using a downscale V-Track test rig.

Transition zones are related to transition radiation. In this context it is worthwhile to mention [19,20].

As for the maintenance, ballast tamping is still the most common way of restoring ballast bed geometry. In [21], a novel test device is used that allows real-time evaluation of the mechanical qualities of ballast beds, which is essential for improving tamping efforts in railway maintenance. The evaluation is based on the impact force on the tamper, which this new device can measure. The results have also been validated numerically using a coupled model where the ballast is approached using the discrete element method.

In [22], closed-form solution is derived by the residue theorem and Green’s function for the dynamic response of an infinite beam resting on a viscoelastic foundation under a harmonic moving load. However, the closed-form results have already been presented [23]. In [22], these results are used to define the most appropriate span length for an analogous finite beam to best represent the deflection of an infinite beam for practical reasons. As shown in [24], this can be a problem for a constant moving force and the solution to this problem is proposed there, but not for a moving inertial load, as has been shown in previous works by the author of this paper, where the validation was always carried out on finite beams, the only precaution being to start the load a little further from the support, as explained in [25]. With the initial conditions, this distance must be longer, since the entire relevant part of initial shape [26] must be contained within the beam at zero time.

Similar models are in [27], but only the steady state vibrations are analyzed there, so that the possibility of instability is a priory removed. However, the model is more detailed because rail roughness and contact springs are introduced. Rail corrugation effect is also included in [28].

It is also necessary to mention other works analyzing the limitations of classical beam theories [29,30,31], large deformations [32], the non-uniform motion of the inertial object [33] or analyses concerned with tracks other than classical ballasted ones [34].

The stability of a moving inertial object on an elastic structure with railway applications has already been investigated by considering different representations of the vehicle (moving object) and the track (supporting elastic structure or guiding structure). Simple vehicle models include a moving mass [35,36], a two-mass oscillator [23,37] and a three-mass oscillator [38]. These models have only one contact point, so no dynamic interaction can be observed, as when considering a set of moving inertial objects in the form of two moving masses [25,39] or a sequence of two-mass oscillators [40,41,42,43]. A more complex vehicle model can also be used to study stability issues [44]. The guiding structure is modeled by an Euler-Bernoulli infinite beam on a continuous foundation with one elastic layer [36,38,39,40] or two [25,45] or three elastic layers [2]. Instead of the Euler–Bernoulli beam, the Timoshenko beam can also be used [38,46]. In order to take into account the influence of the sleeper bay, an infinite beam model on a foundation with discrete supports should be considered for the track [47].

The problem of instability of a moving bogie has already been analyzed by other researchers, however, none of these published works present the complete vibration shapes, closed-form time domain analytical solutions, alternatives to the D-decomposition method, validation on finite beams, etc. The instability of a moving bogie passing through a single-layer model is analyzed in [48] from several aspects. The paper also provides a comparison with a moving oscillator and other physical insights, in particular an analysis of the critical primary stiffness of the suspension. In [49], the authors intend to demonstrate that for a more complex foundation structure the critical stiffness changes its value, however only some isolated values are given and no detailed comparison with [48] is provided. The critical velocity for resonance is obviously much different than for the one-layer model, so the comparison should be organized differently. The same reference input data of the moving bogie are used as in [48], thus due to the large distance between the wheels the dynamic interaction between them is suppressed and the instability is only detectable for realistic values in the supercritical velocity range. In [50] the authors introduced more complicated moving objects, but again, results are restricted to few discrete values of the critical primary suspension stiffness. The aggravation by foundation damping and instability in subcritical velocity range have not been identified in any of these works.

In the present paper, a moving bogie is analyzed with the aim of clarifying the possibility of its instability in the subcritical velocity region. The results are presented in as much as analytical way possible, therefore the previous developments presented in [23,36,39] will be exploited. It has recently been shown that instability of two moving proximate masses passing through a layered track model can occur in the subcritical range of velocities [25,39]. This contradicts the generally accepted idea that unstable velocities are always higher than the lowest critical velocity for resonance. While for a single moving inertial object, where there is only one contact point between the object and the supporting structure, the possibility of instability is governed by the equivalent stiffness at the contact point, for a set of moving inertial objects this is no longer true due to their dynamic interaction. Although it is mathematically impossible to obtain unstable behavior below the critical velocity for a single moving object, for multiple moving objects this can already occur. Moreover, it has been shown that the situation is generally worsened by increased foundation damping, which is not expected. Thus, the novel contribution of this paper lies in the analysis of how the above behavior changes by combining the moving masses into a single bogie, which is more realistic.

The paper is organized in the following way. First, the model is described in Section 2 together with the governing equations and their solution, exploiting the new method from [36] for presenting the results in closed-form formulas. Dimensionless parameters are used to cover a wide range of possible scenarios and to allow the results to be used without having to redo the calculations. The results Section 3 is divided into several parts, specifically the induced frequencies are defined and then the results validation by the analysis on finite beams is described. Furthermore, the changes caused by switching from the Euler-Bernoulli beam to the Timoshenko-Rayleigh beam are derived and a set of selected cases to be presented further is given, together with the corresponding induced frequencies. Subsequently, the conditions for unstable cases in the subcritical range of velocities of two moving proximate masses are summarized and compared with the cases of a two-axle bogie. Then, the effect of the suspension damping, the beam theory and other input parameters are presented. The results part is concluded by comparison with other published works, namely [48]. Same cases are included in Appendix A due to the length of the paper. The paper is concluded by results discussion and finally conclusions are drawn.

2. The Model and the Governing Equations

The model considered in the present paper is depicted in Figure 1. In part (a), two proximate masses pass through the one-layer model, and in part (b), a single bogie consisting of a rigid bogie bar of finite length, primary suspension, and two point masses representing the wheels, moves across the model.

For simplicity and models comparison, it is assumed that the moving masses have the same value and are acted upon by the same constant forces. The primary suspension of the bogie is modeled by spring-damper elements and each of these supports is connected to a point mass representing the wheel, which is in permanent contact with the beam. For simplicity, the supports are considered identical and symmetrically positioned with respect to the mass center of the bar.

The problem of two moving masses has already been solved in previous work, [39]. The problem of moving bogie is governed in fixed coordinates by the following equations:

$$EI{w}_{,xxxx}+m{w}_{,tt}+c_f{w}_{,t}+k_{f}w= p\left(x,t\right)$$

(1)

$$\begin{array}{l}p\left(x,t\right)=\left(P-M{w}_{01,tt}\left(t\right)-k_v\left(w_{01}-u_v+\theta {\displaystyle \frac{d}{2}}\right)-c_v\left(w_{01,t}-u_{v,t}+{\theta }_{,t}{\displaystyle \frac{d}{2}}\right)\right)\delta \left(x-x_1\right)\\ +\left(P-M{w}_{02,tt}\left(t\right)-k_v\left(w_{02}-u_v-\theta {\displaystyle \frac{d}{2}}\right)-c_v\left(w_{02,t}-u_{v,t}-{\theta }_{,t}{\displaystyle \frac{d}{2}}\right)\right)\delta \left(x-x_2\right)\end{array}$$

(2)

$$M_b{u}_{v,tt}\left(t\right)-k_v\left(w_{01}+w_{02}-2u_v\right)-c_v\left(w_{01,t}+w_{02,t}-2u_{v,t}\right)=0$$

(3)

$$J_b{\theta }_{,tt}+k_v\left(w_{01}-w_{02}+\theta d\right){\displaystyle \frac{d}{2}}+c_v\left(w_{01,t}-w_{02,t}+{\theta }_{,t}d\right){\displaystyle \frac{d}{2}}=0$$

(4)

Here, Equation (1) represents the beam equilibrium, Equation (2) describes the loading term, Equation (3) is the vertical equilibrium of the bogie bar and Equation (4) is the moment equilibrium of the bogie bar around its mass center. Dependance on spatial coordinate $x$ and time $t$ is only used for clarity in the loading term $p\left(x,t\right)$ and is omitted in other terms for simplicity. Unknown fields are: the beam vertical displacement $w$, bogie bar mass center vertical displacement $u$ and rotation $\theta$ with positive orientations as indicated in Figure 1. Derivatives are designated by variables in subscript position preceded by a comma.

The model consists of a beam characterized by its bending stiffness $EI$ and mass per unit length $m$ and obeying the Euler-Bernoulli theory. The foundation is described by its vertical and shear stiffness, $k_f$ and $k_s$, and the viscous damping coefficient $c_f$. The moving bogie has a constant velocity $v$ and is composed of a rigid bar with mass $M_b$ and mass moment of inertia $J_b$, spring-damper elements of the primary suspension are characterized by $k_v$ and $c_v$ and two point masses $M$ distanced by $d$ in a symmetrical arrangement with respect to the bogie bar mass center. The model in Figure 1(a) uses only $M$ and $d$, but when appropriate, the moving masses are recalculated to maintain the same total mass as the bogie, for comparison purposes. $\delta$ is the Dirac delta function. The initial conditions are assumed to be homogeneous. The boundary conditions state that the deflection and slope vanish for $x\to \pm \infty$.

Furthermore, it is assumed that the point masses are in permanent contact with the beam. $w_{0j}\left(t\right), j=1,2$ are the displacements of the contact points, $j=1$ corresponds to the rear mass and $j=2$ to the front mass. To eliminate these additional unknowns, the correct expansion of the time derivatives must be used:

$$\begin{gathered} w_{01,t}\left(t\right)=v{w}_{,x}\left(x,t\right)+w_{,t}\left(x,t\right), w_{01,tt}\left(t\right)=v^2{w}_{,xx}\left(x,t\right)+2v{w}_{,xt}\left(x,t\right)+w_{,tt}\left(x,t\right) \\ \mathrm{with} x=x_1 \end{gathered}$$

(5)

$$\begin{gathered} w_{02,t}\left(t\right)=v{w}_{,x}\left(x,t\right)+w_{,t}\left(x,t\right), w_{02,tt}\left(t\right)=v^2{w}_{,xx}\left(x,t\right)+2v{w}_{,xt}\left(x,t\right)+w_{,tt}\left(x,t\right) \\ \mathrm{with} x=x_2 \end{gathered}$$

(6)

where $x_1=vt$, $x_2=vt+d$. Then, Equations (2), (3) and (4) have the following forms:

$$\begin{array}{l}p\left(x,t\right)=\left(P-M\left(w_{,tt}+2v{w}_{,xt}+v^2{w}_{,xx}\right)-k_v\left(w-u_v+\theta {\displaystyle \frac{d}{2}}\right)-c_v\left(v{w}_{,x}+w_{,t}-u_{v,t}+{\theta }_{,t}{\displaystyle \frac{d}{2}}\right)\right)\delta \left(x-vt\right)\\ +\left(P-M\left(w_{,tt}+2v{w}_{,xt}+v^2{w}_{,xx}\right)-k_v\left(w-u_v-\theta {\displaystyle \frac{d}{2}}\right)-c_v\left(v{w}_{,x}+w_{,t}-u_{v,t}-{\theta }_{,t}{\displaystyle \frac{d}{2}}\right)\right)\delta \left(x-vt-d\right)\end{array}$$

(7)

$$\begin{array}{l}{M}_b{u}_{v,tt}\left(t\right)-k_v\left(w\left(vt,t\right)+w\left(vt+d,t\right)-2u_v\right)\\ -c_v\left(v{w}_{,x}\left(vt,t\right)+w_{,t}\left(vt,t\right)+v{w}_{,x}\left(vt+d,t\right)+w_{,t}\left(vt+d,t\right)-2u_{v,t}\right)=0\end{array}$$

(8)

$$\begin{array}{l}{J}_b{\theta }_{,tt}+k_v\left(w\left(vt,t\right)-w\left(vt+d,t\right)+\theta d\right){\displaystyle \frac{d}{2}}\\ +c_v\left(v{w}_{,x}\left(vt,t\right)+w_{,t}\left(vt,t\right)-v{w}_{,x}\left(vt+d,t\right)-w_{,t}\left(vt+d,t\right)+{\theta }_{,t}d\right){\displaystyle \frac{d}{2}}=0\end{array}$$

(9)

Then it is convenient to switch to moving coordinates, $r=x-vt$, fixing the zero position on the rear mass and, similarly as before, it is necessary to properly expand the time derivatives

$$\begin{gathered} w_{,t}\left(x,t\right)=-v{w}_{,r}\left(r,t\right)+w_{,t}\left(r,t\right), \\ w_{,tt}\left(x,t\right)=v^2{w}_{,rr}\left(r,t\right)-2v{w}_{,rt}\left(r,t\right)+w_{,tt}\left(r,t\right) \end{gathered}$$

(10)

The system to be solved has the following final form:

$$EI{w}_{,rrrr}+m\left(w_{,tt}-2v{w}_{,rt}+v^2{w}_{,rr}\right)+c_f\left(w_{,t}-v{w}_{,r}\right)+k_{f}w=p\left(r,t\right)$$

(11)

$$\begin{array}{l}p\left(r,t\right)=\left(P-M{w}_{,tt}-k_v\left(w-u_v+\theta {\displaystyle \frac{d}{2}}\right)-c_v\left(w_{,t}-u_{v,t}+{\theta }_{,t}{\displaystyle \frac{d}{2}}\right)\right)\delta \left(r\right)\\ +\left(P-M{w}_{,tt}-k_v\left(w-u_v-\theta {\displaystyle \frac{d}{2}}\right)-c_v\left(w_{,t}-u_{v,t}-{\theta }_{,t}{\displaystyle \frac{d}{2}}\right)\right)\delta \left(r-d\right)\end{array}$$

(12)

$$M_b{u}_{v,tt}\left(t\right)-k_v\left(w\left(vt,t\right)+w\left(vt+d,t\right)-2u_v\right)-c_v\left(w_{,t}\left(vt,t\right)+w_{,t}\left(vt+d,t\right)-2u_{v,t}\right)=0$$

(13)

$$2J_b{\theta }_{,tt}+k_v\left(w\left(vt,t\right)-w\left(vt+d,t\right)+\theta d\right)d+c_v\left(w_{,t}\left(vt,t\right)-w_{,t}\left(vt+d,t\right)+{\theta }_{,t}d\right)d=0$$

(14)

Next, dimensionless parameters are introduced in the following way:

$$\chi =\sqrt[4]{{\displaystyle \frac{k_f}{4EI}}}, v_{ref}=\sqrt[4]{{\displaystyle \frac{4k_{f}EI}{m^2}}}={\displaystyle \frac{1}{\chi }}\sqrt{{\displaystyle \frac{k_f}{m}}}, w_{st}={\displaystyle \frac{P\chi }{2k_f}}, {\eta }_P={\displaystyle \frac{P}{P}}=1$$

(15)

$$\xi =\chi r, \tilde{d}=d\chi , \tau =\chi v_{ref}t=\sqrt{{\displaystyle \frac{k_f}{m}}}t,\tilde{w}={\displaystyle \frac{w}{w_{st}}},\tilde{u}={\displaystyle \frac{u_v}{w_{st}}}, \tilde{\theta }={\displaystyle \frac{\theta }{w_{st}\chi }}, \alpha ={\displaystyle \frac{v}{v_{ref}}}$$

(16)

$${\eta }_M={\displaystyle \frac{M\chi }{m}}, {\mu }_v={\displaystyle \frac{M_b\chi }{m}}, {\tilde{J}}_b={\displaystyle \frac{{\chi }^3}{m}}{J}_b, {\kappa }_v={\displaystyle \frac{k_v\chi }{k_f}}, {\eta }_v={\displaystyle \frac{c_v\chi }{2\sqrt{m{k}_f}}}, {\eta }_f={\displaystyle \frac{{\mathrm{c}}_f}{2\sqrt{m{k}_f}}}$$

(17)

This results in the system of three equations

$$\begin{array}{l}{\tilde{w}}_{,\xi \xi \xi \xi }+4{\alpha }^2{\tilde{w}}_{,\xi \xi }+4{\tilde{w}}_{,\tau \tau }-8\alpha {\tilde{w}}_{,\xi \tau }+8{\eta }_b\left({\tilde{w}}_{,\tau }-\alpha {\tilde{w}}_{,\xi }\right)+4\tilde{w}=\\ 4\left(2{\eta }_P-{\eta }_M{\tilde{w}}_{,\tau \tau }- {\kappa }_v\left(\tilde{w}-\tilde{u}+{\displaystyle \frac{1}{2}}\tilde{\theta }\tilde{d}\right)-2 {\eta }_v\left({\tilde{w}}_{,\tau }-{\tilde{u}}_{,\tau }+{\displaystyle \frac{1}{2}}{\tilde{\theta }}_{,\tau }\tilde{d}\right)\right)\delta \left(\xi \right)\\ +4\left(2{\eta }_P-{\eta }_M{\tilde{w}}_{,\tau \tau }- {\kappa }_v\left(\tilde{w}-\tilde{u}-{\displaystyle \frac{1}{2}}\tilde{\theta }\tilde{d}\right)-2 {\eta }_v\left({\tilde{w}}_{,\tau }-{\tilde{u}}_{,\tau }-{\displaystyle \frac{1}{2}}{\tilde{\theta }}_{,\tau }\tilde{d}\right)\right)\delta \left(\xi -\tilde{d}\right) \end{array}$$

(18)

$${\mu }_v{\tilde{u}}_{,\tau \tau }\left(\tau \right)-{\kappa }_v\left(\tilde{w}\left(\alpha \tau ,\tau \right)+\tilde{w}\left(\alpha \tau +\tilde{d},\tau \right)-2\tilde{u}\right)-2{\eta }_v\left({\tilde{w}}_{,\tau }\left(\alpha \tau ,\tau \right)+{\tilde{w}}_{,\tau }\left(\alpha \tau +\tilde{d},\tau \right)-2{\tilde{u}}_{,\tau }\right)=0$$

(19)

$$2{\tilde{J}}_b{\tilde{\theta }}_{,tt}+{\kappa }_v\left(\tilde{w}\left(\alpha \tau ,\tau \right)-\tilde{w}\left(\alpha \tau +\tilde{d},\tau \right)+\tilde{\theta }\tilde{d}\right)\tilde{d}+2{\eta }_v\left({\tilde{w}}_{,\tau }\left(\alpha \tau ,\tau \right)-{\tilde{w}}_{,\tau }\left(\alpha \tau +\tilde{d},\tau \right)+{\tilde{\theta }}_{,\tau }\tilde{d}\right)\tilde{d}=0$$

(20)

Solution is obtained by integral transforms, where firstly the Laplace transform is applied

$$\tilde{F}\left(\xi ,\overline{q}\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(\xi ,\tau \right){\mathrm{e}}^{-\overline{q}\tau }d\tau }\mathrm{with} \overline{q}=\mathrm{i}q$$

(21)

and after that the Fourier transform is applied on the beam equilibrium

$$F\left(p,\overline{q}\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\tilde{F}\left(\xi ,\overline{q}\right){\mathrm{e}}^{-\mathrm{i}p\xi }\mathrm{d}\xi }$$

(22)

leading to

$$\begin{array}{l}W\left(p,q\right)={\displaystyle \frac{1}{D\left(p,q\right)}}\left({\displaystyle \frac{8{\eta }_P}{\mathrm{i}q}}+4\left({\eta }_M{q}^2-\left({\kappa }_v+2\mathrm{i}{\eta }_{v}q\right)\right)\tilde{W}\left(0,q\right)+ 2\left({\kappa }_v+2\mathrm{i}{\eta }_{v}q\right)\left(2\tilde{U}-\tilde{\Theta }\tilde{d}\right)\right.\\ \left.+\left({\displaystyle \frac{8{\eta }_P}{\mathrm{i}q}}+4\left({\eta }_M{q}^2-\left({\kappa }_v+2\mathrm{i}{\eta }_{v}q\right)\right)\tilde{W}\left(\tilde{d},q\right)+ 2\left({\kappa }_v+2\mathrm{i}{\eta }_{v}q\right)\left(2\tilde{U}+\tilde{\Theta }\tilde{d}\right)\right){\mathrm{e}}^{-\mathrm{i}p\tilde{d}}\right)\end{array}$$

(23)

where

$$D\left(p,q\right)= p^4-4p^2{\alpha }^2-4q^2+8\alpha pq+8\mathrm{i}{\eta }_f\left(q-\alpha p\right)+4$$

(24)

and the other two equations are

$$\left({\mu }_v{q}^2-2\left( {\kappa }_v+2{\eta }_v\mathrm{i}q\right)\right)\tilde{U}+\left( {\kappa }_v+2{\eta }_v\mathrm{i}q\right)\left(\tilde{W}\left(0,q\right)+\tilde{W}\left(\tilde{d},q\right)\right)=0$$

(25)

$$\left(-2{\tilde{J}}_b{q}^2+\left({\kappa }_v+2{\eta }_v\mathrm{i}q\right){\tilde{d}}^2\right)\tilde{\Theta }+\left({\kappa }_v+2{\eta }_v\mathrm{i}q\right)\left(\tilde{W}\left(0,q\right)-\tilde{W}\left(\tilde{d},q\right)\right)\tilde{d}=0$$

(26)

The inverse Fourier transform

$$\tilde{F}\left(\xi ,\mathrm{i}q\right)={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}F\left(p,\mathrm{i}q\right){\mathrm{e}}^{\mathrm{i}p\xi }\mathrm{d}p}$$

(27)

affects only Equation (23) and yields

$$\begin{array}{l}\tilde{W}\left(\xi ,q\right)={\displaystyle \frac{1}{2\pi }}\left({\displaystyle \frac{8{\eta }_P}{\mathrm{i}q}}+4\left({\eta }_M{q}^2- \left({\kappa }_v+2\mathrm{i}q{\eta }_v\right)\right)\tilde{W}\left(0,q\right)+ 2\left({\kappa }_v+2\mathrm{i}q{\eta }_v\right)\left(2\tilde{U}-\tilde{\Theta }\tilde{d}\right)\right)K\left(\xi ,q\right)\\ +{\displaystyle \frac{1}{2\pi }}\left({\displaystyle \frac{8{\eta }_P}{\mathrm{i}q}}+4\left({\eta }_M{q}^2- \left({\kappa }_v+2\mathrm{i}q{\eta }_v\right)\right)\tilde{W}\left(\tilde{d},q\right)+2\left({\kappa }_v+2\mathrm{i}q{\eta }_v\right)\left(2\tilde{U}+\tilde{\Theta }\tilde{d}\right)\right)K\left(\xi -\tilde{d},q\right)\end{array}$$

(28)

where

$$K\left(\xi ,q\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{{\mathrm{e}}^{\mathrm{i}p\xi }}{D\left(p,q\right)}\mathrm{d}p}$$

(29)

The last step corresponds to the inverse Laplace transform and can be performed as analytically as possible, by identifying the poles of the function to be integrated and performing the contour integration as a sum of residues. To do this, firstly unknowns $\tilde{W}\left(0,q\right)$, $\tilde{W}\left(\tilde{d},q\right)$, $\tilde{U}$ and $\tilde{\Theta }$ must be solved and substituted.

The following system can be formed when $\xi =0$ and $\xi =\tilde{d}$ are substituted into Equation (28) and Equations (25) and (26) are added.

$$\left[\begin{array}{cccc}{A}_{11}& A_{12}& A_{13}& A_{14}\\ A_{21}& A_{22}& A_{23}& A_{24}\\ A_{31}& A_{32}& A_{33}& A_{34}\\ A_{41}& A_{42}& A_{43}& A_{44}\end{array}\right]\cdot \left\{\begin{array}{c}\tilde{W}\left(0,q\right)\\ \tilde{W}\left(\tilde{d},q\right)\\ \tilde{U}\\ \tilde{\Theta }\end{array}\right\}=\left\{\begin{array}{c}{\displaystyle \frac{4{\eta }_P}{\mathrm{i}q}}K\left(0,q\right)+{\displaystyle \frac{4{\eta }_P}{\mathrm{i}q}}K\left(-\tilde{d},q\right)\\ {\displaystyle \frac{4{\eta }_P}{\mathrm{i}q}}K\left(\tilde{d},q\right)+{\displaystyle \frac{4{\eta }_P}{\mathrm{i}q}}K\left(0,q\right)\\ 0\\ 0\end{array}\right\}$$

(30)

where

$$\begin{array}{l}{A}_{11}=\pi -2\left({\eta }_M{q}^2- {\kappa }_v-2{\eta }_v\mathrm{i}q\right)K\left(0,q\right)\\ A_{12}=-2\left({\eta }_M{q}^2- {\kappa }_v-2{\eta }_v\mathrm{i}q\right)K\left(-\tilde{d},q\right)\\ A_{13}=-2\left({\kappa }_v+2{\eta }_v\mathrm{i}q\right)K\left(0,q\right)-2\left({\kappa }_v+2{\eta }_v\mathrm{i}q\right)K\left(-\tilde{d},q\right)\\ A_{14}=\left({\kappa }_v+2{\eta }_v\mathrm{i}q\right)K\left(0,q\right)\tilde{d}- \left({\kappa }_v+2{\eta }_v\mathrm{i}q\right)K\left(-\tilde{d},q\right)\tilde{d}\end{array}$$

(31)

$$\begin{array}{l}{A}_{21}=-2\left({\eta }_M{q}^2- {\kappa }_v-2{\eta }_v\mathrm{i}q\right)K\left(\tilde{d},q\right)\\ A_{22}=\pi -2\left({\eta }_M{q}^2- {\kappa }_v-2{\eta }_v\mathrm{i}q\right)K\left(0,q\right)\\ A_{23}=-2\left({\kappa }_v+2{\eta }_v\mathrm{i}q\right)K\left(\tilde{d},q\right)- 2\left({\kappa }_v+2{\eta }_v\mathrm{i}q\right)K\left(0,q\right)\\ A_{24}=\left({\kappa }_v+2{\eta }_v\mathrm{i}q\right)K\left(\tilde{d},q\right)\tilde{d}-\left({\kappa }_v+2{\eta }_v\mathrm{i}q\right)K\left(0,q\right)\tilde{d}\end{array}$$

(32)

$$\begin{array}{l}{A}_{31}= {\kappa }_v +2{\eta }_v\mathrm{i}q\\ A_{32}= {\kappa }_v +2{\eta }_v\mathrm{i}q\\ A_{33}={\mu }_v{q}^2-2{\kappa }_v -4{\eta }_v\mathrm{i}q\\ A_{34}=0\end{array}$$

(33)

$$\begin{array}{l}{A}_{41}={\kappa }_v\tilde{d}+2{\eta }_v\mathrm{i}q\tilde{d}\\ A_{42}=-{\kappa }_v\tilde{d}-2{\eta }_v\mathrm{i}q\tilde{d}\\ A_{43}=0\\ A_{44}=-2{\tilde{J}}_b{q}^2+{\kappa }_v{\tilde{d}}^2+2{\eta }_v\mathrm{i}q{\tilde{d}}^2\end{array}$$

(34)

After solving system (30), it is possible to substitute back into Equation (28) and perform the inverse Laplace transform. It is convenient to switch coordinates as shown below and, as already mentioned, use the method of contour integration so that the final result can be expressed as a sum of residues:

$$\tilde{w}\left(\xi ,\tau \right)={\displaystyle \frac{1}{2\pi \mathrm{i}}}\underset{T\to \infty }{\lim }{\displaystyle \underset{a-\mathrm{i}T}{\overset{a+\mathrm{i}T}{\int }}\tilde{W}\left(\xi ,\overline{q}\right){\mathrm{e}}^{\overline{q}\tau }\mathrm{d}\overline{q}}={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{-\mathrm{i}a-\infty }{\overset{-\mathrm{i}a+\infty }{\int }}\tilde{W}\left(\xi ,q\right){\mathrm{e}}^{\mathrm{i}q\tau }\mathrm{d}q}={\displaystyle \sum _j\mathrm{res}\left(\mathrm{i}\tilde{W}\left(\xi ,q\right){\mathrm{e}}^{\mathrm{i}q\tau },q_{M_j}\right)}+I_{bc}$$

(35)

where $a$ is a small real number, introduced so that all poles and discontinuities lie on the right-hand side of the integration line in the complex $q$-plane or after switching, above the integration line. $q_{M_j}$ are the induced frequencies and $I_{bc}$ represents the value obtained by numerical integration along the branch cuts. It has been shown that there are step discontinuities in the $K$-function, which must be removed by branch cuts. These discontinuities are along the line where the imaginary part of $q$ equals to ${\eta }_f$. More details in [23]. $I_{bc}$ must be calculated numerically, however its contribution is mostly negligible. Sometimes it is important to match the initial conditions exactly at the very initial times and thus can be neglected in most cases. Another situation is related to cases where some of the expected induced frequencies are missing.

To obtain the closed-form rear mass deflection (after neglecting $I_{bc}$), it is necessary to evaluate and sum up the residues as indicated in Equation (36)

$$\tilde{w}\left(0,\tau \right)={\displaystyle \sum _{j}res\left(\frac{4B_2\left[2B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)+\pi B_3\left(K_{-}+K_0\right)\right]{\mathrm{e}}^{\mathrm{i}q\tau }}{q{E}_{chr}},q_{M_j}\right)}$$

(36)

where:

$$\begin{array}{l}{B}_1= {\eta }_M{q}^2- {\kappa }_v-2\mathrm{i}q{\eta }_v; B_2={\mu }_v{q}^2-2{\kappa }_v -4\mathrm{i}q{\eta }_v; B_3=-2{\tilde{J}}_b{q}^2+{\kappa }_v{\tilde{d}}^2+2\mathrm{i}q{\eta }_v{\tilde{d}}^2;\\ B_4={\kappa }_v +2\mathrm{i}q{\eta }_v; K_0=K\left(0,q\right); K_{+}=K\left(\tilde{d},q\right); K_{-}=K\left(-\tilde{d},q\right);\\ B_5=K_{+}{K}_{-}-K_0^2; B_6=K_{+}+K_{-}+2K_0; B_7=K_{+}+K_{-}-2K_0;\end{array}$$

(37)

and

$$E_{chr}=4B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)\left(2B_4^2-B_1{B}_2\right)+\pi {\tilde{d}}^2{B}_4^2{B}_2{B}_7+2\pi B_4^2{B}_3{B}_6+\pi B_2{B}_3\left(\pi -4K_0{B}_1\right)$$

(38)

The denominator in the expression of Equation (36) corresponds to the simplified determinant of system (30). Its roots indicate the poles. There is always one obvious solution $q_M=0$ which identifies the steady-state part ${\tilde{w}}_{sd}\left(0\right)$ of the full solution, which is

$${\tilde{w}}_{sd}\left(0,\tau \right)={\tilde{w}}_{sd}\left(0\right)={\left.{\displaystyle \frac{4B_2\left[2B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)+\pi B_3\left(K_{-}+K_0\right)\right]}{E_{chr}}}\right|}_{q=0}$$

(39)

with $q=0$ substituted.

Then, if $q_M$ is a root of $E_{chr}=0$, then also the complex conjugate of the opposite value $\left(-q_M\right)*$ is a root. These pairs form harmonic vibration, named as unsteady harmonic part of the full solution ${\tilde{w}}_{usd}\left(0,\tau \right)$. The adjective “unsteady” is used to indicate that this part belongs to the transient solution. $I_{bc}$ is then called the truly transient part to distinguish it. The transient solution is therefore the sum of ${\tilde{w}}_{usd}\left(0,\tau \right)$ and $I_{bc}$.

$${\tilde{w}}_{usd}\left(0,\tau \right)={{\displaystyle \sum _j\left.\frac{4B_2\left[2B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)+\pi B_3\left(K_{-}+K_0\right)\right]{\mathrm{e}}^{\mathrm{i}q\tau }}{q{E}_{chr ,q}}\right|}}_{q=q_{M_j}}$$

(40)

As shown in Equation (40), it is necessary to evaluate the derivative of $E_{chr}$ with respect to $q$. This naturally implies derivatives of $K_0$, $K_{+}$ and $K_{-}$. Each of these values (without derivative) is also evaluated by contour integration. When using the derivative, all poles are doubled, but there is also an analytical formula for evaluating residues at double poles. Thus, the full evaluation is in fact analytical.

Similar formula as for ${\tilde{w}}_{usd}\left(0,\tau \right)$ can be deduced for the front mass deflection

$${\tilde{w}}_{usd}\left(\tilde{d},\tau \right)={\displaystyle \sum _{j}res\left(\frac{4B_2\left[2B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)+\pi B_3\left(K_{+}+K_0\right)\right]{\mathrm{e}}^{\mathrm{i}q\tau }}{q{E}_{chr}},q_{M_j}\right)}$$

(41)

for the mass center of the bogie bar vertical displacement

$${\tilde{u}}_{usd}\left(0,\tau \right)={\displaystyle \sum _{j}res\left(\frac{-4B_4\left[4B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)+\pi B_3{B}_6\right]{\mathrm{e}}^{\mathrm{i}q\tau }}{q{E}_{chr}},q_{M_j}\right)}$$

(42)

and for the rotation

$${\tilde{\theta }}_{usd}\left(0,\tau \right)={\displaystyle \sum _{j}res\left(\frac{4\pi \tilde{d}{B}_2{B}_4\left(K_{+}-K_{-}\right){\mathrm{e}}^{\mathrm{i}q\tau }}{q{E}_{chr}},q_{M_j}\right)}$$

(43)

The way $q$ is introduced implicates that the complete solution becomes unstable as soon as the imaginary part of $q_M$ becomes negative. $I_{bc}$ can never destabilize the solution because the branch cuts are located in the part of the complex $q$-plane with non-negative imaginary part.

After defining all the displacements at the active points, it is possible to substitute back into Equation (28) and determine the full deflection shapes of the beam at any time by closed form formula.

$${\tilde{w}}_{sd}\left(\xi ,\tau \right)={\tilde{w}}_{sd}\left(\xi \right)={\left.{\displaystyle \frac{8{\eta }_P{E}_{chr}+E_1}{2\pi E_{chr}}}\right|}_{q=0}K\left(\xi ,0\right)+{\left.{\displaystyle \frac{8{\eta }_P{E}_{chr}+E_2}{2\pi E_{chr}}}\right|}_{q=0}K\left(\xi -\tilde{d},0\right)$$

(44)

where

$$E_1=4B_1{B}_2\left(8B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)+4\pi B_3\left(K_{-}+K_0\right)\right)-2B_3\left[8B_3\left(4B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)+\pi B_3{B}_6\right)-4\pi {\tilde{d}}^2{B}_2{B}_4\left(K_{-}-K_{+}\right)\right]$$

(45)

$$E_2=4B_1{B}_2\left(8B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)+4\pi B_3\left(K_{+}+K_0\right)\right)-2B_3\left[8B_3\left(4B_5\left({\tilde{d}}^2{B}_4^2+B_1{B}_3\right)+\pi B_3{B}_6\right)+4\pi {\tilde{d}}^2{B}_2{B}_4\left(K_{-}-K_{+}\right)\right]$$

(46)

and

$${\tilde{w}}_{usd}\left(\xi ,\tau \right)={\displaystyle \sum _j\left[{\left.\frac{8{\eta }_P{E}_{chr}+E_1}{2\pi q{E}_{chr ,q}}\right|}_{q=q_{M_j}}K\left(\xi ,q_{M_j}\right)+{\left.\frac{8{\eta }_P{E}_{chr}+E_2}{2\pi q{E}_{chr ,q}}\right|}_{q=q_{M_j}}K\left(\xi -\tilde{d},q_{M_j}\right)\right]}{\mathrm{e}}^{\mathrm{i}{q}_{M_j}\tau }$$

(47)

3. Results

Several cases are selected to demonstrate in which cases instability can occur in the subcritical velocity range. However, this is not the only objective. It is also shown how the obtained results are related to the case of two moving masses, how they are affected by the input parameters and how they are affected by considering a Timoshenko-Rayleigh beam instead of Euler-Bernoulli. Several other cases are shown to demonstrate the ease of using the closed-form formulas. The results are validated by analysis on long finite beams, where the solution is obtained by modal expansion. Finally, the results are compared with other published works.

3.1. Induced Frequencies

In order to use the closed-form formulas deduced in Section 2, it is necessary to solve the complex roots of the complex characteristic equation

$$E_{chr}=0$$

(48)

These roots are called induced frequencies. Following the author’s previously published works, this can be done quite simply by implementing the iterative procedures described there. In general, there are four pairs of roots. As was explained above, it is enough for one of them to be unstable for the entire solution to be unstable.

3.2. Validation on Finite Beams

Similar to the previous works, the validation is performed on long finite beams by modal expansion. Since the number of modes must be significant, specifically 500 modes were implemented in all tested cases, it is necessary to reorder the terms in governing matrices to reduce the computational time. It has also been found that 500 modes are sufficient in all cases tested.

To solve the problem, the governing equations are not converted to moving coordinates and can be solved in original parameters and converted to dimensionless only for comparison of results. As in previous works, without loss of generality, the beam can be considered simply supported, since then the vibration modes are described by a sine function and the modal frequencies and numbers are analytical. Also, in accordance with the previous works, the loading must start slightly from the left support to eliminate its influence.

The beam equilibrium equations (1) and (7) are multiplied by a mode shape and integrated over the beam length. In the other equations (8) and (9), only the modal expansion of the beam deflection is introduced and then all the equations are grouped into matrix form

$$M\left(t\right)\cdot q_{,tt}\left(t\right)+C\left(t\right)\cdot q_{,t}\left(t\right)+K\left(t\right)\cdot q\left(t\right)=p\left(t\right)$$

(49)

The trick is to introduce additional variables

$$\begin{array}{l}{Q}_{1,tt}=P-{\displaystyle \sum _{k}M{w}_k\left(vt\right)q_{k,tt}}-{\displaystyle \sum _k\left(2Mv{w}_{k,x}\left(vt\right)+c_v{w}_k\left(vt\right)\right)}{q}_{k,t}\\ -{\displaystyle \sum _k\left(M{v}^2{w}_{k,xx}\left(vt\right)+c_{v}v{w}_{k,x}\left(vt\right)+k_v{w}_k\left(vt\right)\right)}{q}_k\end{array}$$

(50)

$$\begin{array}{l}{Q}_{2,tt}=P-{\displaystyle \sum _{k}M{w}_k\left(vt+d\right)q_{k,tt}}-{\displaystyle \sum _k\left(2Mv{w}_{k,x}\left(vt+d\right)+c_v{w}_k\left(vt+d\right)\right)}{q}_{k,t}\\ -{\displaystyle \sum _k\left(M{v}^2{w}_{k,xx}\left(vt+d\right)+c_{v}v{w}_{k,x}\left(vt+d\right)+k_v{w}_k\left(vt+d\right)\right)}{q}_k\end{array}$$

(51)

Then for the vector of unknowns $q\left(t\right)={\left[q_1,\dots q_n,u_v,\theta ,Q_1,Q_2\right]}^T$ the loading term is $p\left(t\right)={\left[0,\mathrm{...0},0,0,P,P\right]}^T$ and mass $M\left(t\right)$, damping $C\left(t\right)$ and stiffness $K\left(t\right)$ matrices have the following form

$$M\left(t\right)=\left[\begin{array}{ccccc}{\delta }_{jk}& 0& 0& -w_j\left(vt\right)& -w_j\left(vt+d\right)\\ 0& M_b& 0& 0& 0\\ 0& 0& J_b& 0& 0\\ M{w}_k\left(vt\right)& 0& 0& 1& 0\\ M{w}_k\left(vt+d\right)& 0& 0& 0& 1\end{array}\right]$$

(52)

$$C\left(t\right)=\left[\begin{array}{ccccc}{\delta }_{jk}{\displaystyle \frac{c_b}{m}}& -c_v{w}_j\left(vt\right)-c_v{w}_j\left(vt+d\right)& {\displaystyle \frac{1}{2}}d{c}_v{w}_j\left(vt\right)-{\displaystyle \frac{1}{2}}d{c}_v{w}_j\left(vt+d\right)& 0& 0\\ -c_v{w}_k\left(vt\right)-c_v{w}_k\left(vt+d\right)& 2c_v& 0& 0& 0\\ {\displaystyle \frac{d}{2}}\left(c_v{w}_k\left(vt\right)-c_v{w}_k\left(vt+d\right)\right)& 0& {\displaystyle \frac{1}{2}}{c}_v{d}^2& 0& 0\\ 2Mv{w}_{k,x}\left(vt\right)+c_v{w}_k\left(vt\right)& 0& 0& 0& 0\\ 2Mv{w}_{k,x}\left(vt+d\right)+c_v{w}_k\left(vt+d\right)& 0& 0& 0& 0\end{array}\right]$$

(53)

$$K\left(t\right)=\left[\begin{array}{ccccc}{\delta }_{jk}{\omega }_j^2& -k_v{w}_j\left(vt\right)-k_v{w}_j\left(vt+d\right)& {\displaystyle \frac{d}{2}}{k}_v{w}_j\left(vt\right)-{\displaystyle \frac{d}{2}}{k}_v{w}_j\left(vt+d\right)& 0& 0\\ \begin{array}{l}-k_v\left(w_k\left(vt\right)+w_k\left(vt+d\right)\right)\\ -c_{v}v\left(w_{k,x}\left(vt\right)+w_{k,x}\left(vt+d\right)\right)\end{array}& 2k_v& 0& 0& 0\\ \begin{array}{l}{\displaystyle \frac{d}{2}}{k}_v\left(w_k\left(vt\right)-w_k\left(vt+d\right)\right)\\ +{\displaystyle \frac{d}{2}}{c}_{v}v\left(w_{k,x}\left(vt\right)-w_{k,x}\left(vt+d\right)\right)\end{array}& 0& {\displaystyle \frac{1}{2}}{k}_v{d}^2& 0& 0\\ M{v}^2{w}_{k,xx}\left(vt\right)+k_v{w}_k\left(vt\right)+c_{v}v{w}_{k,x}\left(vt\right)& 0& 0& 0& 0\\ \begin{array}{l}+M{v}^2{w}_{k,xx}\left(vt+d\right)+k_v{w}_k\left(vt+d\right)\\ +c_{v}v{w}_{k,x}\left(vt+d\right)\end{array}& 0& 0& 0& 0\end{array}\right]$$

(54)

3.3. Timoshenko-Rayleigh Beam

When Timoshenko beam is considered, then Equation (1) is replaced by a set of equations

$$-G\overline{A}{w}_{,xx}+G\overline{A}{\phi }_{,x}+m{w}_{,tt}+c_f{w}_{,t}+k_{f}w= p\left(x,t\right)$$

(55)

$$-G\overline{A}{w}_{,x}+G\overline{A}\phi +J{\phi }_{,tt}-EI{\phi }_{,xx}= 0$$

(56)

where $G$ is the shear modulus, $\overline{A}$ is the reduced cross-sectional area, $\phi$ is the angle of rotation of the cross-section and $J$ is the mass moment of inertia of the cross-section which can be replaced by $m{r}^2$ with $r$ being the radius of gyration.

By switching to moving coordinates, one gets

$$-G\overline{A}{w}_{,rr}+G\overline{A}{\phi }_{,r}+m\left(w_{,tt}-2v{w}_{,rt}+v^2{w}_{,rr}\right)+c_f\left(w_{,t}-v{w}_{,r}\right)+k_{f}w=p\left(r,t\right)$$

(57)

$$-G\overline{A}{w}_{,r}+G\overline{A}\phi +m{r}^2\left({\phi }_{,tt}-2v{\phi }_{,rt}+v^2{\phi }_{,rr}\right)-EI{\phi }_{,rr}= 0$$

(58)

After introduction of additional dimensionless parameters

$$\tilde{\phi }={\displaystyle \frac{\phi }{w_{st}\chi }}, \tilde{S}={\displaystyle \frac{G\overline{A}}{2\sqrt{k_{f}EI}}}, \tilde{r}=r\chi$$

(59)

this leads into

$$-4\tilde{S}{\tilde{w}}_{,\xi \xi }+4\tilde{S}{\tilde{\phi }}_{,\xi }+4{\alpha }^2{\tilde{w}}_{,\xi \xi }+4{\tilde{w}}_{,\tau \tau }-8\alpha {\tilde{w}}_{,\xi \tau }+8{\eta }_f\left({\tilde{w}}_{,\tau }-\alpha {\tilde{w}}_{,\xi }\right)+4\tilde{w}=\tilde{p}\left(\xi ,\tau \right)$$

(60)

$$-4\tilde{S}{\tilde{w}}_{,\xi }+4\tilde{S}\tilde{\phi }+4{\tilde{r}}^2\left({\alpha }^2{\tilde{\phi }}_{,\xi \xi }-2\alpha {\tilde{\phi }}_{,\xi \tau }+{\tilde{\phi }}_{,\tau \tau }\right)-{\tilde{\phi }}_{,\xi \xi }= 0$$

(61)

After application of integral transforms, the image of $\tilde{\phi }$ can be solved

$$\tilde{\Phi }= {\displaystyle \frac{4\mathrm{i}p\tilde{S}}{4\tilde{S}-4{\tilde{r}}^2{\left(\alpha p-q\right)}^2+p^2}}\tilde{W}$$

(62)

and introduced into the first equation. Further steps follow the previous derivations.

3.4. Tested Cases

Several cases were analyzed to fulfill the objectives of this paper. In order to define some admissible range of parameters, it is necessary to consider a wide range of possible real parameters. As for the supporting structure, the values of $EI$ and $m$ are contained in a rather narrow interval, approximately as $EI\in 〈4.7;6.4〉{\mathrm{MNm}}^2$ and $m\in 〈54;60〉\mathrm{kg}/\mathrm{m}$, respectively. $k_s$ is sometimes omitted in similar models, so that the value of zero is allowed and the major variety is assigned to $k_f\in 〈0.22;1000〉{\mathrm{MN}/\mathrm{m}}^2$. Regarding the moving bogie, it would not be appropriate to mix all possible values, since the mass of the bogie bar should be related to the wheel masses and, together with the distance between the wheels, should determine approximately its mass moment of inertia. For this reason, two additional coefficients, $c{f}_1$ and $c{f}_2$, are introduced to link these values as

$${\mu }_v=c{f}_1{\eta }_M, {\tilde{J}}_b=c{f}_2{\eta }_M{\tilde{d}}^2$$

(63)

In fact, using Equations (16) and (17), this also means that

$$M_b=c{f}_{1}M, J_b=c{f}_{2}M{d}^2$$

(64)

There is a wide range of values for primary suspension stiffness, thus in the following examples two values have been used: 0.02 for soft and 0.2 for stiffer case. Damping parameters ${\eta }_f$, ${\eta }_v$ and velocity ratio $\alpha$ vary according to the effects being analyzed. The input data of all tested cases are listed in Table 1. Additionally, Table 2 summarizes the induced frequencies as they were used for residues evaluation.

The general tendency of the frequency lines is shown in Figure 2 using the parameters of Case 1, but for the entire range of velocities. Since the values are quite far apart, they are plotted separately. In general, there are four pairs of induced frequencies. Their evolution is interrupted when their imaginary part reaches ${\eta }_f$, as was demonstrated in [23].

Contrary to other more complicated models, here the critical velocity for resonance, meaning the critical velocity of one constant force traversing at constant velocity the one-layer model is well-defined by analytical formula.

$${\alpha }_{cr}^{E-B}={\alpha }_{ref}\sqrt{1+{\eta }_s}\mathrm{with} {\alpha }_{ref}=1$$

(65)

For the Timoshenko-Rayleigh beam the critical velocity is presented for instance in [51] or [48], with dimensionless parameters used here it is more explicit how this value relates to Euler-Bernoulli beam

$${\alpha }_{cr}^{T-R}={\alpha }_{ref}{\displaystyle \frac{\sqrt{{\tilde{r}}^2-\tilde{S}-8{\tilde{r}}^2{\tilde{S}}^2+4{\tilde{S}}^{3/2}\sqrt{4{\tilde{r}}^4\tilde{S}-{\tilde{r}}^2+\tilde{S}}}}{2\left(\tilde{S}-{\tilde{r}}^2\right)}}$$

(66)

3.5. Comparison with Moving Masses

Unlike a single moving inertial object, i.e. an object with a single contact point with the supporting structure, where instability in the subcritical velocity range is mathematically impossible, when there are more moving objects, their dynamic interaction does not respect this condition. In fact, instability in the subcritical velocity range occurs very often, but for objects with unrealistically high mass. This has already been shown in [39]. Here, the negative effect of foundation damping is summarized by the representation of the instability lines in Figure 3. To demonstrate the features being analyzed and for perception of general tendencies, the vertical axis is extended up to ${\eta }_M=200$, which is already only an academic value. $\alpha$ expands up to 1.2 and $\alpha =1$ is the critical velocity for resonance. For simplicity and further correspondence with moving bogie, two equal moving masses acted by the same constant forces are considered. Then, in dimensionless parameters, when the effect of shear stiffness is neglected, the only parameters to vary are the distance between mases and foundation damping.

Unlike in all other works published by other researchers (e.g. [48,49,50]) where the instability is determined by the D-decomposition method, here the so-called lines of instability are traced in the plane of velocity-moving mass ratios by identifying the real-valued induced frequencies, which in a damped case obligatorily alter by one the number of unstable induced frequencies, so that stable and unstable velocity intervals are clearly visible for any specific moving mass ratio. This approach is much more powerful and simpler than the D-decomposition method.

Figure 3. Instability lines of two moving proximate masses of the same value. Different colors stand for different distances $\tilde{d}=1:0.25:2.25$ as indicated in the legend. Black dotted line marks the critical velocity. Parts (a) to (f) correspond to ${\eta }_f=0.05:0.05:0.3$.

Figure 3. Instability lines of two moving proximate masses of the same value. Different colors stand for different distances $\tilde{d}=1:0.25:2.25$ as indicated in the legend. Black dotted line marks the critical velocity. Parts (a) to (f) correspond to ${\eta }_f=0.05:0.05:0.3$.

It can be concluded that with increasing damping, all branches of instability shift to lower moving masses and at the same time become straighter, so that lower masses occupy larger velocity intervals. It has thus been proven that contrary to what is generally expected, foundation damping worsens the situation. However, it is necessary to note that in closed intervals the degree of instability is usually much lower than in semi-closed intervals.

In Figure 3 only relevant distances between mases are included, that is those that still have part of the branch within the defined scale in subcritical velocity range. By analyzing a larger range of distances between masses, it can be concluded that the starting point of the first branch of instability shifts to the left with increasing distance, but at the same time the moving masses increase.

When the moving bogie is considered, there are more parameters that can vary. To make the comparison more realistic, the total mass should be the same. The case with $c{f}_1=1.7$ is considered in Figure 4 for comparison and thus the previous instability lines from Figure 3 are recalculated accordingly. Besides $c{f}_1=1.7$, the moving bogie is characterized by $c{f}_2=0.01$ and ${\kappa }_v=0.02$ to make its influence as less as possible.

It can be concluded that even for weak primary suspension stiffness and low mass moment of inertia of the bogie bar, all cases with moving bogie perform better, meaning the instability lines are above the ones for moving masses, except for the lowest tested distance $\tilde{d}=1$. Instability lines in Figure 4 are represented only up to the critical velocity. The relevant instability lines are those that have a part of a branch within the chosen scale in the subcritical velocity range. Then the following rule applies: if $\tilde{d}$ is relevant for moving masses but not for a moving bogie, then $\tilde{d}+0.25$ is not plotted for moving masses, even if it is relevant.

To analyze in detail the only case of moving bogie that performs worse than moving masses, Case 1 from Table 1 is selected. This case is unstable. Nevertheless, as it is seen from Table 2, the unstable frequency has quite a low imaginary part (in absolute value). This means that the unstable frequency is gaining importance only in higher times, while in the beginning this case resembles a highly damped case. In Figure 5 time series are shown, obtained by closed-form formulas and also on finite beam for validation. The high times in Figure 6 are calculated using closed-form formulas only, it would be computationally untenable to do this on a finite beam due to the excessive number of vibration modes and also using numerical integration due to numerical precision.

3.6. Effect of the Suspension Damping

In this section it will be shown that the damping of the suspension has a very positive effect on the bogie stability, moreover this effect is helped by the higher mass moment of inertia of the bogie bar. To demonstrate this, $c{f}_1=1.7$ and the other coefficient varies between $c{f}_2=0.01$ and $c{f}_2=0.1$. Weak foundation damping of ${\eta }_f=0.05$ and stiffer primary suspension stiffness of ${\kappa }_v=0.2$ are implemented.

Figure 7 clearly shows that with a higher mass moment of inertia of the bogie bar, practically all instability lines of the cases with primary suspension damping are removed from the subcritical velocity range.

3.6. Effect of the Beam Theory

As mentioned earlier, it is also necessary to see how the results are affected by the adopted beam theory. Two cases are shown in Figure 8. The additional parameters of the Timoshenko-Rayleigh beam theory should not be chosen independently, as is clear from Equation (59). Therefore, $\tilde{S}=100$, $\tilde{r}=0.03$ are selected for a softer foundation and $\tilde{S}=10$, $\tilde{r}=0.09$ for a stronger one. The other parameters are: $c{f}_1=1.7$, $c{f}_2=0.1$. Weak foundation damping of ${\eta }_f=0.05$ and a stiffer primary suspension stiffness of ${\kappa }_v=0.2$ with no damping are implemented.

Only the subcritical velocity range is shown, However, since the difference in critical velocity in the Timoshenko-Rayleigh case using Equation (66) is very small, 0.993 and 0.980, respectively, the graphs are plotted up to the value 1, as before.

Figure 8. Instability lines for a moving bogie. Different colors stand for different distances $\tilde{d}=1:0.25:2.25$ as indicated in the legend. Full curves are for the Euler-Bernoulli theory and dashed curves for the Timoshenko-Rayleigh. (a) $\tilde{S}=100$, $\tilde{r}=0.03$; (b) $\tilde{S}=10$, $\tilde{r}=0.09$.

Figure 8. Instability lines for a moving bogie. Different colors stand for different distances $\tilde{d}=1:0.25:2.25$ as indicated in the legend. Full curves are for the Euler-Bernoulli theory and dashed curves for the Timoshenko-Rayleigh. (a) $\tilde{S}=100$, $\tilde{r}=0.03$; (b) $\tilde{S}=10$, $\tilde{r}=0.09$.

It can be concluded that for softer foundation there are practically no differences between results obtained using Euler-Bernoulli and Timoshenko-Rayleigh theories. Some differences are noticeable for a stronger foundation, but they do not seem to be fundamental.

3.7. Time Series of Other Cases from Table 1

The objectives of this subsection are threefold, since in addition to the analyzed effect, it is shown how easy it is to use the closed-form formulas, and it is also verified by analysis on finite beams that the contribution of integration along the branch cuts is negligible.

Cases 2 and 3 are selected to show that the dynamic interaction discussed in previous works for two moving masses is maintained for moving bogie. Slightly higher $c{f}_1=2$ is introduced, but the changes are not significant. Two distances are compared, $\tilde{d}=2$ (Figure 9) and $\tilde{d}=5$ (Figure 10). Other parameters are included in Table 1.

Figure 9 shows a very strong dynamic interaction. In Figure 10 this effect is not as pronounced, however, like for two moving masses, the dynamic interaction persists up to relatively high distances. Only when $K_{+}$ and $K_{-}$ start to have negligible values with respect to $K_0$, then the interaction is negligible. This proves that the function envelopes for $\tilde{\theta }$ still show very strong interaction in Figure 10, because $K_0$ is not included in Equation (43).

Case 4 is selected for comparison with Case 3, to see the effect of $c{f}_1=1.7$ with respect to $c{f}_1=2$. Besides the time series (Figure 11) it is also interesting to see differences in induced frequencies in Table 2.

Cases 2-4 are undamped; therefore, all amplitudes remain unchanged, the induced frequencies are real (Table 2), and thus these cases also serve to demonstrate how easy it is to deal with undamped cases without numerical problems.

Cases 5 and 6 are chosen such that Case 5 is unstable for subcritical velocity, but by adding damping to the primary suspension and leaving the other parameters unchanged, it becomes stable – Case 6. Time series of Case 5 are reported in Figure 12. During the initial stages there is no indication of instability as the rate of instability is very low (Table 2). The unstable frequency is overpassing the stable ones only at higher times (Figure 13).

Case 6 is already stable. Times series are shown in Figure 14 for initial stages and in Figure 15 for higher times.

The other cases from Table 1 are summarized in Appendix A.

3.8. Comparison with Other Published Works

While in [49] only isolated cases are shown, in [48] an evolution over a range of velocities is presented, which provides a better insight into the problem. Nevertheless, the reference case corresponds to a case with negligible dynamic interaction between the wheels, which is assumed as 15m, so it is not an actual bogie but a simplified vehicle with two axels. The evolution of the critical stiffness of the primary suspension is chosen for comparison. The results are shown in Figure 16 using the input data from [48].

It can be concluded that there is a very good agreement after recalculation to dimensionless data as in [48]. Some minor differences can be justified by the fact that nowadays it is possible to use symbolic software allowing much higher accuracy and faster evaluation than at the time of [48] was written.

In [48,49,50] the distance between wheels is so large that the dynamic interaction is not significant, this way also the instability in subcritical velocity range does not occur.

4. Discussion

This paper presents a detailed analysis of the instability of a bogie passing through a one-layer model, focusing on the possibility that this may occur in the subcritical velocity range. All presented results were obtained (semi)analytically, quickly and with very high accuracy.

The results of two moving proximate masses are used as a reference case and it is shown that in general, the more realistic the moving object, the lower the risk of instability at such low velocities. This was confirmed even under the conditions of neglected suspension damping, reduced suspension stiffness and small mass moment of inertia of the bogie bar. On the other hand, it was confirmed that the negative effect of foundation damping found in the two moving mass cases is preserved for the moving bogie. This effect is in contrast to the effect of foundation damping in the case of a single moving inertial object, where increased damping always shifts the onset of instability to higher velocities.

It has been demonstrated that instability lines related to moving bogies have more branches than those related to moving masses, which complicates the analysis.

It has been shown that suspension damping improves the unstable behavior in all cases, but it is much more efficient with a higher mass moment of inertia of the bogie bar.

Regarding the distances between the masses: it is not true that the lowest means worse. Larger distances shift the starting point of the first instability branch to lower speeds, but at the same time the smallest moving mass ratio of such branch increases, therefore, it is necessary to balance these two tendencies.

It has been proven that there are cases with realistic parameters where instability occurs in the subcritical velocity range, and therefore this article serves as a warning that instability at supercritical velocities cannot be taken for granted a priori.

Unlike in all other works published by other researchers (e.g. [48,49,50]) where the instability is determined by the D-decomposition method, here the so-called lines of instability are traced in the plane of velocity-moving mass ratios by identifying the real-valued induced frequencies, which in a damped case obligatorily alter by one the number of unstable induced frequencies, so that stable and unstable velocity intervals are clearly visible for any specific moving mass ratio. This approach is much more powerful and simpler than the D-decomposition method.

The presented analysis can be easily extended to other types of supporting structures, since all expressions can be used for more complicated foundation models, simply by replacing Equation (24) with an equation reflecting the new foundation model. There is only one complication, namely that in other foundation models the analytical formulas have limited use, since discontinuities in the K-function do not allow to resolve all the necessary induced frequencies for all velocities. This is only relevant for the time series; the instability lines are not affected. As an alternative for the time series, either numerical integration or finite models can be used.

5. Conclusions

In this paper the possibility of instability in subcritical range of velocities was discussed in detail. The results presented are very easy to use, as the dimensionless parameters cover a wide range of possible scenarios, and the instability lines provide an overall view of the onset or aggravation of unstable behavior or the recovery of stable behavior. Such comprehensive information is very useful, as readers do not have to repeat the calculations and can simply calculate their dimensionless parameters for the case they are analyzing.

The validation of the results on finite volumes is clearly described and easily applicable to other researchers. This is simpler than the numerical integration associated with the inverse Laplace transform.

The unstable induced frequencies compared to the stable ones provide an initial estimate of the degree of instability. Indeed, it can be seen in Table 2 that the unstable frequencies of the selected cases have a relatively low imaginary part (in absolute value). As has been shown, the effect of foundation damping worsens the instability, however, during the initial times it has a very strong impact on the stable frequencies, so that the vibrations initially resemble the heavily damped cases. The unstable frequency only surpasses the stable ones at higher times. This is important for assessing the severity of instability and can be used in vibration mitigation measures. Since the analysis allows the determination of the vibration amplitude, together with the degree of instability, it identifies the time required for the unstable harmonics to overcome the stable ones. For this reason, it is important to have analytical formulas that are not affected by numerical approximations during numerical integration and the time response can be extended as needed.