Establishment of Dynamic Model for Axle Box Bearing of High-Speed Trains Under Variable Speed Conditions

In this study, a dynamic model for the bearing rotor system of a high-speed train under variable speed conditions is established. In contrast to previous studies, the contact stress is simplified in the proposed model and the compensation balance excitation caused by the rotor mass eccentricity considered. The angle iteration method is used to overcome the challenge posed by the inability to determine the roller space position during bearing rotation. The simulation results show that the model accurately describes the dynamics of bearings under varying speed profiles that contain acceleration, deceleration, and speed oscillation stages. The order ratio spectrum of the bearing vibration signal indicates that both the single and multiple frequencies in the simulation results are consistent with the theoretical results. Experiments on bearings with outer and inner ring faults under various operating conditions are performed to verify the developed model.


Introduction
In recent years, the rapid development of high-speed trains has led to a rapid increase in the mileage and number of high-speed trains globally. The safe and reliable operation of trains is critical to railway safety [1,2].
As key components of train running systems, the movement modelling and state monitoring of bogie axle box bearings have attracted extensive research attention [3]. Dynamic simulations of bearing systems are commonly used to better understand the mechanism of bearing faults and their effects on movement. In these bearing models, a set of differential equations is established to study the influence of parameters on bearing motion states. These equations can be solved using numerical solution methods.
In the early models, relatively simple nonlinear equations were used to simulate the bearing motion and study the motion law of bearings based on the exact solutions of these equations. For instance, Chinta et al. [4] used dimensionless differential equations to describe the unbalanced response of rotor bearings and Floquet theory to study the stability and bifurcation behavior of bearing systems. Sinou [5] established a dynamic model of a flexible rotor-rolling bearing system under unbalanced excitation, used the harmonic balance method to study the nonlinear motion law of the system, and analyzed the influence of the bearing radial clearance and unbalanced mass on the Hertzian contact force between the roller and raceway.
Because traditional methods for solving nonlinear equations can only generate approximate solutions for nonlinear systems with fewer degrees of freedom (DOF), more advanced methods have been developed for higher-DOF bearing systems. Fukata et al. [6] used a numerical method to solve the differential equations for ball bearing motion for the first time and intuitively analyzed the

Open Access
Chinese Journal of Mechanical Engineering motion state of the bearings. Tiwari et al. [7] established a nonlinear dynamic model for bearing rotor systems and analyzed the influence of radial clearance on motion stability based on Floquet theory. Sawalhi et al. [8] built a dynamic model of a bearing and gear coupling system in which a small sprung mass with relatively high damping is used to represent a typical high frequency. The model was verified experimentally under different fault conditions [9]. Tsuha et al. [10] developed a rolling bearing model with improved fidelity by introducing a set of equivalent contact stiffness and damping values and using the theory of elastic fluid lubrication to calculate the contact force. Song et al. [11] proposed a dynamic model for spindle bearing systems by combining the angular contact and floating bearings. Compared with traditional methods, more factors that affect the bearing motion are considered in these methods. The established bearing motion models are therefore more accurate. At the same time, the numerical solution methods provide a more intuitive solution of the equations as well as greater convenience for further nonlinear analysis. However, these previous works have only focused on bearings under constant operating conditions, which cannot adequately describe dynamic operating conditions. In contrast to bearings operating under relatively consistent operating conditions, high-speed train bearings in running systems often operate under high-speed, heavy-duty, and harsh conditions and are severely affected by strong noises in the surrounding environment. It is difficult to evaluate the bearing quality directly and study the bearing motion states. To address these challenges, dynamic model simulation methods are required for high-speed train bearings. The influence of various factors on the motion states of bearings should be considered in the dynamic modeling method to reduce the research cost and improve research efficiency. Cao et al. [12] established a nonlinear dynamic model for the CRH1 EMU (China Railway High-speed-1, Electric Multiple Units) coupling system in which the loose between the inner ring and main shaft of the axle box bearing in a high-speed train is considered. The motion characteristics of the model were also studied and compared using optimized numerical algorithms under various working conditions. Cao et al. [13] established a dynamic model for high-speed rolling bearings by considering actual factors such as the centrifugal and thermal expansion of the inner ring. The model was then used to predict system parameters such as the contact angle and bearing stiffness and to perform a correlative analysis of the bearing damage mechanism. Yang et al. [14] established a 4 DOF nonlinear dynamic model with classic faults for high-speed train bearing rotor systems. The model was verified using a certain type of high-speed train bearing under outer ring, inner ring, and roller element fault conditions. Liu et al. [15] established a 12 DOF nonlinear dynamic model for high-speed train axle box bearings with outer ring faults. The influence of the bearing speed and fault size on the system motion was studied from a nonlinear perspective. Wang et al. [16] proposed a novel stochastic vehicle-track coupled model to evaluate the dynamic performance of axle box bearings in a high-speed train with unsteady wind loads and random track irregularities. The model was used to study the influence of crosswind speed on the dynamic performance of the bearings.
Despite these achievements in bearing model research, the existing works face the limitation that the failure mechanism and bearing motion were all studied under static bearing speeds. Mishra et al. [17] established a bond graph model of rolling element bearings under special operating conditions and experimentally verified its effectiveness [18]. The main function of this model is to generate a series of nonstationary signals for fault diagnosis algorithm research. However, this model cannot describe the vibration response of each part of the bearing rotor system. Moreover, the previous bearing models have a single structure, which is insufficient to represent the structural characteristics of the tapered roller bearings in high-speed trains. Furthermore, in addition to the bearing mass eccentricity, the shaft mass eccentricity is also very important but has always been ignored in previous studies on bearing rotor systems. To address these limitations, in this study, a dynamic model for the coupling system of bearings and shafts is established to accurately simulate the dynamics of an axle box bearing in a high-speed train under variable speed conditions. The novelties of this study are threefold: First, the shaft and two bearings are regarded as a coupled system in which the shaft is regarded as a rotor and the mass eccentricity of the shaft is considered. Second, the angular iteration method is used to calculate the spatial position of the bearing rollers at each moment under the variable-speed operation of the bearing. Third, the model is verified through a series of simulations and actual experiments with different fault modes.
The remainder of this paper is organized as follows: Section 2 describes the process of building the shaftbearing coupling system model under variable speed conditions and introduces the outer ring fault, inner ring fault, and rolling element fault models. Section 3 presents the simulation results when different types of faults are present in the model under variable speed conditions. Section 4 describes the verification of the model effectiveness through a comparison of the physical and simulation experiments.

Bearing Dynamic Model
The running part of a high-speed train consists of an axle, a left wheel, and a right wheel. According to the law of energy conservation, the axle mass is equivalent to the axle center of mass. Figure 1 shows the structure of the axle box bearing in a high-speed train. O 1 , O 2 and O 3 are the bearing geometry center, shaft section center, and shaft section center of mass, respectively; e is the eccentricity of the shaft cross section; m c is the equivalent mass (kg) of the shaft at the center; K and C are the stiffness and damping of the shaft, respectively; F xL and F yL are the reaction forces (N) received by the left end bearing along the horizontal and vertical directions, respectively; F xR and F yR are the reaction forces (N) received by the right end bearing along the horizontal and vertical directions, respectively; and α and β are the contact and half cone angles in degrees, respectively. Because the rollers of tapered roller bearings change the direction of the contact force and decompose the load on the bearing into the axial and radial components, the half-cone angle should be considered in the model. Figure 2 shows the schematic diagram of the bearing dynamics model. m 1 , m 2 and m b are the masses (kg) of the bearing inner ring, bearing, and unit resonator, respectively; x 1 and y 1 are the displacements (m) of the inner ring along the horizontal and vertical directions, respectively; K 1 and C 1 are the support stiffness (N/m) and support damping (N·s/m) of the inner ring, respectively; x 2 and y 2 are the displacements (m) of the outer ring along the horizontal and vertical directions, respectively; K 2 and C 2 are the support stiffness (N/m) and support damping (N·s/m) of the outer ring, respectively; and y b , K b and C b are the vertical displacement (m), stiffness (N/m), and damping (N·s/m) of the unit resonator, respectively. The bearing model in Figure 2 is coupled with a series of spring mass models and simulated by establishing a set of differential equations according to Newton's second law.
When a fault occurs, the periodic impact due to the fault leads to the natural vibration of the bearing inner and outer rings and other components. These vibrations can be accurately simulated by adjusting the stiffness and damping coefficients of the unit resonator [19,20]. According to Newton's second law, the vibration of the shaft along the horizontal and vertical directions can be described by: The vibration equations for the right-end bearing along the horizontal and vertical directions are expressed as follows: The vibrations of the left-end bearing along the horizontal and vertical directions can be calculated as follows: The vibration equations for the unit resonator along the vertical direction at the right and left bearings can be expressed as: where x c and y c are the displacements of the shaft centroid along the horizontal and vertical directions, respectively; x r1 and y r2 are the displacements of the outer and inner rings at the right end bearing along the horizontal direction, respectively; y r1 and y r2 are the displacements of the outer and inner rings at the right end bearing along the vertical direction, respectively; x l1 , x l2 , y l1 , and y l2 are the corresponding displacements at the left end bearing; y rb and y lb are the displacements of the unit resonator at the right and left ends along the vertical direction, respectively; and F is the axle load of the train.

Bearing Support Reaction Force
The bearing support reaction force is required in this model. However, the magnitude and direction of the force changes with the rotation of the bearing. The force is calculated as follows: • The contact force at each roller is calculated.
• The force is decomposed along the horizontal and vertical directions of the bearing.
• The bearing support reaction forces along the two directions are obtained by combining all the forces in the respective directions. Because the forces on the inner and outer rings of the bearings are transmitted through the rollers during bearing movement, the reaction force is the sum of the resultant forces of the contact forces at each roller. In practice, the motion state and force of the bearings are very complicated. To simplify the calculation, the following assumptions were made. These assumptions are commonly used and do not reduce the accuracy or fidelity of the model. • The outer ring of the bearings is fixed on a rigid element with a rotating velocity of 0 and exhibits horizontal and vertical displacements. • The rolling elements are equidistantly located on the raceway and perform pure rolling (i.e., sliding of the rolling elements is not considered). • The contact stress is modelled using Hooke's law. Figure 3 shows the side view of a double-row tapered roller bearing with N 0 rollers equally separated by an angle of 2π/N 0 . During rotating operation, the angle that the i-th roller (i = 1, 2, …, N 0 ) turns at time t can be calculated as: where γ(t) is the angle through which the roller indexed by 1 turns at time t, and its initial value is γ(t 0 ) = 0. The normal contact deformation of the i-th roller of the rightend bearing with the raceway at angular displacement θ i can be calculated as: where c 0 denotes the bearing radial clearance. The contact force between the ball and the raceway can be expressed as a nonlinear Hertz contact force [21,22]: where K t is the contact stiffness and n = 10/9. When δ i > 0, H i = 1, which indicates that there is a nonlinear Hertz contact force, whereas when δ i < 0, H i = 0, which indicates that there is no nonlinear Hertz contact force. By decomposing the contact force at every roller along the horizontal and vertical directions, the resultant force in the two directions can be obtained. The total contact force of the bearing along the horizontal and vertical directions is the sum of the resultant forces from all the rollers, which is given by: where F xR and F yR are the contact forces along the horizontal and vertical directions, respectively. The corresponding F xL and F yL for the left-end bearing can be similarly obtained.

Model Analysis with Different Fault Modes
To solve the model, the angle γ(t) through which the roller rotates at any time t is required. However, in a high-speed train, the angular rotating velocity of the main shaft, w 0 (t), is not constant. To address this, we simplified and applied the angle iteration method [17] to calculate the angle γ(t) in real time. The angle γ(t) can be calculated as: When the rotational speed changes with time, γ(t) cannot be expressed as a linear function of time, and the total angle through which the roller turns at each moment cannot be determined from the angular velocity.
In the proposed angle iteration method, it is assumed that the bearing speed is piecewise constant, that is, the speed is constant within a short period of time Δt (Δt → 0). Consequently, the angle through which the bearing rotates at time (t + Δt) can be expressed as γ(t + Δt) = γ(t) + w(t)Δt where w(t) is the relative speed of the roller with respect to the location of the fault, which can be on the outer ring, inner ring, or roller. Therefore, w(t) has different values for different fault types. When the fault is on the outer (or inner) ring, w(t) is the speed of the roller relative to the outer (or inner) ring at t. When the fault is on the rolling element, w(t) is the rotation speed of the roller at t. A schematic of this process is shown in Figure 3. The problem is then transformed into the estimation of w(t). In the simulation, Δt represents the step size used in each iteration of the ode45 solver. Figure 4 shows a scenario in which the outer ring has a fault [23,24] located in the lower bearing area where θ c is the location of the fault center, θ ca is the central angle of the arc on which the fault is located, L 0 is the fault width, r is the inner ring radius of the bearing, R is the outer ring radius of the bearing, and r g is the section radius at the center of mass of the tapered roller. When a roller rolls over the fault area, the roller descending depth c outer between the roller and raceway can be expressed as:

Modeling of Outer Ring Fault
where θ ca = 2 arcsin( L 0 2R ). Because the outer ring speed is zero, the relative speed of the roller with respect to the outer ring w outer (t) can be expressed as: where d is the roller diameter, and D is the bearing pitch diameter. Figure 5 shows a scenario in which the inner ring of the bearing has a fault. In contrast to the outer ring fault, the location of the inner ring fault continuously changes as the bearing rotates. Because there is relative movement between the bearing cage and inner ring, the inner ring is used as the reference system. This makes the form of the inner ring fault equivalent to that of the outer ring fault.

Modeling of Inner Ring Fault
(10) When a roller rolls over the fault area on the inner raceway, the roller descending depth c inner between the roller and raceway can be expressed as: where θ ca = 2 arcsin( L 0 2r ). The relative speed of the roller with respect to the inner ring w inner (t) can then be calculated as: Figure 6 shows a scenario in which the rolling element of the bearing has a fault and θ rc is the central angle of the arc on which the fault is located. When the bearing has a rolling element fault, the roller contacts the inner ring raceway and the outer ring raceway once in its own rotation cycle as the inner ring of the bearing rotates. Here, we set the bearing cage as the reference system and assume that the initial angle of the fault area is 0° and that the rollers only rotate (without revolution). When the roller rotates 90° counterclockwise, it contacts the inner ring. When the roller rotates 270° counterclockwise, it contacts the outer ring. During this process, the roller descending depth c roller when the roller rolls over the fault area can be calculated as:

Modeling of Rolling Element Fault
where φ = π/2 and 3π/2 when the inner and outer raceways are contacted, respectively; θ rc = 2 arcsin( L 0 2r 0 ) ; and the rotating speed of the roller w roller (t) is given by: Based on the above analysis, the contact deformation in the fault area under different fault modes can be summarized as follows: where c change is c outer , c inner , or c roller depending on the fault mode. (15)

System Parameters
A series of simulations were performed to verify the proposed model under varying speeds. Tables 1, 2, and 3 list the main parameters of the high-speed EMU axle box bearing, rotor system, and unit resonator used in the simulation, respectively.

Introduction to GATD
To verify the model, the simulation results are compared with the theoretical results. The model was first simulated under variable speed conditions to obtain the time-domain waveform of the bearing vibration acceleration. The order ratio spectrum of the vibration was then calculated for bearing diagnosis. The model is verified to be correct if the order ratio spectrum of the vibration obtained from the simulation is close to the theoretical one. Therefore, the key step is to find a suitable fault diagnosis method to obtain the order ratio spectrum of the bearing vibration acceleration. In 2016, Urbanek et al. [25] proposed a generalized angular temporal deterministic (GATD) signal model for rotating machinery failures under variable speed conditions based on cyclostationary theory. The GATD signal extracted from wind turbine bearings enables the accurate identification of bearing faults. The advantage of the GATD-based method is that it does not rely on angle resampling and can be used to perform bearing fault diagnosis while retaining the original signal characteristics. Therefore, this method is adopted in this study to extract the fault characteristics from the model simulation results. The GATD-based fault diagnosis process is as follows: (1) First, the collected faulty bearing signal x(t) is normalized according to the Z-score standardization method: where φ(t) is the angle through which the roller turns at time t with the initial phase φ 0 , 1 T = max dϕ dx , w T (t) is the Hanning window function, and τ(φ) is the angle-fixed time increment used for positioning the window. μ T (φ) is the local angular-temporal mean value of the signal x(t), and σ T (φ) is the localized angular-temporal standard deviation of x(t), which is calculated as: where λ = t -τ(φ+φ 0 ).
(2) A Fourier transform is performed on the normalized signal: (3) A fast Fourier transform is performed on the squared envelope spectrum of Z T (φ, f; φ 0 ) to obtain the angular-temporal spectrum (ATS) of the signal: where Z T denotes Z T (φ, t; φ 0 ) to simplify the notation, Ω is the frequency related to angle-fixed events (expressed in orders), f denotes the event occurrence frequency in Hz, and Θ is the maximum angle period. The obtained ATS is a family of curves distributed in three-dimensional (3-D) space. The order ratio spectrum of the signal can be obtained by plotting ATS T with Ω as the horizontal axis.

Model Simulation
According to theoretical mechanics, the fault characteristic order ratio of the outer ring, inner ring, and rolling element can be calculated as [ [26]]: For the bearing used in this study, f o = 7.082, f i = 9.917, and f r = 2.871.
The simulation was conducted using MATLAB at a sampling frequency of 512000 Hz. The angular rotation velocity of the main shaft w 0 (t) given by Eq. (25) is shown (22) in Figure 7. This speed was chosen to simulate the entire operating process of the bearing system, which includes accelerating, decelerating, and oscillating processes with a sine law. Three simulations were conducted with the fault size of L 0 = 1 mm for the three types of faults, i.e., an outer ring fault, an inner ring fault, and a rolling element fault.
The order ratio spectrum analysis method in Ref. [25] was used to obtain the corresponding order ratio spectrum from the solutions of Eqs. (1)-(4). The timedomain waveform of the bearing acceleration and the corresponding order ratio spectrum were obtained by simulating the model. Figures 8, 9, 10 show the results.
It can be intuitively seen from the above time-domain graphs that when a fault occurs, the acceleration amplitude of the bearing changes with the shaft angular speed. In addition, the maximum amplitude increases to different extents depending on the type of fault. This implies that the location of the fault has different effects on the vibration of the bearing. The bearing has the most severe vibration when a fault occurs on the rolling elements. Moreover, the results show that the simulated characteristic order ratios of the outer ring (7.085), inner ring (9.9172), and rolling element (2.8694) are close to the calculation results given by Eqs. (22)-(24). Simultaneously, the corresponding frequency multiplication can be observed in the order ratio spectrum, which proves the correctness of the model.
Because the order ratio spectrum is obtained from the signal after two Fourier transforms, there is a linear relationship between the order ratio spectrum amplitude and the original signal amplitude. Let the amplitudes corresponding to a single frequency (ACSFs) in the signal order ratio spectra under an outer ring fault and an inner (25) Figure 11. Figure 11 shows that the ACSF values under the outer ring fault are always larger than those under the inner ring fault, and that the ACSF values in the order-ratio spectrum under these two fault conditions increase as the fault size increases. The latter indicates that the fault size has an important influence on the motion state of the bearing. As the size of the fault increases, more energy is generated in the fault area when the roller collides with the raceway, which causes an increase in amplitude at fault characteristic frequency.

Experimental Verification
To make the operating conditions close to the actual ones, a speed profile in which the bearing rotating speed was increased from 0 to a certain speed, maintained for a period of time, and then decreased to 0 was used in the experiments. It is important to note that in the experiments, no tachometer was used to measure the speed of the bearing and no program was used to control its speed. For fault type verification, the vibrations were obtained from the experiment and the speed profile was not needed. However, the speed profile was used as the model input for the verification of the bearing operating condition. As the model speed profile was unknown, the instantaneous speed of the bearing must be estimated from the bearing vibration signals.

Fault Type Verification
Two EMU axle box bearings were used in the experiments. One bearing had a through-groove fault on the outer ring with a width of L 0 = 1 mm, and the other a pitting fault on the inner ring with L 0 = 0.1 mm. A schematic diagram of the outer ring, fault on the outer ring, and fault on the inner ring is shown in Figure 12.
In the experiments, the vibration signal of the bearing was collected for 60 s at the sampling frequencies of 51 and 200 Hz.
Acceleration sensors were installed near the experimental bearing so that the vibration signal could be accurately measured. Figure 13 shows the experimental platform and horizontal and vertical locations of the sensors.
The experiments were conducted on the bearings with the seeded faults. Figure 14 shows the time-domain signal and order ratio spectrum of the bearing vibration acceleration for the bearing with the outer ring fault. The differences between the experimental and calculation results for the first three-order fault feature order ratio are 5.97%, 5.67% and 6.34%, respectively. Figure 15 shows the time-domain signal and order ratio spectrum of the bearing vibration acceleration for the bearing with an inner ring fault. The differences between the experimental and calculation results for the first three-order fault feature order ratio are 2.59%, 3.95%, and 4.15%, respectively. As mentioned earlier, the bearing speed changed during the experiments. Under these nonlinear operating conditions, the acceleration amplitude of the bearing changed continuously with time, which is typical of variable-speed bearing operation in high-speed trains. The results show that compared to the theoretical values given by Eqs. (22) and (23), the differences in the fault characteristic order ratios under the two experimental conditions are within the allowable range, which proves the effectiveness of the model.

Bearing Operating Conditions Verification
To further prove that the model is also effective under actual working conditions and robust against interference from noise, the bearing speed estimated in the experiment was used as the input of the model, and noise was added to the bearing vibration acceleration signal obtained from the model. Order ratio spectrum analysis was performed on the model output vibration signal with noise. The model is verified to be correct if the model results are consistent with the calculation results.
Because no speed measuring device such as a tachometer was used and the speed profile was not known, the real-time speed of the bearing spindle could not be obtained. Therefore, the instantaneous speed curve was extracted from the signals collected in the two experiments and a time function was used to fit the two working conditions for the model speed input. Simultaneously, on-site noise was added to the model to simulate an actual noisy environment.

Verification Under Outer Ring Fault
The speed curve was extracted from the experimental signals using the method in Ref. [27], and the noise signal extracted using the wavelet noise reduction method [28]. Figure 16 shows the instantaneous speed and experimental site noise extracted from the outer ring fault test signal. The fitted curve function is f(t) = − 0.002t 4 + 0.2535t 3 − 11.7572t 2 +232.4658t -2.0618 × 10 −4 .
It can be observed that the bearing speed experienced a process of increase, slow change, and decrease during the experiment. The noise amplitude also underwent the same process. Figure 17(a) shows the vibration signal obtained from the model, fitted speed curve, and noise. The order ratio spectrum in Figure 17(b) shows that the frequency-singled, frequency-doubled, and frequencytripled values are close to the theoretical values given by Eq. (22). The results demonstrate the effectiveness of the model.

Verification Under Inner Ring Fault
The same methods were used to obtain the time-frequency curve and experimental site noise in the inner ring fault experiment. Figure 18 shows the instantaneous speed and experimental site noise extracted from the inner ring fault test signal. The fitted curve function is f(t) = − 0.0028t 4 + 0.3423t 3 − 14.3913t 2 + 251.987t − 6.0213×10 −4 . The model vibration signal shown in Figure 19(a) was obtained using the proposed model under the fitted speed profile shown in Figure 18(a). The order ratio spectrum in Figure 19(b) shows that the frequency-singled, frequency-doubled, and frequency-tripled values are close to the theoretical values given by Eq. (23). The results demonstrate the effectiveness of the model.
It is observed that the ACSF values under the outer ring fault are smaller. This is because when the bearing speed changed, the instantaneous impact load and amplitude caused by a single fault were reduced. There are also very few messy frequency components in the order ratio spectrum, which proves that the model has good antiinterference ability. Distinct results could be obtained from the simulation model even when experimental site noise was added.

Conclusions
(1) A dynamic model for the bearing system in the running part of high-speed trains under variable speed conditions was established. The angle iteration method is used in the model to address the problem that the roller space position cannot be determined during bearing rotation. This problem was ignored in previous studies. (2) By simulating the actual operating conditions of the train, the vibration response of the system during a process of speed increase, decrease, and oscillation was obtained. The results indicate that the vibration acceleration amplitude of the bearing outer ring is positively correlated with the spindle speed of the bearings. The results of the order spectrum analysis are also consistent with the theoretical results.