# Numerical Investigation into the Critical Speed and Frequency of the Hunting Motion in Railway Vehicle System

Authors: Jianfeng Sun; Maoru Chi; Wubin Cai; Xuesong Jin
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The critical speed and hunting frequency are two basic research objects of vehicle system dynamics and have a significant influence on the dynamic performance. A lateral dynamic model with 17 degrees of freedom was established in this study to investigate the critical speed and hunting frequency of a high-speed railway vehicle. The nonlinearities of wheel/rail contact geometry, creep forces, and yaw damper were all considered. A heuristic nonlinear creep model was employed to estimate the contact force between the wheel and the rail. The Maxwell model, which covers the influence of the stiffness characteristic, is used to simulate the yaw damper. To reflect the blow-off of the yaw damper, the damping coefficient is described by stages. Based on the mathematical model, the combined effects of vehicle parameters on the critical speed in the straight line and hunting frequency of the wheelset were investigated innovatively. The novel phenomenon that the hunting frequency exhibits a sudden increase from a smaller value to a larger value when the blow-off of the yaw damper occurs was discovered during the calculations. The extents to which various parameters affect the critical speed and hunting frequency are clear on the basis of the numerical results. Moreover, all of the parameter values were divided into three sections to determine the sensitive range for the critical speed and hunting frequency. The results show that the first section of values plays the decisive role on both the critical speed and the hunting frequency for all parameters analyzed. The investigation in this paper enriches the study of hunting stability and gives some ideas to probably solve the abnormal vibrations during the actual operation.


---

## Body

## 1. Introduction

Stability analysis, referred to as one of the classical research areas in vehicle system dynamics, is studied throughout the development of railway vehicles and plays an important role in the development of high-speed trains in China [1]. The primary subject of a stability analysis is hunting motion, which is characterized by self-excited lateral-yaw oscillations [2]. The creep forces and wheel conicity are regarded as the main causes of hunting instability [3]. The hunting motion occurs at a certain speed called the critical speed, which determines the actual operating speed of the railway vehicle. Once the vehicle speed exceeds the critical speed, the vibration amplitude increases with increasing speed and eventually causes a violent lateral swing of the wheelset, which deteriorates the ride comfort, makes the vehicle prone to derailment, damages the track, and induces fatigue failure in the vehicle structure. The hunting motion is therefore the foundation of vehicle system dynamics and has attracted increasing attention from researchers and engineers.Since Stephenson found the hunting motion phenomenon, researchers have conducted extensive studies on this subject from different aspects [4]. The objects of the study expanded gradually from a single wheelset or bogie to a vehicle and even to an entire train. With the increase in the degrees of freedom, the model became more complex and the factors considered were more comprehensive. Early works generally used linear models which correspond to systems of linear differential equations [5–7]. The critical speed is defined as the lowest speed at which the system characteristic equation has an unstable root, and the influences of different parameters on the hunting stability thus can be analyzed using roots of the characteristic equation. However, with the development of nonlinear dynamics, it was recognized that linear models offer a limited view of vehicle dynamics and cannot describe the entire characteristics of hunting motion. Therefore, the subsequent studies by Law and Brand [8], Hedrick and Arslan [9], and many other researchers extended the linear models to include the nonlinear factors of the vehicle systems, such as the wheel/rail contact geometry, creep forces, and suspension elements. The bifurcation theory provides the dynamic behavior of nonlinear systems when the system parameters change and plays an important role in the analysis of nonlinear system. Going beyond the classical linear stability theory, True [10–14] calculated the unknown solution after the bifurcation point and plotted the diagram of subcritical Hopf bifurcation. Subsequently, a group of scientists [15–19] further studied bifurcation by using numerical and analytical methods. They considered nonlinear creep-creep forces as well as nonlinear contact geometry at the wheel/rail interface. Several nonlinear dynamic phenomena, including chaotic motion, were observed in their models and the results demonstrated that the nonlinear elements in railway vehicles significantly affected bifurcation features. However, these studies assumed that suspension elements never left their linear range, and their nonlinearities were not involved in their studies. The yaw damper, which is a hydraulic system connecting the bogie frame with the car body longitudinally, is the most important component to prevent unstable oscillations of high-speed vehicles. Considering the nonlinearity of yaw dampers at the primary suspension, Ahmadian and Yang [20, 21] analyzed Hopf bifurcation in a wheelset and a locomotive bogie. Applying the center manifold theorem and the method of normal form, Yan [22] presented the influence of the nonlinear yaw damper on the bifurcation type of the bogie. Polach [23] compared the blow-off forces of a vehicle with yaw dampers and that without a yaw damper and discussed bifurcation characteristics using two bifurcation analysis methods. Uyulan [24] analyzed and compared the Hopf bifurcation behavior of a two-axle railway bogie and a dual wheelset by considering the nonlinearities of creep forces and yaw damper forces. However, their studies ignored or did not analyze the oil stiffness of the yaw damper, and the stiffness behavior of the yaw damper actually had a considerable influence on the dynamic characteristics of the yaw damper [25–27].The bogie hunting motion, regarded as a special model in the vehicle system, has a significant influence on the vehicle stability. The damping ratio of this model largely determines whether the hunting motion converges to the equilibrium or diverges to a stable limit cycle. In addition to the damping ratio, the frequency of the hunting motion should also be given enough attention, as it could affect the abnormal vibration of the vehicle on many situations. Both car body hunting [28] and car body shaking [29], observed in daily operation in China, are relevant to the hunting frequency. Furthermore, the traditional formula of the hunting wavelength derived by Klingle [30] is no longer suitable for the analysis of high-speed vehicles. Apart from the equivalent conicity, many other parameters in the vehicle also have large influence on the hunting frequency. The nonlinearities of the wheel/rail contact geometry, creep forces, and suspension elements are also needed to be taken into consideration during the calculation of hunting frequency. However, the corresponding study on the hunting frequency is rare at present.In this paper, a lateral dynamic model of the high-speed vehicle was established, and the nonlinearities of wheel/rail contact geometry, creep forces, and yaw damper forces were taken into consideration. The parameters involved in the nonlinear wheel/rail contact relationship were attained using polynomial interpolation. The nonlinear relationship between creepage and creep force based on a heuristic nonlinear creep model was employed for estimating the contact forces between the wheel and the rail. To integrate the stiffness behavior, the yaw damper was modeled as a Maxwell element, whose damping coefficient is described by stages. The influences of the vehicle parameters on the critical speed and hunting frequency were studied together. The ranges of the parameters studied were divided into three sections and the influence proportion of each section was subsequently compared.


## 2. Mathematical Model of a High-Speed Railway Vehicle

### 2.1. Description of the System

Figure1 shows a typical high-speed vehicle schematic: (a) denotes the top view and (b) denotes the front view. The vehicle is modeled as a four-axle mass-spring-damping multi-rigid-body system including a car body, two bogie frames, four wheelsets, and two-stage suspensions. The wheelsets have 2 DOFs, including motions of traversey and yawφ. The bogie frames and the car body both have 3 DOFs, including motions of traversey, yawφ, and rollΦ, with respect to their mass centroids. Therefore, the vehicle system has 17 DOFs. The vehicle proceeds at a constant speedv on a straight track.Figure 1
Model of a high-speed vehicle. (a) Top view. (b) Front view.
(a) (b)The subscriptsc,t, andw represent the car body, bogie frame, and wheelset, respectively, and the subscriptsp ands denote the primary and secondary suspensions, respectively. It should be noted that physical quantities, such as Kpx, Kpy, Cpx, Cpy, etc., are defined in the nomenclature displayed in the Appendix.This study used the common way of representing a yaw damper: a linear spring in series with a viscous damper.Cd denotes the damping generated as a result of the resistance of the fluid passing through the various damping valves, andKd denotes the series stiffness representing an elastic characteristic which includes the flexibility related to the connection joints of the damper and oil compressibility.It is convenient to derive the equations of motion for the lateral dynamic vehicle system by applying the d’ Alembert-Lagrange principle.


### 2.2. Governing Equations of Motion

The governing equations of motion for the lateral displacementyc, yaw angleφc, and roll angleΦc of the car body are given, respectively, by the following:(1)Mcy¨c=Fsy1+Fsy2(2)Jczφ¨c=Msz1+Msz2+Fsy1-Fsy2l(3)JcxΦ¨c=Msx1+Msx2+Msy1+Msy2+MgcNote that subscripts 1 and 2 in Equations (1)–(3) indicate that the corresponding physical properties relate to the front and rear bogie frames, respectively. Furthermore, the double dots above the quantities on the left-hand side of Equations (1)–(3) denote differentiation with respect to the time variablet.As shown in Figure1, the suspension forces acting on the car body in the lateral direction,Fsy1 andFsy2, the suspension moments acting on the car body in the vertical direction,Msz1 andMsz2, and the suspension moments acting on the car body in the longitudinal direction,Msx1,Msx2,Msy1,Msy2, andMgc, are given by(4)Fsy1=2Ksyyt1+h3Φt1-yc-lφc+h1Φc+2Csyy˙t1+h5Φ˙t1-y˙c-lφ˙c+h2Φ˙c(5)Fsy2=2Ksyyt2+h3Φt2-yc-lφc+h1Φc+2Csyy˙t2+h5Φ˙t2-y˙c-lφ˙c+h2Φ˙c(6)Msz1=2Ksxb22φt1-φc+2b4Fd1(7)Msz2=2Ksxb22φt2-φc+2b4Fd2whereFd1 andFd2 are the damping forces produced by the yaw dampers of the front and rear bogie frames, and they will be listed in the next section. Consider(8)Msx1=2b22KszΦt1-Φc+2b32CszΦ˙t1-Φ˙c(9)Msx2=2b22KszΦt2-Φc+2b32CszΦ˙t2-Φ˙c(10)Msy1=-2Ksyh1yt1+h3Φt1-yc-lφc+h1Φc-2Csyh2y˙t1+h5Φ˙t1-y˙c-lφ˙c+h2Φ˙c(11)Msy2=-2Ksyh1yt2+h3Φt2-yc-lφc+h1Φc-2Csyh2y˙t2+h5Φ˙t2-y˙c-lφ˙c+h2Φ˙c(12)Mgc=Mcgh1ΦcThe governing equations of motion for the lateral displacementyti, yaw angleφti, and roll angleΦti of the bogie frames are given, respectively, by(13)Mtiy¨ti=-Fsyi+Fpyi1+Fpyi2(14)Jtziφ¨ti=-Mszi+Mpzi1+Mpzi2+Fpyi1-Fpyi2l1(15)JtxiΦ¨ti=-Msxi+Mpxi+Msi+Fpyi1+Fpyi2h4+MgtiNote that subscripts i=1,2 in Equations (13)–(15) indicate the corresponding physical properties related to the front and rear bogie frame, and the subscripts 1 and 2 denote the corresponding properties related to the front and rear wheelset in a bogie.Meanwhile, with respect to the bogie frames, the suspension forces acting in the lateral direction,Fpy1 andFpy2, the suspension moments acting in the vertical direction,Mpz1 andMpz2, the suspension moments acting in the longitudinal direction,Mpxi,Msi, andMgti, are given by(16)Fpyi1=2Kpyyw1-yt1-l1φt1+h4Φt1+2Cpyy˙w1-y˙t1-l1φ˙t1+h4Φ˙t1i=12Kpyyw3-yt2-l1φt2+h4Φt2+2Cpyy˙w3-y˙t2-l1φ˙t2+h4Φ˙t2i=2(17)Fpyi2=2Kpyyw2-yt1+l1φt1+h4Φt1+2Cpyy˙w2-y˙t1+l1φ˙t1+h4Φ˙t1i=12Kpyyw4-yt2+l1φt2+h4Φt2+2Cpyy˙w4-y˙t2+l1φ˙t2+h4Φ˙t2i=2(18)Mpzi1=2Kpxb12φw1-φt1+2Cpyl12φ˙w1-φ˙t1i=12Kpxb12φw3-φt2+2Cpyl12φ˙w3-φ˙t2i=2(19)Mpzi2=2Kpxb12φw2-φt1+2Cpyl12φ˙w2-φ˙t1i=12Kpxb12φw4-φt2+2Cpyl12φ˙w4-φ˙t2i=2(20)Mpxi=-2Kpzb12Φti-2Cpzb12Φ˙ti(21)Msi=2Ksyh3yc+lφc-yt1-h3Φt1-h1Φc+2Csyh5y˙c+lφ˙c-y˙t1-h5Φ˙t1-h2Φ˙ci=12Ksyh3yc-lφc-yt2-h3Φt2-h1Φc+2Csyh5y˙c-lφ˙c-y˙t2-h5Φ˙t2-h2Φ˙ci=2(22)Mgti=Mtgh4ΦtiThe governing equations of motion for the lateral displacementywi and yaw angle φwi of the wheelsets are given, respectively, by(23)Mwy¨wij=-Fpyij+Fryij-Fgij(24)Jwzφ¨wij=-2Kpxb12φwij-φti+Mrzij+Mgijwhere subscripts i=1,2 in Equations (23) and (24) are defined as above, and subscriptsj=1, 2 denote the corresponding physical properties related to the front and rear wheelset. The nonlinear terms, the lateral creep forceFryij, the gravity restoring forceFgij, the moment in the vertical direction produced by the longitudinal creep force and spin momentMrzij, and the moment in the vertical direction produced by the gravity restoring forceMgijare listed below.


### 2.3. Nonlinear Terms

#### 2.3.1. Nonlinearity of Creep Force

A single wheelset model consisting of 2 DOFs is established, and the front view of the wheelset is depicted in Figure2. In the figure,a denotes half of the rolling cycle gauge;v is the forward speed of the wheelset; andrL andrR denote the radius of the left and right wheels at the contact points on the rail, respectively. The creepages of the left and right wheels in the longitudinal, lateral, and spin directions of the contact plane can be derived as follows:(25)ζLt=vcos⁡φ+y˙sin⁡φ-aφ˙-rLΩvζLa=y˙cos⁡φ-vsin⁡φcos⁡δLvζLs=φ˙cos⁡δL-Ωsin⁡δLvζRt=vcos⁡φ+y˙sin⁡φ+aφ˙-rRΩvζRa=y˙cos⁡φ-vsin⁡φcos⁡δRvζRs=φ˙cos⁡δR+Ωsin⁡δRvwhereΩis the rolling angular velocity of the wheelset relative to its own axle; the subscriptsL andR represent the left and right wheels, respectively, while the subscriptst,a, ands represent the longitudinal and lateral directions, respectively.Figure 2
Front view of the wheelset.By using Kalker’s linear creep theory [31], the longitudinal creep force, lateral creep force, and spin creep moment can be written in the following forms:(26)Fkt=-f11ζktFka=-f22ζka-f23ζksMkz=f23ζka-f33ζkswhere f11,f22,f23, and f33 are the creep coefficients and can be calculated referring to Shabana et al. [32] Note that subscripts k=L,R in Equation (26) indicate the corresponding physical properties related to the left and right wheels. A heuristic nonlinear creep model [33] that combines Kalker’s linear creep theory with a creep force saturation representation is used for computing creep force. The creep forces at each wheel are given by(27)Fkt1=αFktFka1=αFkawhere the saturation constant α is given by(28)α=1ββ-β23+β327forβ≤31βforβ>3where β=βL+βR/2,βk=Fkt2+Fka2/μN,μ is the wheel/rail friction coefficient, andN is the normal force to the contact plane. The contact plane creep forces and moment are then transformed to the wheelset coordinates by(29)FLy=FLa1cos⁡δL,MLz1=MLzcos⁡δLFRy=FRa1cos⁡δR,MRz1=MRzcos⁡δRThe forces and moments acting on the track coordinates can be obtained as follows:(30)Fry=FLy+FRycos⁡φ+FLt1+FRt1sin⁡φFg=Wg2tan⁡δR-tan⁡δLMrz=FRt1-FLt1a+MLz1+MRz1Mg=Wgasin⁡φ2tan⁡δL+tan⁡δRcorresponding to each wheelset as(31)Fryij=Fryywij,y˙wij,φwij,φ˙wijFgij=Fgywij,y˙wij,φwij,φ˙wijMrzij=Mrzywij,y˙wij,φwij,φ˙wijMgij=Mgywij,y˙wij,φwij,φ˙wij


#### 2.3.2. Nonlinearity of Wheel/Rail Contact

With respect to the rolling radiirL andrR and the contact anglesδL andδR, they can be taken from tables as functions of excursiony [34]. However, the cost of computation is high when using this method. Function fitting is a good alternative and has been adopted by some authors [18, 22]. The relationship between the equivalent radius of the wheel and the lateral displacement is given by(32)rL=r0+λ1y+λ2y2+λ3y3+···rR=r0+λ1-y+λ2-y2+λ3-y3+···wherer0 is the distance between the contact point and the center of the wheel without displacement,λ1 is the linear gradient on the conical tread, and λi(i⩾2) is the nonlinear correction factors of the geometry relationship between the wheel and the rail.The relationship between the contact angles and the lateral displacement is given by(33)δL=d0+d1y+d4y4+d5y5δR=d0+d1-y+d4-y4+d5-y5One more prerequisite mentioned in [35] has also been applied in this study as(34)tan⁡δR-tan⁡δL=∑s=0ne2sy2s+1The sine and cosine functions are expanded in the Taylor series, and only the first two terms are considered and all higher order derivatives are ignored. Similarly, the tangent function in the Taylor series is expanded to the first two terms in the calculation ofMg.The wheel profile type of S1002CN and rail profile type of UIC60 are selected in this study, and the parameters ofλi(i=1,2,3…) and e2s(s=0,1…n) are selected from the results in Yan and Zeng [22]. Based on the relationship between the lateral displacementy and the contact anglesδL orδR, the values ofd0=0.0544,d1=-9.45,d4=3.8e7, andd5=-3.58e9 can be obtained using fifth-order polynomial interpolation. A comparison between the real data and the result of polynomial interpolation is given in Figure 3.Figure 3
Relationship between lateral displacementy and contact angle.


#### 2.3.3. Nonlinearity of Damping Force

The yaw damper considered here is modeled as a Maxwell element, which is used as a part model introduced in Pracny et al. [25] and can describe the effects of hysteresis in a given frequency range. As mentioned previously, it is based on a linear spring and a viscous damper connected in series, as shown in Figure 4(a), whereFd denotes the damping force,x denotes the total displacement, andxd denotes the displacement of the damper. To consider the curve negotiating property of the vehicle and protect the damper, a yaw damper possessing the blow-off characteristic is applied, and the damping can be described by stages, as shown in Figure 4(b), where Cd1 is the first stage damping before blow-off, Cd2 is the second stage damping after blow-off,vu is the blow-off velocity, andFu is the blow-off force.Figure 4
Model and nonlinearity of yaw damper. (a) Maxwell model. (b) Damping characteristic curve.
(a) (b)The relationship for the spring force can be written as follows:(35)Fd=Kdx-xdAs the damper is assumed to be nonlinear, the rate of this displacement is given by(36)x˙d=FdCd1x˙d≤vuFd-FuCd2+vux˙d>vuDifferentiating Equation (35) with respect to time and replacing the damper velocity with Equation (36), we obtain the following differential equation:(37)F˙d=Kdx˙-KdCd1Fdx˙d≤vuKdx˙-KdFd-FuCd2+vux˙d>vuwhich can be integrated in closed form. The analytical solution corresponding to the x˙t=0 for t⩽0 is given by (38)Fdt=Kd∫0te-t-sKd/Cd1x˙sdsx˙d≤vuKd∫0te-t-sKd/Cd1x˙s+FuCd2-vudsx˙d>vuEquation (38) shows that the current force F(t) is a function of the entire time history of the excitation between s=0 and the actual time s=t. Thus, replacing the integral with the sum and discretization of the time tk=k∆t and si=i∆t, we obtain(39)Fdtk≈KdΔt∑i=0ke-k-iΔtKd/Cd1x˙iΔtx˙d≤vuKdΔt∑i=0ke-k-iΔtKd/Cd1x˙iΔt+FuCd2-vux˙d>vuThus far, all the terms on the right-hand side of Equations (1)–(24) are formulated, and the time responses of the vehicle components can be obtained by numerical calculation.


## 2.1. Description of the System

Figure1 shows a typical high-speed vehicle schematic: (a) denotes the top view and (b) denotes the front view. The vehicle is modeled as a four-axle mass-spring-damping multi-rigid-body system including a car body, two bogie frames, four wheelsets, and two-stage suspensions. The wheelsets have 2 DOFs, including motions of traversey and yawφ. The bogie frames and the car body both have 3 DOFs, including motions of traversey, yawφ, and rollΦ, with respect to their mass centroids. Therefore, the vehicle system has 17 DOFs. The vehicle proceeds at a constant speedv on a straight track.Figure 1
Model of a high-speed vehicle. (a) Top view. (b) Front view.
(a) (b)The subscriptsc,t, andw represent the car body, bogie frame, and wheelset, respectively, and the subscriptsp ands denote the primary and secondary suspensions, respectively. It should be noted that physical quantities, such as Kpx, Kpy, Cpx, Cpy, etc., are defined in the nomenclature displayed in the Appendix.This study used the common way of representing a yaw damper: a linear spring in series with a viscous damper.Cd denotes the damping generated as a result of the resistance of the fluid passing through the various damping valves, andKd denotes the series stiffness representing an elastic characteristic which includes the flexibility related to the connection joints of the damper and oil compressibility.It is convenient to derive the equations of motion for the lateral dynamic vehicle system by applying the d’ Alembert-Lagrange principle.


## 2.2. Governing Equations of Motion

The governing equations of motion for the lateral displacementyc, yaw angleφc, and roll angleΦc of the car body are given, respectively, by the following:(1)Mcy¨c=Fsy1+Fsy2(2)Jczφ¨c=Msz1+Msz2+Fsy1-Fsy2l(3)JcxΦ¨c=Msx1+Msx2+Msy1+Msy2+MgcNote that subscripts 1 and 2 in Equations (1)–(3) indicate that the corresponding physical properties relate to the front and rear bogie frames, respectively. Furthermore, the double dots above the quantities on the left-hand side of Equations (1)–(3) denote differentiation with respect to the time variablet.As shown in Figure1, the suspension forces acting on the car body in the lateral direction,Fsy1 andFsy2, the suspension moments acting on the car body in the vertical direction,Msz1 andMsz2, and the suspension moments acting on the car body in the longitudinal direction,Msx1,Msx2,Msy1,Msy2, andMgc, are given by(4)Fsy1=2Ksyyt1+h3Φt1-yc-lφc+h1Φc+2Csyy˙t1+h5Φ˙t1-y˙c-lφ˙c+h2Φ˙c(5)Fsy2=2Ksyyt2+h3Φt2-yc-lφc+h1Φc+2Csyy˙t2+h5Φ˙t2-y˙c-lφ˙c+h2Φ˙c(6)Msz1=2Ksxb22φt1-φc+2b4Fd1(7)Msz2=2Ksxb22φt2-φc+2b4Fd2whereFd1 andFd2 are the damping forces produced by the yaw dampers of the front and rear bogie frames, and they will be listed in the next section. Consider(8)Msx1=2b22KszΦt1-Φc+2b32CszΦ˙t1-Φ˙c(9)Msx2=2b22KszΦt2-Φc+2b32CszΦ˙t2-Φ˙c(10)Msy1=-2Ksyh1yt1+h3Φt1-yc-lφc+h1Φc-2Csyh2y˙t1+h5Φ˙t1-y˙c-lφ˙c+h2Φ˙c(11)Msy2=-2Ksyh1yt2+h3Φt2-yc-lφc+h1Φc-2Csyh2y˙t2+h5Φ˙t2-y˙c-lφ˙c+h2Φ˙c(12)Mgc=Mcgh1ΦcThe governing equations of motion for the lateral displacementyti, yaw angleφti, and roll angleΦti of the bogie frames are given, respectively, by(13)Mtiy¨ti=-Fsyi+Fpyi1+Fpyi2(14)Jtziφ¨ti=-Mszi+Mpzi1+Mpzi2+Fpyi1-Fpyi2l1(15)JtxiΦ¨ti=-Msxi+Mpxi+Msi+Fpyi1+Fpyi2h4+MgtiNote that subscripts i=1,2 in Equations (13)–(15) indicate the corresponding physical properties related to the front and rear bogie frame, and the subscripts 1 and 2 denote the corresponding properties related to the front and rear wheelset in a bogie.Meanwhile, with respect to the bogie frames, the suspension forces acting in the lateral direction,Fpy1 andFpy2, the suspension moments acting in the vertical direction,Mpz1 andMpz2, the suspension moments acting in the longitudinal direction,Mpxi,Msi, andMgti, are given by(16)Fpyi1=2Kpyyw1-yt1-l1φt1+h4Φt1+2Cpyy˙w1-y˙t1-l1φ˙t1+h4Φ˙t1i=12Kpyyw3-yt2-l1φt2+h4Φt2+2Cpyy˙w3-y˙t2-l1φ˙t2+h4Φ˙t2i=2(17)Fpyi2=2Kpyyw2-yt1+l1φt1+h4Φt1+2Cpyy˙w2-y˙t1+l1φ˙t1+h4Φ˙t1i=12Kpyyw4-yt2+l1φt2+h4Φt2+2Cpyy˙w4-y˙t2+l1φ˙t2+h4Φ˙t2i=2(18)Mpzi1=2Kpxb12φw1-φt1+2Cpyl12φ˙w1-φ˙t1i=12Kpxb12φw3-φt2+2Cpyl12φ˙w3-φ˙t2i=2(19)Mpzi2=2Kpxb12φw2-φt1+2Cpyl12φ˙w2-φ˙t1i=12Kpxb12φw4-φt2+2Cpyl12φ˙w4-φ˙t2i=2(20)Mpxi=-2Kpzb12Φti-2Cpzb12Φ˙ti(21)Msi=2Ksyh3yc+lφc-yt1-h3Φt1-h1Φc+2Csyh5y˙c+lφ˙c-y˙t1-h5Φ˙t1-h2Φ˙ci=12Ksyh3yc-lφc-yt2-h3Φt2-h1Φc+2Csyh5y˙c-lφ˙c-y˙t2-h5Φ˙t2-h2Φ˙ci=2(22)Mgti=Mtgh4ΦtiThe governing equations of motion for the lateral displacementywi and yaw angle φwi of the wheelsets are given, respectively, by(23)Mwy¨wij=-Fpyij+Fryij-Fgij(24)Jwzφ¨wij=-2Kpxb12φwij-φti+Mrzij+Mgijwhere subscripts i=1,2 in Equations (23) and (24) are defined as above, and subscriptsj=1, 2 denote the corresponding physical properties related to the front and rear wheelset. The nonlinear terms, the lateral creep forceFryij, the gravity restoring forceFgij, the moment in the vertical direction produced by the longitudinal creep force and spin momentMrzij, and the moment in the vertical direction produced by the gravity restoring forceMgijare listed below.


## 2.3. Nonlinear Terms

### 2.3.1. Nonlinearity of Creep Force

A single wheelset model consisting of 2 DOFs is established, and the front view of the wheelset is depicted in Figure2. In the figure,a denotes half of the rolling cycle gauge;v is the forward speed of the wheelset; andrL andrR denote the radius of the left and right wheels at the contact points on the rail, respectively. The creepages of the left and right wheels in the longitudinal, lateral, and spin directions of the contact plane can be derived as follows:(25)ζLt=vcos⁡φ+y˙sin⁡φ-aφ˙-rLΩvζLa=y˙cos⁡φ-vsin⁡φcos⁡δLvζLs=φ˙cos⁡δL-Ωsin⁡δLvζRt=vcos⁡φ+y˙sin⁡φ+aφ˙-rRΩvζRa=y˙cos⁡φ-vsin⁡φcos⁡δRvζRs=φ˙cos⁡δR+Ωsin⁡δRvwhereΩis the rolling angular velocity of the wheelset relative to its own axle; the subscriptsL andR represent the left and right wheels, respectively, while the subscriptst,a, ands represent the longitudinal and lateral directions, respectively.Figure 2
Front view of the wheelset.By using Kalker’s linear creep theory [31], the longitudinal creep force, lateral creep force, and spin creep moment can be written in the following forms:(26)Fkt=-f11ζktFka=-f22ζka-f23ζksMkz=f23ζka-f33ζkswhere f11,f22,f23, and f33 are the creep coefficients and can be calculated referring to Shabana et al. [32] Note that subscripts k=L,R in Equation (26) indicate the corresponding physical properties related to the left and right wheels. A heuristic nonlinear creep model [33] that combines Kalker’s linear creep theory with a creep force saturation representation is used for computing creep force. The creep forces at each wheel are given by(27)Fkt1=αFktFka1=αFkawhere the saturation constant α is given by(28)α=1ββ-β23+β327forβ≤31βforβ>3where β=βL+βR/2,βk=Fkt2+Fka2/μN,μ is the wheel/rail friction coefficient, andN is the normal force to the contact plane. The contact plane creep forces and moment are then transformed to the wheelset coordinates by(29)FLy=FLa1cos⁡δL,MLz1=MLzcos⁡δLFRy=FRa1cos⁡δR,MRz1=MRzcos⁡δRThe forces and moments acting on the track coordinates can be obtained as follows:(30)Fry=FLy+FRycos⁡φ+FLt1+FRt1sin⁡φFg=Wg2tan⁡δR-tan⁡δLMrz=FRt1-FLt1a+MLz1+MRz1Mg=Wgasin⁡φ2tan⁡δL+tan⁡δRcorresponding to each wheelset as(31)Fryij=Fryywij,y˙wij,φwij,φ˙wijFgij=Fgywij,y˙wij,φwij,φ˙wijMrzij=Mrzywij,y˙wij,φwij,φ˙wijMgij=Mgywij,y˙wij,φwij,φ˙wij


### 2.3.2. Nonlinearity of Wheel/Rail Contact

With respect to the rolling radiirL andrR and the contact anglesδL andδR, they can be taken from tables as functions of excursiony [34]. However, the cost of computation is high when using this method. Function fitting is a good alternative and has been adopted by some authors [18, 22]. The relationship between the equivalent radius of the wheel and the lateral displacement is given by(32)rL=r0+λ1y+λ2y2+λ3y3+···rR=r0+λ1-y+λ2-y2+λ3-y3+···wherer0 is the distance between the contact point and the center of the wheel without displacement,λ1 is the linear gradient on the conical tread, and λi(i⩾2) is the nonlinear correction factors of the geometry relationship between the wheel and the rail.The relationship between the contact angles and the lateral displacement is given by(33)δL=d0+d1y+d4y4+d5y5δR=d0+d1-y+d4-y4+d5-y5One more prerequisite mentioned in [35] has also been applied in this study as(34)tan⁡δR-tan⁡δL=∑s=0ne2sy2s+1The sine and cosine functions are expanded in the Taylor series, and only the first two terms are considered and all higher order derivatives are ignored. Similarly, the tangent function in the Taylor series is expanded to the first two terms in the calculation ofMg.The wheel profile type of S1002CN and rail profile type of UIC60 are selected in this study, and the parameters ofλi(i=1,2,3…) and e2s(s=0,1…n) are selected from the results in Yan and Zeng [22]. Based on the relationship between the lateral displacementy and the contact anglesδL orδR, the values ofd0=0.0544,d1=-9.45,d4=3.8e7, andd5=-3.58e9 can be obtained using fifth-order polynomial interpolation. A comparison between the real data and the result of polynomial interpolation is given in Figure 3.Figure 3
Relationship between lateral displacementy and contact angle.


### 2.3.3. Nonlinearity of Damping Force

The yaw damper considered here is modeled as a Maxwell element, which is used as a part model introduced in Pracny et al. [25] and can describe the effects of hysteresis in a given frequency range. As mentioned previously, it is based on a linear spring and a viscous damper connected in series, as shown in Figure 4(a), whereFd denotes the damping force,x denotes the total displacement, andxd denotes the displacement of the damper. To consider the curve negotiating property of the vehicle and protect the damper, a yaw damper possessing the blow-off characteristic is applied, and the damping can be described by stages, as shown in Figure 4(b), where Cd1 is the first stage damping before blow-off, Cd2 is the second stage damping after blow-off,vu is the blow-off velocity, andFu is the blow-off force.Figure 4
Model and nonlinearity of yaw damper. (a) Maxwell model. (b) Damping characteristic curve.
(a) (b)The relationship for the spring force can be written as follows:(35)Fd=Kdx-xdAs the damper is assumed to be nonlinear, the rate of this displacement is given by(36)x˙d=FdCd1x˙d≤vuFd-FuCd2+vux˙d>vuDifferentiating Equation (35) with respect to time and replacing the damper velocity with Equation (36), we obtain the following differential equation:(37)F˙d=Kdx˙-KdCd1Fdx˙d≤vuKdx˙-KdFd-FuCd2+vux˙d>vuwhich can be integrated in closed form. The analytical solution corresponding to the x˙t=0 for t⩽0 is given by (38)Fdt=Kd∫0te-t-sKd/Cd1x˙sdsx˙d≤vuKd∫0te-t-sKd/Cd1x˙s+FuCd2-vudsx˙d>vuEquation (38) shows that the current force F(t) is a function of the entire time history of the excitation between s=0 and the actual time s=t. Thus, replacing the integral with the sum and discretization of the time tk=k∆t and si=i∆t, we obtain(39)Fdtk≈KdΔt∑i=0ke-k-iΔtKd/Cd1x˙iΔtx˙d≤vuKdΔt∑i=0ke-k-iΔtKd/Cd1x˙iΔt+FuCd2-vux˙d>vuThus far, all the terms on the right-hand side of Equations (1)–(24) are formulated, and the time responses of the vehicle components can be obtained by numerical calculation.


## 2.3.1. Nonlinearity of Creep Force

A single wheelset model consisting of 2 DOFs is established, and the front view of the wheelset is depicted in Figure2. In the figure,a denotes half of the rolling cycle gauge;v is the forward speed of the wheelset; andrL andrR denote the radius of the left and right wheels at the contact points on the rail, respectively. The creepages of the left and right wheels in the longitudinal, lateral, and spin directions of the contact plane can be derived as follows:(25)ζLt=vcos⁡φ+y˙sin⁡φ-aφ˙-rLΩvζLa=y˙cos⁡φ-vsin⁡φcos⁡δLvζLs=φ˙cos⁡δL-Ωsin⁡δLvζRt=vcos⁡φ+y˙sin⁡φ+aφ˙-rRΩvζRa=y˙cos⁡φ-vsin⁡φcos⁡δRvζRs=φ˙cos⁡δR+Ωsin⁡δRvwhereΩis the rolling angular velocity of the wheelset relative to its own axle; the subscriptsL andR represent the left and right wheels, respectively, while the subscriptst,a, ands represent the longitudinal and lateral directions, respectively.Figure 2
Front view of the wheelset.By using Kalker’s linear creep theory [31], the longitudinal creep force, lateral creep force, and spin creep moment can be written in the following forms:(26)Fkt=-f11ζktFka=-f22ζka-f23ζksMkz=f23ζka-f33ζkswhere f11,f22,f23, and f33 are the creep coefficients and can be calculated referring to Shabana et al. [32] Note that subscripts k=L,R in Equation (26) indicate the corresponding physical properties related to the left and right wheels. A heuristic nonlinear creep model [33] that combines Kalker’s linear creep theory with a creep force saturation representation is used for computing creep force. The creep forces at each wheel are given by(27)Fkt1=αFktFka1=αFkawhere the saturation constant α is given by(28)α=1ββ-β23+β327forβ≤31βforβ>3where β=βL+βR/2,βk=Fkt2+Fka2/μN,μ is the wheel/rail friction coefficient, andN is the normal force to the contact plane. The contact plane creep forces and moment are then transformed to the wheelset coordinates by(29)FLy=FLa1cos⁡δL,MLz1=MLzcos⁡δLFRy=FRa1cos⁡δR,MRz1=MRzcos⁡δRThe forces and moments acting on the track coordinates can be obtained as follows:(30)Fry=FLy+FRycos⁡φ+FLt1+FRt1sin⁡φFg=Wg2tan⁡δR-tan⁡δLMrz=FRt1-FLt1a+MLz1+MRz1Mg=Wgasin⁡φ2tan⁡δL+tan⁡δRcorresponding to each wheelset as(31)Fryij=Fryywij,y˙wij,φwij,φ˙wijFgij=Fgywij,y˙wij,φwij,φ˙wijMrzij=Mrzywij,y˙wij,φwij,φ˙wijMgij=Mgywij,y˙wij,φwij,φ˙wij


## 2.3.2. Nonlinearity of Wheel/Rail Contact

With respect to the rolling radiirL andrR and the contact anglesδL andδR, they can be taken from tables as functions of excursiony [34]. However, the cost of computation is high when using this method. Function fitting is a good alternative and has been adopted by some authors [18, 22]. The relationship between the equivalent radius of the wheel and the lateral displacement is given by(32)rL=r0+λ1y+λ2y2+λ3y3+···rR=r0+λ1-y+λ2-y2+λ3-y3+···wherer0 is the distance between the contact point and the center of the wheel without displacement,λ1 is the linear gradient on the conical tread, and λi(i⩾2) is the nonlinear correction factors of the geometry relationship between the wheel and the rail.The relationship between the contact angles and the lateral displacement is given by(33)δL=d0+d1y+d4y4+d5y5δR=d0+d1-y+d4-y4+d5-y5One more prerequisite mentioned in [35] has also been applied in this study as(34)tan⁡δR-tan⁡δL=∑s=0ne2sy2s+1The sine and cosine functions are expanded in the Taylor series, and only the first two terms are considered and all higher order derivatives are ignored. Similarly, the tangent function in the Taylor series is expanded to the first two terms in the calculation ofMg.The wheel profile type of S1002CN and rail profile type of UIC60 are selected in this study, and the parameters ofλi(i=1,2,3…) and e2s(s=0,1…n) are selected from the results in Yan and Zeng [22]. Based on the relationship between the lateral displacementy and the contact anglesδL orδR, the values ofd0=0.0544,d1=-9.45,d4=3.8e7, andd5=-3.58e9 can be obtained using fifth-order polynomial interpolation. A comparison between the real data and the result of polynomial interpolation is given in Figure 3.Figure 3
Relationship between lateral displacementy and contact angle.


## 2.3.3. Nonlinearity of Damping Force

The yaw damper considered here is modeled as a Maxwell element, which is used as a part model introduced in Pracny et al. [25] and can describe the effects of hysteresis in a given frequency range. As mentioned previously, it is based on a linear spring and a viscous damper connected in series, as shown in Figure 4(a), whereFd denotes the damping force,x denotes the total displacement, andxd denotes the displacement of the damper. To consider the curve negotiating property of the vehicle and protect the damper, a yaw damper possessing the blow-off characteristic is applied, and the damping can be described by stages, as shown in Figure 4(b), where Cd1 is the first stage damping before blow-off, Cd2 is the second stage damping after blow-off,vu is the blow-off velocity, andFu is the blow-off force.Figure 4
Model and nonlinearity of yaw damper. (a) Maxwell model. (b) Damping characteristic curve.
(a) (b)The relationship for the spring force can be written as follows:(35)Fd=Kdx-xdAs the damper is assumed to be nonlinear, the rate of this displacement is given by(36)x˙d=FdCd1x˙d≤vuFd-FuCd2+vux˙d>vuDifferentiating Equation (35) with respect to time and replacing the damper velocity with Equation (36), we obtain the following differential equation:(37)F˙d=Kdx˙-KdCd1Fdx˙d≤vuKdx˙-KdFd-FuCd2+vux˙d>vuwhich can be integrated in closed form. The analytical solution corresponding to the x˙t=0 for t⩽0 is given by (38)Fdt=Kd∫0te-t-sKd/Cd1x˙sdsx˙d≤vuKd∫0te-t-sKd/Cd1x˙s+FuCd2-vudsx˙d>vuEquation (38) shows that the current force F(t) is a function of the entire time history of the excitation between s=0 and the actual time s=t. Thus, replacing the integral with the sum and discretization of the time tk=k∆t and si=i∆t, we obtain(39)Fdtk≈KdΔt∑i=0ke-k-iΔtKd/Cd1x˙iΔtx˙d≤vuKdΔt∑i=0ke-k-iΔtKd/Cd1x˙iΔt+FuCd2-vux˙d>vuThus far, all the terms on the right-hand side of Equations (1)–(24) are formulated, and the time responses of the vehicle components can be obtained by numerical calculation.


## 3. Numerical Analysis of the Critical Speed and Hunting Frequency of the Vehicle System

### 3.1. Hopf Bifurcation Characteristics

A bifurcation diagram is a powerful tool for analyzing railway nonlinear dynamics. The railway system expresses the typical Hopf bifurcation, whose diagram is illustrated in Figure5 [17]. In Figure 5, the abscissa represents the vehicle speed while the ordinate represents the amplitude of the limit cycle. The solid and dashed lines denote the stable and unstable limit cycles, respectively. PointA is the Hopf bifurcation point, and the corresponding speedVlin at pointA is defined as the linear critical speed of railway vehicles. PointB is the turning point from a stable equilibrium point to a periodic solution, and the corresponding speed Vnl at point B is defined as the nonlinear critical speed of railway vehicles.Figure 5
Typical types of Hopf bifurcations of a railway vehicle system.For the subcritical Hopf bifurcation, system vibration will result in a stable equilibrium position for any external excitations when the vehicle speedV is less thanVnl. WhenV is larger thanVlin, the system will exhibit a stable periodic motion with a large amplitude regardless of the amplitude of the disturbance. In the intervalVnl <V <Vlin, however, the vibration of the system depends on the amplitude of the external excitation, and the interval is an uncertain region. For the supercritical Hopf bifurcation (b), the system always exhibits a stable periodic motion whenV achieves the linear critical speedVlin. With the increase inV, the amplitude of the hunting motion increases gradually. The supercritical Hopf bifurcation (a) has both subcritical Hopf bifurcation and the supercritical Hopf bifurcation characteristics. The uncertain region ofV extends toVnl ~Vc, in which the solution will converge to equilibrium or a limit cycle with a small amplitude at a small track excitation and diverge to a limit cycle with a large amplitude at a large track excitation. The bifurcation diagram can be calculated using a path following method or a set of numerical simulations [36–38]. The path following method is used in this paper.


### 3.2. Critical Hunting Speed

To investigate the stability of a railway vehicle system in practical engineering, the nonlinear critical hunting speedVnl is generally considered. From a mechanical viewpoint, a system can be regarded as stable if the oscillation of the system decays after a discontinuation excitation. Increasing the vehicle speed slowly, when the system displays a periodic motion at a constant amplitude, the corresponding speed can be regarded as the critical speedVcr, namely,Vnl. The lateral displacement of the front wheelset is used as the output value to judge the critical speed in this study. With respect to excitations, three types of excitations are applied widely to railway vehicle system [15], and type (b) is chosen to calculate the critical speed here:(a)
No excitation, running on an ideal track, starting from the limit cycle, and reducing the speed until a stable wheelset motion is achieved(b)
Excitation by a singular irregularity, followed by an ideal track (or with a short irregularity sequence followed by an ideal track), with or without variation in the excitation amplitude(c)
Excitation by a stochastic (measured) track irregularity, and applying the criteria used during the vehicle acceptance test.


### 3.3. Hunting Frequency

The hunting frequency of the wheelset is an important factor in a vehicle dynamic system, whose value can significantly affect the occurrences of car body hunting and car body shaking phenomena. According to the Klingel theory [30], there is a direct relationship between the hunting frequency of the wheelset and the forward speed, as well as with the equivalent conicity of the wheel tread for a simple wheelset with coned wheels. The author’s model is completely different and complex, and the influences of other components and nonlinearities on the hunting frequency of the wheelset are worth investigating.In this study, the calculation of hunting frequency is divided into two parts. When the vehicle speed is less thanVcr, the system vibration will result in a stable equilibrium position and the nonlinearities will have no influences on hunting frequency. For a linear system, the hunting frequency can be determined using the roots of the characteristic equation, whose imaginary part is the angular frequency of the wheelset hunting motion. When the vehicle speed exceedsVcr, the system vibration will diverge to a limit cycle with a large amplitude and the nonlinearities of the system should be taken into consideration. In this way, a numerical method based on the characteristics of time histories of vehicle components is proposed to calculate the hunting frequency of the wheelset.The nonlinear governing equations of motion of the vehicle can be solved using the fourth-order Runge–Kutta method. The linear interpolation is subsequently used to change the variable step time series into a fixed step. The wheelset motion signals are collected over a finite timeTc and consist of a discrete number of points obtained at a selected sampling frequency. The data can then be modeled as a sum of the sine and cosine functions of timet with the period being an integer submultiple ofTc. The fast Fourier transform (FFT) is then used to calculate the frequency spectrum, and the hunting frequency can thus be obtained by selecting the main frequency in the spectrum. During the calculation of hunting frequency, the track irregularity is removed and the wheelset is set to an initial displacement.


## 3.1. Hopf Bifurcation Characteristics

A bifurcation diagram is a powerful tool for analyzing railway nonlinear dynamics. The railway system expresses the typical Hopf bifurcation, whose diagram is illustrated in Figure5 [17]. In Figure 5, the abscissa represents the vehicle speed while the ordinate represents the amplitude of the limit cycle. The solid and dashed lines denote the stable and unstable limit cycles, respectively. PointA is the Hopf bifurcation point, and the corresponding speedVlin at pointA is defined as the linear critical speed of railway vehicles. PointB is the turning point from a stable equilibrium point to a periodic solution, and the corresponding speed Vnl at point B is defined as the nonlinear critical speed of railway vehicles.Figure 5
Typical types of Hopf bifurcations of a railway vehicle system.For the subcritical Hopf bifurcation, system vibration will result in a stable equilibrium position for any external excitations when the vehicle speedV is less thanVnl. WhenV is larger thanVlin, the system will exhibit a stable periodic motion with a large amplitude regardless of the amplitude of the disturbance. In the intervalVnl <V <Vlin, however, the vibration of the system depends on the amplitude of the external excitation, and the interval is an uncertain region. For the supercritical Hopf bifurcation (b), the system always exhibits a stable periodic motion whenV achieves the linear critical speedVlin. With the increase inV, the amplitude of the hunting motion increases gradually. The supercritical Hopf bifurcation (a) has both subcritical Hopf bifurcation and the supercritical Hopf bifurcation characteristics. The uncertain region ofV extends toVnl ~Vc, in which the solution will converge to equilibrium or a limit cycle with a small amplitude at a small track excitation and diverge to a limit cycle with a large amplitude at a large track excitation. The bifurcation diagram can be calculated using a path following method or a set of numerical simulations [36–38]. The path following method is used in this paper.


## 3.2. Critical Hunting Speed

To investigate the stability of a railway vehicle system in practical engineering, the nonlinear critical hunting speedVnl is generally considered. From a mechanical viewpoint, a system can be regarded as stable if the oscillation of the system decays after a discontinuation excitation. Increasing the vehicle speed slowly, when the system displays a periodic motion at a constant amplitude, the corresponding speed can be regarded as the critical speedVcr, namely,Vnl. The lateral displacement of the front wheelset is used as the output value to judge the critical speed in this study. With respect to excitations, three types of excitations are applied widely to railway vehicle system [15], and type (b) is chosen to calculate the critical speed here:(a)
No excitation, running on an ideal track, starting from the limit cycle, and reducing the speed until a stable wheelset motion is achieved(b)
Excitation by a singular irregularity, followed by an ideal track (or with a short irregularity sequence followed by an ideal track), with or without variation in the excitation amplitude(c)
Excitation by a stochastic (measured) track irregularity, and applying the criteria used during the vehicle acceptance test.


## 3.3. Hunting Frequency

The hunting frequency of the wheelset is an important factor in a vehicle dynamic system, whose value can significantly affect the occurrences of car body hunting and car body shaking phenomena. According to the Klingel theory [30], there is a direct relationship between the hunting frequency of the wheelset and the forward speed, as well as with the equivalent conicity of the wheel tread for a simple wheelset with coned wheels. The author’s model is completely different and complex, and the influences of other components and nonlinearities on the hunting frequency of the wheelset are worth investigating.In this study, the calculation of hunting frequency is divided into two parts. When the vehicle speed is less thanVcr, the system vibration will result in a stable equilibrium position and the nonlinearities will have no influences on hunting frequency. For a linear system, the hunting frequency can be determined using the roots of the characteristic equation, whose imaginary part is the angular frequency of the wheelset hunting motion. When the vehicle speed exceedsVcr, the system vibration will diverge to a limit cycle with a large amplitude and the nonlinearities of the system should be taken into consideration. In this way, a numerical method based on the characteristics of time histories of vehicle components is proposed to calculate the hunting frequency of the wheelset.The nonlinear governing equations of motion of the vehicle can be solved using the fourth-order Runge–Kutta method. The linear interpolation is subsequently used to change the variable step time series into a fixed step. The wheelset motion signals are collected over a finite timeTc and consist of a discrete number of points obtained at a selected sampling frequency. The data can then be modeled as a sum of the sine and cosine functions of timet with the period being an integer submultiple ofTc. The fast Fourier transform (FFT) is then used to calculate the frequency spectrum, and the hunting frequency can thus be obtained by selecting the main frequency in the spectrum. During the calculation of hunting frequency, the track irregularity is removed and the wheelset is set to an initial displacement.


## 4. Results and Discussion

Simulation studies of self-excitation vibrations and stability of the nonlinear vehicle model were carried out based on the fourth-order Runge–Kutta method [39]. The results for the lateral displacement of the wheelset in this study are based on the first bogie leading wheelset.To determine the critical speed, the Dichotomy method was applied to the simulation. The simulation was started with a speed of 300 km/h, and the speed was increased or decreased in intervals of 10 km/h to determine another boundary of the interval. The calculation precision was 1 km/h, and the time series of the lateral displacement and the state-space plot of the wheelset lateral shift velocity versus displacement were determined. In the state-space plot, the data of the initial 3 s are removed. The results obtained for the last two speeds are shown in Figure6. It can be seen that once the vehicle speed is less than 322 km/h, for instance, 321 km/h, the lateral displacements of the vehicle components tend toward the equilibrium point after the track irregularities disappear. When the vehicle speed reaches 322 km/h, the lateral displacements of the vehicle components diverge to a stable limit cycle and a sustained lateral hunting motion occurs. Therefore, it can be concluded that the critical hunting speed of this vehicle is 321 km/h.Figure 6
Calculated wheelset lateral movements. (a) Time history at 321 km/h. (b) Shift velocity versus displacement at 321 km/h. (c) Time history at 322 km/h. (d) Shift velocity versus displacement at 322 km/h.
(a) (b) (c) (d)The bifurcation diagram of this vehicle is shown in Figure7(a). It is a typical supercritical Hopf bifurcation of type (a), and the linear critical speed is 440 km/h while the nonlinear critical speed is 321 km/h. PointP is on the unstable limit cycle with an amplitude of 3 mm, and the corresponding speed is 375 km/h. It means that, under the initial excitation of 3 mm, when the vehicle speed is less than 375 km/h the solution will converge to equilibrium, and when the vehicle speed is greater than 375 km/h the solution will diverge to a limit cycle.Figure 7
Calculation results. (a) Bifurcation diagram. (b) Hunting frequency.
(a) (b)Figure7(b) presents the variation in the hunting frequency of the wheelset with increasing vehicle speed. The initial wheelset amplitude is 3 mm. It is understandable that the hunting frequency increases with the increase in speed. However, there is a jump from pointQ1 to pointQ2 at 375 km/h which just corresponds to pointP in Figure 7(a). That is to say, the hunting frequency exhibits a sudden change when the solution of the movement changes from equilibrium to the limit cycle. It is partly caused by the discontinuity of the damping coefficient of the yaw damper, and the evidence can be seen in Figure 8, which shows the time histories of the wheelset displacement and yaw damper velocityvd. The red line represents the blow-off velocity of the yaw damper. It can be seen that, whenvd is lower than the blow-off velocity, the cycle of the wheelset displacement is 0.25 s, namely, a frequency of 4 Hz. Whenvd is greater than the blow-off velocity, the cycle of the wheelset displacement is 0.19 s, namely, a frequency of 5.3 Hz. The influence of the damping coefficient on the hunting frequency needs to be explained, and it is described below. Thus, it can be concluded that the hunting frequency exhibits a sudden increase once the blow-off of the yaw damper occurs, and the hunting frequency can be divided into two parts: before blow-off and after blow-off.Figure 8
Time histories of (a) wheelset lateral displacement and (b) yaw damper velocity.
(a) (b)To illustrate the consequences of this phenomenon, the time-frequency analysis for a section of the measured accelerations on the car body floor is given in Figure9(a). It can be seen that the car body acceleration around 10 Hz exists for almost the whole-time history, which corresponds to the elastic mode of the car body, while the acceleration around 1 Hz is caused by the rigid mode of the car body. The car body acceleration during the period of 39220–39240 s, whose frequency is around 10Hz, is significantly larger than other sections. It is a typical phenomenon of car body shaking [27], and the local enlargement of the acceleration in the resonance section is shown in Figure 9(b). It can be seen that there exists a harmonic vibration in the resonance section, and the acceleration amplitude increases dramatically. In general, the hunting frequency of the bogie can hardly reach 10 Hz, only if the yaw damper blows off and a large jump occurs in the hunting frequency. The situation outlined above exists in daily operation when the equivalent conicity reaches a high level and the vehicle loses stability.Figure 9
Measured signal in time and frequency domains. (a) STFT spectrum of lateral acceleration on the back floor. (b) Time history of local lateral acceleration of vehicle floor.
(a) (b)Predictably, the choice of key parameters largely determines the critical speed of the vehicle. Furthermore, the resonance problem such as those pertaining to the hunting and shaking of the car body, which are closely related to the hunting frequency, should be focused on during the design of the vehicle. For this reason, the combined effects of key parameters in the vehicle system on critical speed and hunting frequency are analyzed immediately. These parameters include the first stage dampingCd1, the second stage damping Cd2, the series stiffness of the yaw damperKd, the coefficient of frictionu, longitudinal stiffness of primary suspensionKpx, creep force coefficientsf11 andf22, and the lateral damping of secondary suspensionCsy.Figure10 shows the influences of yaw damper parameters on the critical speed and hunting frequency. During the calculation of hunting frequency, the vehicle speed is taken as 300 km/h. In the figures, the black solid line represents the critical speedVcr, the blue solid line represents the hunting frequencyf, and the blue dash line represents the jump of the hunting frequency. Figure 10(a) depicts the results of Cd1. It can be seen thatVcr increases logarithmically with the increase in Cd1 and finally stabilizes at 330 km/h. Thef exhibits a sudden change at the value of 140 kN·s/m whereVcr reaches 300 km/h. Before that,f maintains a constant value, and once Cd1 exceeds 140 kN·s/m,f decreases with the increase in Cd1. We can thus draw the conclusion that Cd1 only affects the hunting frequency before blow-off of the yaw damper. Figure 10(b) illustrates the results of Cd2. It can be seen thatVcr virtually increases linearly in two parts with the increase in Cd2. The dividing point of these two parts is located at the value of 15 kN·s/m. Similar to Cd1,f also exhibits a sudden change at the value of 15 kN·s/m whereVcr reaches 300 km/h. The difference is thatf decreases with the increase in Cd2 before the sudden change, and, once Cd2 exceeds 15 kN·s/m,f maintains a constant value. It indicates that Cd2 only affects the hunting frequency after blow-off of the yaw damper. Furthermore, it can be seen from Figure 10(c) thatVcr firstly increases with increasingKdand then decreases and finally tends to a stable value. Thef decreases with the increase inKd overall, and the jump off also exists at the critical speed of 300 km/h.Figure 10
Change in critical speed and hunting frequency with (a) increase ofCd1, (b) increase ofCd2, and (c) increase ofKd.
(a) (b) (c)Figure11 shows the change ofVcr andf with the increase ofu,f11,Kpx, andCsy. During the calculation off, the vehicle speed is taken as 350 km/h. The lateral creep force coefficientf22 is set equal tof11. They both affectVcr andf by changing the creep forces. TheVcr decreases with the increase inu, as shown in Figure 11(a). In addition,f increases with the increase inu whenVcr is less than 350 km/h, andf maintains a constant value whenVcr is larger than 350 km/h. It indicates thatu affects the hunting frequency only if the vehicle loses stability. Figure 11(b) shows the results off11 (f22), and it can be seen that, with the increasing off11 (f22),Vcr decreases gradually,f increases overall, and the jump off exists at the critical speed of 350 km/h. Apart from the yaw damper parameters,Kpx andCsy are another two important parameters in the vehicle suspension system.Figure 11
Change in critical speed and hunting frequency with (a) increase ofu, (b) increase off11, (c) increase ofKpx, and (d) increase ofCsy.
(a) (b) (c) (d)Figure11(c) illustrates the influence ofKpx onVcr andf. The results show thatVcr increases logarithmically with the increase inKpx and finally stabilizes at 322 km/h. Thef decreases gradually with the increase inKpx and finally stabilizes at 4.9 Hz. As forCsy, the influences on bothVcr andf are coincident with the results ofKpx in the trend, as shown in Figure 11(d).To further compare the influence of the same parameter whose value is taken in different sections on the critical speed and hunting frequency, the parameters analyzed above were divided into three sections, as shown in Table1. Within the range of parameters, the change of critical speed is defined as ΔVcr, while the change of hunting frequency is defined as Δf. The comparison results of different parameters and different sections with regard to ΔVcr are shown in Figure 12(a). It can be seen that parameters ofu, Cd2, Kd,f11, and Cd1 have a significant influence on the critical speed. Within the range of these parameters, ΔVcr all exceed 100 km/h. The parameters ofKpx andCsy, by contrast, have less influence on the critical speed, and ΔVcr are within 50 km/h. As for the proportion of each section, the first section plays a decisive role inCd1,u,Kpx, andCsy. The first section is close to the second section inf11, and the third section accounts for a small proportion in all parameters. Figure 12(b) indicates the results of Δf, and the influence of frequency jump is removed during the calculation. It can be seen that the parameters ofu and Cd2 process the greatest impact on the hunting frequency, and the parameters of Cd1,Kd, andf11 come next. However, it is worth noting that the parameters ofu andCd2 only work when the vehicle loses stability, and parameters of Cd1 only work when the vehicle does not lose stability. The parameters ofKpx andCsy, by contrast, have less influence on the hunting frequency, and Δf are within 0.8 Hz. As for the proportion of each section, all of the parameters show the same regularity; that is, the first section comes first, the second section comes next, and the third section is the least.Table 1
Value division of analyzed parameters.
Parameter First section Second section Third section Unit C d 1 20–300 300–600 600~1000 kN·s/m K d 0.1–6 6–12 12–20 MN/m C d 2 0–40 40–80 80–120 kN·s/m u 0.05–0.2 0.2–0.35 0.35–0.5 - K p x 10–30 30–50 50–80 MN/m f 11 2–5 5–8 8–11 MN C s y 0–20 20–40 40–60 kN·s/mFigure 12
Comparison of the influences of different parameters and different sections on (a) critical speed and (b) hunting frequency.
(a) (b)The parameters analyzed in the previous section significantly influenceVcr without changing the bifurcation type. However, the blow-off forceFu affects not onlyVcr but also the bifurcation type, as shown in Figure 13. With the increase inFu, the type of Hopf bifurcation transforms from subcritical to supercritical Hopf bifurcation of type (a), and, whenFu tends to infinity, namely, with a linear damping coefficient, the bifurcation type transforms to supercritical Hopf bifurcation of type (b). Furthermore,Fu does not affect the linear critical speed, which is decided by Cd1, and increasingFu can increaseVcr until the system reaches the linear critical speed.Figure 13
Variation of bifurcation type with increase inFu.


## 5. Conclusions

In this study, the lateral dynamic model of a vehicle with 17 degrees of freedom was established, and the nonlinearities of wheel/rail contact geometry, creep forces, and yaw damper were considered. The wheel profile type of S1002CN and the rail profile type of UIC 60 were selected, and the bifurcation type was supercritical Hopf bifurcation of type (a). In addition to the critical speed, the hunting frequency of the wheelset was considered together during the calculation. The following conclusions can be drawn:(1)
Owing to the discontinuity of the damping coefficient of the yaw damper, the hunting frequency of the wheelset increases suddenly when blow-off of the yaw damper occurs.(2)
Parametersu, Cd2, andKd have the greatest influence on the critical speed, followed by parameters Cd1 andf11, and parametersKpx andCsy have the least influence. Parametersu and Cd2 have a significant influence on the hunting frequency; however, they only work when the vehicle loses stability. Conversely, Cd1 only affects the hunting frequency when the vehicle does not lose stability.(3)
As for the influence proportion of each section on the critical speed and hunting frequency, all of the parameters analyzed show the same regularity that the first section comes first, the second section comes next, and the third section is the least.(4)
Fu affects not only the critical speed but also the bifurcation type. With the increase in the blow-off force, the Hopf bifurcation transforms from subcritical to supercritical of type (a), and then to supercritical of type (b).It is noted that the results presented herein are based on calculation performed considering the straight line, and further research regarding the influence of vehicle parameters on the stability when a vehicle traverses a curve is needed. This will be a topic in our future work.



---

*Source: 7163732-2019-04-21.xml*