13162 FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 22, No 4, 2024, pp. 773 - 785 https://doi.org/10.22190/FUME241003045S © 2024 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper DYNAMIC STABILITY ANALYSIS OF A COUPLED MOVING BOGIE SYSTEM Vladimir Stojanović1,2, Lubomir Dimitrov3, Pancho Tomov3, Deli Li2, Vlastimir Nikolić1 1Univerity of Niš, Faculty of Mechanical Engineering, Serbia 2Lakehead University, Ontario, Canada 3Technical University of Sofia, Bulgaria ORCID iDs: Vladimir Stojanović https://orcid.org/0000-0003-1011-4979 Lubomir Dimitrov https://orcid.org/0000-0001-8584-5958 Pancho Tomov https://orcid.org/0000-0003-1649-6947 Deli Li https://orcid.org/0009-0005-5375-402X Vlastimir Nikolić https://orcid.org/0009-0004-5752-5083 Abstract. This study explores the dynamic behavior and vibration stability of a complex moving mechanical oscillator coupled with a three-part viscoelastically connected continuous beam-foundation system. Elastic waves generated by the oscillator can destabilize the system, a scenario common in high-speed trains and vehicles. The paper presents an improved analytical approach, focusing on the effects of a variable primary stiffness suspension on vibration stability. Using the D-decomposition method and the principle of the argument, instability intervals are determined for varying stiffness configurations, highlighting the benefits of nontraditional suspension systems in improving local stability. Key words: Instability, Critical velocity, D-decomposition method, Principle of argument 1. INTRODUCTION The rapid expansion of road and railway transport has intensified the need for a detailed analysis of the dynamic effects of physical parameters on the stability of mechanical systems. Vibration instability, often resulting from energy transfer through anomalous Doppler waves [1, 2], is a particular concern when high-speed trains exceed the propagation speed of elastic waves within railway structures. This phenomenon has been extensively studied, particularly in simple mechanical models, leading to further investigations into more complex systems under moving loads. Michaltsos [3] analyzed the linear dynamic response of a simply supported beam subjected to a moving load with constant magnitude and variable velocity. The study examines the effects of acceleration and deceleration on the beam's behavior under both a single-axle load and a two-axle vehicle model, with damping considered for the latter. Subsequent studies, such as those by Kim [4] and Śniady [5], incorporated additional factors like rotary inertia and axial compressive forces into the Received: October 03, 2024 / Accepted November 06, 2024 Corresponding author: Vladimir Stojanović Faculty of Mechanical Engineering, University of Niš, A. Medvedeva 14, 18000 Niš, Serbia E-mail: vladimir.stojanovic@masfak.ni.ac.rs https://orcid.org/0000-0003-1011-4979 https://orcid.org/0000-0001-8584-5958 https://orcid.org/0000-0003-1649-6947 https://orcid.org/0009-0005-5375-402X https://orcid.org/0009-0004-5752-5083 774 V. STOJANOVIĆ, L. DIMITROV, P. TOMOV, D. LI, V. NIKOLIĆ stability analyses, enriching the understanding of these dynamic systems. Basu and Rao [6] derived analytical solutions for the steady-state response of an infinite beam on a viscoelastic foundation under a constant-velocity moving load. Similarly, Froio et al. [7] explored the steady-state response of a uniform, infinite Euler-Bernoulli beam on a Pasternak foundation subjected to a constant-velocity moving load. They obtained a universal closed-form solution using a Fourier transform, capable of representing the response for any set of beam and foundation parameters. Mazilu presented significant results in [8], addressing the interaction between a vehicle, simplified to a two-mass oscillator, and a slab track modeled as an infinite structure composed of elastically supported double Euler-Bernoulli beams. To carry out the time-domain analysis, a semi-analytical method was developed, leveraging the remarkable properties of the slab track's time-domain Green's functions. Meanwhile, Karimi and Ziaei- Rad [9] studied the nonlinear coupled vibration analysis of a beam with moving supports under the action of a moving mass. Researchers such as Verichev and Metrikine [10, 11, 12] successfully applied the D-decomposition method to define stability regions for systems with moving masses along beams. Similarly, Dimitrovová and Varandas [13] examined the transient dynamic response of a beam supported on a foundation with a sudden change in stiffness and subjected to a constant-velocity moving force, utilizing finite integral transformations. In a related study, Dimitrovová [14] proposed a repetitive procedure for determining oscillation frequencies induced by a moving mass. Additionally, Erbas et al. [15] performed a three-dimensional analysis of near-resonant regimes of a point load moving uniformly along the surface of a coated elastic half-space. Svedholm et al. [16] contributed valuable insights into the behavior of loads moving on beams of finite length, while Lee and Renshaw [17] introduced a new solution technique for solving the moving mass problem in nonconservative, linear, distributed parameter systems using complex eigenfunction expansions. All these works underscored the importance of advanced analytical methods, such as complex eigenfunction expansions and finite integral transformations, for understanding the interactions between continuous and discrete systems. These methods are critical for analyzing the dynamic challenges posed by high-speed trains, where the interaction between the moving oscillator (train) and the supporting structure (rail) produces complex responses. Stojanović et al. [18] conducted an in-depth analysis of the vibration stability of a coupled bogie system moving uniformly along a flexibly supported infinite high-order shear deformable coupled beam system on a viscoelastic base. This paper presents a model of a complex coupled system with four supports, featuring varying suspension configurations. The stability of vibrations is analyzed by examining the interaction between an infinite beam and a complex mechanical oscillator moving uniformly. Building upon the work of Stojanović and Petković [19], which explored the stochastic stability of complex finite beam and beam-foundation systems, this study extends the analysis to a high-order shear deformable model of an infinite beam system subjected to a moving oscillator. Stability regions for the complexly coupled mechanical oscillator are determined for different stiffness moduli of the primary suspension, providing a comprehensive analysis of vibration stability in relation to varying suspension configurations. 2. MODELING OF THE COMPLEX MOVING OSCILLATOR The model of the complex foundation system is illustrated in Fig. 1. In this configuration, the rail is represented as a supported infinite beam, while the slab and base are modeled as Dynamic Stability Analysis of a Coupled Moving Bogie System 775 additional beams with free boundary conditions. These elements are interconnected by viscoelastic spring-damping components, which are crucial for capturing the real behavior of the track structure. These viscoelastic connections account for both stiffness and damping, significantly influencing the system's dynamic response to moving loads. To accurately capture the behavior of the multilayer beam system, the high-order shear deformable beam theory (also known as the Reddy-Bickford beam theory) is employed. This advanced theory extends beyond the classical beam theory by incorporating the effects of shear deformation and rotational inertia. The Reddy-Bickford theory allows for a more realistic and detailed representation of the coupled mechanical interactions between the beam layers in the system. The interaction between the infinite coupled beam-layer system is described by the following set of equations. Ψ1+kf3(w1-w2)+νf3(ẇ1-ẇ2)=p(x,t), (1) Φ1=0, (2) Ψ2+kf2(w2-w3)-kf3(w1-w2)+νf2(ẇ2-ẇ3)-νf3(ẇ1-ẇ2)=0, (3) Φ2=0, (4) Ψ3+kf1w3-kf2(w2-w3)+νf1ẇ3-νf2(ẇ2-ẇ3)=0, (5) Φ3=0, (6) where Ψi=Ci1ẅi-Ci2(φi '+wi '')-Ci3ẅi ''+Ci4wi (4) +Ci5φ̈ i '+Ci6φ i ''', Φi=Di1φ̈ i -Di2φ i ''+Di3(φi +wi ')+Di4ẅi '+Di5wi '' ', i=1, 2, 3, and where Cik and Dij are constants obtained and documented in the paper by Stojanović and Petković [20] in the form for one beam (appropriate indices need to be added for each beam individually) and can be utilized for all others. E and G are the Young’s modulus and the shear modulus of the beam material; ρ is the mass density of the beam material; kfi and νfi are the stiffness moduli and the viscosity shear moduli of the foundation per unit length. wi(x, t) are vertical deflections of the beams and φ i (x,t) are the angles of rotation of the cross-section of the beams. The following notation was used to write the expression ẅi= ∂ 2 wi ∂t2 , wi '= ∂wi ∂x , wi ''= ∂ 2 wi ∂x2 , φ̈ i = ∂ 2 wi ∂t2 , φ i '= ∂wi ∂x , φ i ''= ∂ 2 wi ∂x2 . Right hand side of Eq. (1) represents the reaction of the bogie system and has the following form p(x,t)=-δ(x-vt-d3) (m d 2 w2 01 dt2 +L̃214) -δ(x-vt-d2) (m d 2 w2 02 dt2 +L̃223) -δ(x-vt-d1) (m d 2 w1 01 dt2 +L̃112) -δ(x-vt) (m d 2 w1 02 dt2 +L̃121) . (7) w1 01(t), w1 02(t), w2 01(t), w2 02(t), are the vertical displacements of the masses (wheels); w1 0(t), w2 0(t) are the vertical displacements of the centers of mass of the bars that represent the bogies’ bodies, and w0 0(t) is the vertical displacement of the bar that represents the car’s 776 V. STOJANOVIĆ, L. DIMITROV, P. TOMOV, D. LI, V. NIKOLIĆ body. θ0, θ1, and 𝜃2 are the angles of rotation of the bars around their centers of mass. M0, M1 and J0, J1 are the masses and moments of inertia of the bars; k0 and ϵ0 are the stiffness modulus and the damping factor of the primary suspension; k1 and ϵ1 are the stiffness modulus and the damping factor of the secondary suspension. The Dirac delta- function is designated with δ (∙). The factors ζ1, ζ2, ζ3, and ζ4 are multipliers for the stiffness modulus of the primary suspension. Fig. 1 Complexly coupled mechanical oscillator Together with Eqs. (1-3) the following set of Eqs. (8-13) represent the full form of the dynamic equilibrium of complex moving oscillator M1 d 2 w2 0 dt2 -L̃214-L̃223+L̃2 =0, (8) (J1+M 1 r1 2) d 2 θ2 dt2 +L̃214 ζ4 ζ3+ζ4 dL-L̃223 ζ3 ζ3+ζ4 dL=0, (9) M1 d 2 w1 0 dt2 -L̃112-L̃121+L̃1=0, (10) (J1+M 1 r2 2) d 2 θ1 dt2 +L̃112 ζ2 ζ1+ζ2 dL-L̃121 ζ1 ζ1+ζ2 dL=0, (11) Dynamic Stability Analysis of a Coupled Moving Bogie System 777 M0 d 2 w0 0 dt2 -L̃1-L̃2=0, (12) J0 d 2 θ0 dt2 +(L̃2-L̃1) dL0 2 =0, (13) w2 01(t)=w(vt+d3,t), w2 02(t)=w(vt+d2,t), w1 01(t)=w(vt+d1,t), w1 02(t)=w(vt,t), lim |x-vt|→∞ wi(x,t)=0, lim |x-vt|→∞ φ i (x,t)=0, (14) where L̃ijp=ζpk 0 (wi 0j -wi 0+θilij(-1) j+1)+ϵ0 d dt (wi 0j -wi 0+θilij(-1) j+1), L̃1=k 1 (w1 0-w0 0-θ0 dL0 2 -r2θ1)+ϵ1 d dt (w1 0-w0 0-θ0 dL0 2 -r2θ1) , L̃2=k 1 (w2 0-w0 0+θ0 dL0 2 +r1θ2)+ϵ1 d dt (w2 0-w0 0+θ0 dL0 2 +r1θ2) , l11= ζ2 ζ1+ζ2 dL,l12= ζ1 ζ1+ζ2 dL,l21= ζ4 ζ3+ζ4 dL, l22= ζ3 ζ3+ζ4 dL, r1=l22- dL 2 , r2=l11- dL 2 . By using a moving reference system ξ=x-vt the previous equations can be rewritten. It is now possible to analyze the characteristic equation of vibrations by applying integral transforms to the modified system of equations Wsi(ξ,s)=∫ wi(ξ,t) e-stdt ∞ 0 , Φsi(ξ,s)=∫ φ i (ξ,t) e-stdt, ∞ 0 Ŵk,si(k,s)=∫ Wsi(ξ,s)e-ikξdξ ∞ -∞ , Φ̂si(k,s)=∫ Φsi(ξ,s)e-ikξdξ ∞ -∞ W1 01(s)=∫ w1 01(t)e-stdt ∞ 0 , W1 02(s)=∫ w1 02(t)e-stdt ∞ 0 , W2 01(s)=∫ w2 01(t)e-stdt ∞ 0 , W2 02(s)=∫ w2 02(t)e-stdt ∞ 0 , W0 0(s)=∫ w0 0(t)e-stdt, ∞ 0 W1 0(s)=∫ w1 0(t)e-stdt, ∞ 0 W2 0(s)=∫ w2 0(t)e-stdt, ∞ 0 Θ0=Θ 0 (s)=∫ θ0(t)e -stdt, ∞ 0 Θ1=Θ 1 (s)=∫ θ1(t)e -stdt, ∞ 0 Θ2=Θ 2 (s)=∫ θ2(t)e -stdt. ∞ 0 778 V. STOJANOVIĆ, L. DIMITROV, P. TOMOV, D. LI, V. NIKOLIĆ By using the procedure given in Stojanović et al. [18] we are presuming that the initial conditions of the problem Eqs. (1-3) and Eqs. (8-13) are insignificant (the initial conditions do not have an impact on the system stability because of the fact that the governing equations are linear) and with keeping in mind that: W1 02(s)= W(0,s), W1 01(s)= W(d1,s) ,W2 02(s)= W(d2,s) ,W2 01(s)= W(d3,s), the application of the transforms results in the system of algebraic equations after elimination Φ̂si(k,s) in the present form F(k,s) Ŵk,s(k,s)=e-id1k(-ζ2k0 W(d1,s) -ms2 W(d1,s) -sϵ0 W(d1,s) +ζ2k0W1 0(s)+sϵ0W1 0(s)- ζ2 2 Θ1k0dL ζ1+ζ2 - ζ2Θ1sϵ0dL ζ1+ζ2 ) +e-id2k(-ζ3k0 W(d2,s) -ms2 W(d2,s) -sϵ0 W(d2,s)+ζ3k0W2 0(s) +sϵ0W2 0(s)+ ζ3 2 Θ2k0dL ζ3+ζ4 + ζ3Θ2sϵ0dL ζ3+ζ4 ) +e-id3k(-ζ4k0W (d3,s)-ms2W(d3,s)-sϵ0W(d3,s) +ζ4k0W2 0(s)+sϵ0W2 0(s)- ζ4 2 Θ2k0dL ζ3+ζ4 - ζ4Θ2sϵ0dL ζ3+ζ4 ) -ms2W(0,s)-sϵ0W(0,s)+ζ1k0W1 0(s)+sϵ0W1 0(s)+ ζ1 2 Θ1k0dL ζ1+ζ2 + ζ1Θ1sϵ0dL ζ1+ζ2 , (15) M1s2W2 0(s)-ζ3k0 W(d2,s) -ζ4k0 W(d3,s) -sϵ0 W(d2,s) -sϵ0 W(d3,s) +ζ3k0W2 0(s)+ζ4k0W2 0(s)-k1W0 0(s)+k1W2 0(s)-sϵ1W0 0(s)+2sϵ0W2 0(s) +sϵ1W2 0(s)+Θ2k1 ( ζ3dL ζ3+ζ4 - dL 2 )+ ζ3 2 Θ2k0dL ζ3+ζ4 - ζ4 2 Θ2k0dL ζ3+ζ4 + 1 2 Θ0k1dL0+Θ2sϵ1 ( ζ3dL ζ3+ζ4 - dL 2 )+ ζ3Θ2sϵ0dL ζ3+ζ4 - ζ4Θ2sϵ0dL ζ3+ζ4 + 1 2 Θ0sϵ1dL0=0, (16) Θ2J1s2+Θ2M1s2 ( ζ3dL ζ3+ζ4 - dL 2 ) 2 - ζ3 2 k0dL W(d2,s) ζ3+ζ4 + ζ4 2 k0dL W(d3,s) ζ3+ζ4 - ζ3sϵ0dLW(d2,s) ζ3+ζ4 + ζ4sϵ0dLW(d3,s) ζ3+ζ4 + ζ3 2 k0dLW2 0(s) ζ3+ζ4 - ζ4 2 k0dLW2 0(s) ζ3+ζ4 + ζ3sϵ0dLW2 0(s) ζ3+ζ4 - ζ4sϵ0dLW2 0 (s) ζ3+ζ4 + ζ3 3 Θ2k0dL 2 (ζ3+ζ4) 2 + ζ4 3 Θ2k0dL 2 (ζ3+ζ4) 2 + ζ3 2 Θ2sϵ0dL 2 (ζ3+ζ4) 2 + ζ4 2 Θ2sϵ0dL 2 (ζ3+ζ4) 2 =0, (17) M1s2W1 0(s)-ζ2k0 W(d1,s) -sϵ0 W(d1,s) -ζ1k0 W(0,s) -sϵ0 W(0,s) +ζ1k0W1 0(s)+ζ2k0W1 0(s)-k1W0 0(s)+k1W1 0(s)-sϵ1W0 0(s)+2sϵ0W1 0(s) Dynamic Stability Analysis of a Coupled Moving Bogie System 779 +sϵ1W1 0(s)-Θ1k1 ( ζ2dL ζ1+ζ2 - dL 2 )+ ζ1 2 Θ1k0dL ζ1+ζ2 - ζ2 2 Θ1k0dL ζ1+ζ2 - 1 2 Θ0k1dL0 -Θ1sϵ1 ( ζ2dL ζ1+ζ2 - dL 2 )+ ζ1Θ1sϵ0dL ζ1+ζ2 - ζ2Θ1sϵ0dL ζ1+ζ2 - 1 2 Θ0sϵ1dL0=0, (18) Θ1J1s2- ζ1 2 k0dLW(0,s) ζ1+ζ2 + ζ2 2 k0dLW(d1,s) ζ1+ζ2 - ζ1sϵ0dLW(0,s) ζ1+ζ2 + ζ2sϵ0dLW(d1,s) ζ1+ζ2 + ζ1 2 k0dLW1 0(s) ζ1+ζ2 - ζ2 2 k0dLW1 0(s) ζ1+ζ2 + ζ1sϵ0dLW1 0(s) ζ1+ζ2 - ζ2sϵ0dLW1 0(s) ζ1+ζ2 + ζ1 3 Θ1k0dL 2 (ζ1+ζ2) 2 + ζ2 3 Θ1k0dL 2 (ζ1+ζ2) 2 +Θ1M1s2 ( ζ2dL ζ1+ζ2 - dL 2 ) 2 + ζ1 2 Θ1sϵ0dL 2 (ζ1+ζ2) 2 + ζ2 2 Θ1sϵ0dL 2 (ζ1+ζ2) 2 =0, (19) M0s2W0 0(s)+2k1W0 0(s)-k1W1 0(s)-k1W2 0(s)+2sϵ1W0 0(s)-sϵ1W1 0(s) -sϵ1W2 0(s)+Θ1k1 ( ζ2dL ζ1+ζ2 - dL 2 ) -Θ2k1 ( ζ3dL ζ3+ζ4 - dL 2 )+Θ1sϵ1 ( ζ2dL ζ1+ζ2 - dL 2 ) -Θ2sϵ1 ( ζ3dL ζ3+ζ4 - dL 2 )=0, (20) Θ0J0s2+ 1 4 dL0(k1+sϵ1) [ 2(ζ1+ζ2)(ζ3+ζ4)(W2 0(s)-W1 0(s))+dL((ζ1+ζ2)(ζ3-ζ4)Θ2-(ζ1-ζ2)(ζ3+ζ4)Θ1) (ζ1+ζ2)(ζ3+ζ4) +2Θ0dL0]=0. (21) The beam-foundation reaction is mathematically represented by the function F(k,s), which appears on the left-hand side of Eq. (15). We obtain the following system of algebraic equations with respect to the Laplace-displacements of the contact points ws(0,s), ws(d1,s), ws(d2,s), ws(d3,s) after using the inverse Fourier transform over k to the system of Eqs. (34-40), and eliminating W0 0, W 1 0 , W2 0, Θ0, Θ 1 , Θ2 and subsequently setting ξ=0, ξ=d1, ξ=d2, ξ=d3 𝐆 ∙ws[(0, d1, d2, d3), s]=0 ⇔ [ g 11 ⋯ g 14 ⋮ ⋱ ⋮ g 41 ⋯ g 44 ] { ws(0,s) ws(d1,s) ws(d2,s) ws(d3,s)} ={ 0 0 0 0 }. (22) We get the system of integrals that can be seen in the Eq. (22) as a result of the application of the inverse Fourier transform. System of algebraic Eqs. (15-21) possesses a non-trivial solution provided that the determinant of the matrix G satisfies the condition det G=0, (23) In the following sections the model stability is analyzed by investigating the eigenvalues of the characteristic Eq. (23). 780 V. STOJANOVIĆ, L. DIMITROV, P. TOMOV, D. LI, V. NIKOLIĆ 3. SOLUTION METHOD AND STABILITY CRITERIA Having at least one root of the characteristic Eq. (23) with a positive real part can lead to instability in the vibration of the complex moving oscillator. This equation is an integral equation corresponding to the Laplace variables and the issue of finding its complex roots is crucial for the model under study. To simplify the solution, the D-composition method, presented in Stojanović et al. [18], was used. In this approach, k0 is treated as a complex parameter, and the straight lines = iΩ is mapped onto the complex plane k0. The D- decomposition divides the complex plane k0 into regions with different numbers of roots for Eq. (23), focusing on those with Re(s)>0. By using the D-decomposition method we can determine a relative variation of the number of the unstable roots possessing the stiffness k0, however not the number itself. In order to find the number of unstable roots for any single value of stiffness k0 we are applying the argument principle, explained in Stojanović et al. [18]. Parametric analyses are performed for the super-critical case when the complex moving oscillator system’s velocity is higher than the certain critical velocity vcr, because in sub-critical case (v< vcr), there seems not to be any crossing point of the D- decomposition curves and the positive part of the real axes. This may indicate that in the sub-critical case the number of unstable roots (roots with a positive real part) of the charactieristic equation does not change with the stiffness of the supports of the complex moving oscillator k0. On the other hand, in the super-critical case, the number of unstable roots starts varying the moment when k0 passes the critical values of the stiffness. 4. RESULTS AND DISCUSSION The D-decomposition curves were derived using the same model parameters as in Stojanović et al. [18]. The referent case (referred to as Case 1) represents the traditional primary suspension system, where the stiffness k0 is defined with factors ζ1=ζ2=ζ3=ζ4=1. Variations in primary suspension stiffness for other cases are outlined in Table 1. For different values of k0, the D-decomposition curves and the corresponding number of unstable roots were computed at supercritical velocity. The dynamic stability of the system was thoroughly analyzed, focusing on the instability that arises due to supercritical velocity. Fig. 2 illustrates the D-decomposition curves for the referent case, while Table 2 provides details on the number of unstable roots for different intervals of k0. Table 1 Varying non-dimensional parameter of the primary suspension Varying factor ζ 1 ζ 2 ζ 3 ζ 4 Case 1 1 1 1 1 Case 2 8 1 1 1 Case 3 1 8 1 8 Case 4 1 8 1 4 Dynamic Stability Analysis of a Coupled Moving Bogie System 781 In this case a stable region (with no unstable roots) was identified within the stiffness range k0∈(3.067, 5.074)×107 N/m. It is important to note that for the selected supercritical velocity, other intervals of k0 exhibit varying numbers of unstable roots, from which the instability intervals were calculated. This reference case represents the classical coupled moving oscillator model. Table 2 Review of instability regions (v=1.12vcr) k0[107N/m] 0 2.573 3.067 5.074 5.076 +∞ Case 1 4 2 0 2 4 Fig. 2 D-decomposition curves (Case 1) In Case 2, where the first suspension has a higher stiffness (ζ1=8, while the other factors remaining ζ2=ζ3=ζ4=1), the stable region shifts to k0∈(2.8834, 5.0793)×107 N/m (see Table 3). The factor ζ1 was selected after numerous numerical experiments showed no significant changes in dynamic stability within the range 1<ζ1<8. Complete stable region is determined according to obtained D-decomposition curves depicted in Fig.3. The complete stable region was determined based on the D-decomposition curves shown in Fig. 3. These results demonstrate that adjusting the stiffness of the suspension can expand the stability range and shift it toward lower stiffness values. This comparison highlights the advantages of the proposed varying stiffness model in controlling stability regions. Table 3 Review of instability regions (v=1.12vcr) k0[107N/m] 0 0.4691 1.5851 2.8365 2.8834 5.0793 7.1094 +∞ Case 2 4 2 4 2 0 2 4 782 V. STOJANOVIĆ, L. DIMITROV, P. TOMOV, D. LI, V. NIKOLIĆ Fig. 3 D-decomposition curves (Case 2) In Case 3, where the stiffness factors were chosen symmetrically (ζ1=ζ3=8 and ζ2=ζ4=1), the stable regions are distributed across two distinct intervals: k0∈(0.5077, 1.6516)×107N/m and k0∈(2.6090, 6.1575)×107 N/m. Compared to Case 1, where the stable region was limited k0∈(3.067, 5.074)×107 N/m, Case 3 shows a significantly larger stability range, as depicted in Fig. 4. Finally, in Case 4 the last suspension has reduced stiffness, with the ζ - factors defined as ζ1=1, ζ2=8, ζ3=1, ζ4=4. The results for this case are presented in Table 5 and illustrated in Fig. 5. In this configuration, the stable regions are found in the intervals k0∈(0.9368, 1.9278)×107 N m⁄ and k0∈(2.3176, 5.7494)×107 N m⁄ . Compared to the reference Case 1, there is a clear expansion of the stability range. However, in comparison to Case 3, the overall stable region is slightly smaller. Nevertheless, the stability zones are shifted to lower values of the stiffness coefficient k0, indicating that stability can be achieved at reduced stiffness. These findings emphasize the potential benefits of varying the stiffness coefficients in the primary suspension system. The results suggest that further exploration and deeper investigation into this approach could yield additional insights and improvements in dynamic stability. Table 4 Review of instability regions (v=1.12vcr) k0[107N/m] 0 0.4394 0.5077 1.6516 2.6090 6.1575 6.2844 +∞ Case 3 4 2 0 2 0 2 4 Dynamic Stability Analysis of a Coupled Moving Bogie System 783 Fig. 4 D-decomposition curves (Case 3) Table 5 Review of instability regions (v=1.12vcr) k0[107N/m] 0 0.4709 0.9368 1.9278 2.3176 5.7494 6.2270 +∞ Case 3 4 2 0 2 0 2 4 Fig. 5 D-decomposition curves (Case 4) 784 V. STOJANOVIĆ, L. DIMITROV, P. TOMOV, D. LI, V. NIKOLIĆ 5. CONCLUSIONS The study presented a detailed analysis of the stability of a complex coupled mechanical system, specifically focusing on the effects of varying stiffness in the primary suspension system under moving loads. Through the use of D-decomposition curves, the stability regions were identified for different configurations of suspension stiffness, revealing significant insights into the dynamics of the system. The results showed that varying the stiffness of the primary suspension significantly affects the stability of the system. In the referent case (case 1), a stable region was identified within the stiffness range k0∈(3.067, 5.074)×107 N/m. However, when the stiffness of the first suspension was increased in case 2, the stability region shifted slightly and expanded, demonstrating the potential benefits of tuning the system’s stiffness parameters. Case 3, where the suspension was symmetrically adjusted, exhibited a significantly larger stability range compared to Case 1, with stability occurring in two distinct regions. This indicates that symmetric variations in stiffness can lead to broader stability domains. And in Case 4, where the stiffness of the last suspension was reduced, the stable regions were shifted to lower stiffness values, allowing for stability at reduced stiffness values. Though the overall stable area was slightly smaller than in Case 3, the findings suggest that varying suspension stiffness can optimize stability even at lower stiffness values. Overall, this research highlights the importance of carefully adjusting suspension stiffness to enhance the dynamic stability of mechanical systems under moving loads. The results suggest that further exploration of suspension stiffness variations could offer valuable strategies for improving the stability and performance of complex mechanical systems in real- world applications, particularly in the context of transportation infrastructure. Acknowledgement: We acknowledge the financial support from Natural Sciences and Engineering Research Council of Canada through discovery grants and from the Government of Canada’s New Frontiers in Research Fund (NFRF) [ NFRFR-2021-00262 ] and research was also financially supported by the Ministry of Science, Technological Development and Innovation of the Republic of Serbia (Contract No. 451-03-47/2023-01/ 200109). REFERENCES 1. Fryba, L., 1999, Vibration of Solids and Structures Under Moving Loads, Telford, London. 2. Stojanović, V., Kozić, P., Petković, M.D., 2017, Dynamic instability and critical velocity of a mass moving uniformly along a stabilized infinity beam, International Journal of Solids and Structures, 108, pp. 164–174. 3. Michaltsos, G.T., 2002, Dynamic behavior of a single-span beam subjected to loads moving with variable speeds, Journal of Sound and Vibration, 258(2), pp.359–372. 4. Kim, S.-M., 2005, Stability and dynamic response of Rayleigh beam–columns on an elastic foundation under moving loads of constant amplitude and harmonic variation, Engineering Structures, 27(6), pp. 869–880. 5. Śniady, P., 2008, Dynamic Response of a Timoshenko Beam to a Moving Force, ASME Journal of Applied Mechanics, 75(2), 024503. 6. Basu, D., Kameswara Rao, N.S.V., 2012, Analytical solutions for Euler–Bernoulli beam on visco-elastic foundation subjected to moving load, International Journal for Numerical and Analytical Methods in Geomechanics, 37(8), pp. 945–960. 7. Froio, D., Rizzi, E., Simões, F.M.F., Pinto Da Costa, A., 2018, Universal analytical solution of the steady- state response of an infinite beam on a Pasternak elastic foundation under moving load, International Journal of Solids and Structures, 132-133, pp. 245–263. 8. Mazilu, T., 2009, Interaction between a moving two-mass oscillator and an infinite homogeneous structure: Green’s functions method, Archive of Applied Mechanics, 80, pp. 909–927. 9. Karimi, A.H., Ziaei-Rad, S., 2015, Vibration analysis of a beam with moving support subjected to a moving mass travelling with constant and variable speed, 29,pp. 372–390. Dynamic Stability Analysis of a Coupled Moving Bogie System 785 10. Verichev, S.N., Metrikine, A.V., 2000, Dynamic rigidity of a beam in a moving contact, Journal of Applied Mechanics and Technical Physics, 41, pp. 1111–1117. 11. Metrikine, A.V., Verichev, S.N., 2001, Instability of vibrations of a moving two-mass oscillator on a flexibly supported Timoshenko beam, Archive of Applied Mechanics, 71(9), pp. 613–624. 12. Verichev, S.N., Metrikine, A.V., 2002, Instability of a railway vehicle moving on a flexibly supported Timoshenko beam, Journal of Sound and Vibration, 253(3), pp. 653–668. 13. Dimitrovová, Z., Varandas, J.N., 2009, Critical velocity of a load moving on a beam with a sudden change of foundation stiffness: Applications to high-speed trains, Computers & Structures, 87(19-20), pp. 1224–1232. 14. Dimitrovová, Z., 2017, New semi-analytical solution for a uniformly moving mass on a beam on a two- parameter visco-elastic foundation, International Journal of Mechanical Sciences, 127, pp. 142–162. 15. Erbas, B., Kaplunov, J., Prikazchikov, D.A., Sahin, O., 2016, The near-resonant regimes of a moving load in a three-dimensional problem for a coated elastic half-space, Mathematics and Mechanics of Solids, 22, pp. 89–100. 16. Svedholm, C., Zangeneh, A., Pacoste, C., François, S., Karoumi, R., 2016, Vibration of damped uniform beams with general end conditions under moving loads, Engineering Structures, 126,pp. 40–52. 17. Lee, K.Y., Renshaw, A.A., 2000, Solution of the Moving Mass Problem Using Complex Eigenfunction Expansions, ASME Journal of Applied Mechanics, 67(4), 823. 18. Stojanović, V., Deng, J., Milić, D., Petković, M.D., 2023, Dynamics of moving coupled objects with stabilizers and unconventional couplings, Journal of Sound and Vibration, 570, 118020. 19. Stojanović, V., Petković, M., 2013, Moment Lyapunov Exponents and Stochastic Stability of a Three- Dimensional System on Elastic Foundation Using a Perturbation Approach, ASME Journal of Applied Mechanics, 80(5), 051009. 20. Stojanović, V., Petković, M.D., 2016, Nonlinear dynamic analysis of damaged Reddy–Bickford beams supported on an elastic Pasternak foundation, Journal of Sound and Vibration, 385, pp. 239–266.