Abstract

A three-dimensional deformable finite element model of wheel-rail system is established. The vertical vibration natural frequencies are calculated by modal analysis. The improved central difference method by ABAQUS explicit nonlinear dynamics module is chosen to calculate the vertical vibration of wheel axle. Through the vibration measurement experiment of wheel by China Academy of Railway Sciences, the contact parameter is optimized and the model is modified. According to the nonlinear simulation results, the important influence of vertical track irregularity on the wheel set high-frequency vibration is discovered. The advantage frequencies through spectrum are close to the vertical vibration natural frequencies of seventh, eighth, ninth, and tenth mode. When the velocity is 350 km/h, system’s nonlinear dynamic characteristic is higher than the results at 200 km/h. The critical wavelength of vertical track irregularity on the wheel vertical vibration is about 0.8 m. In addition, a particular attention should be paid to one situation that the main dominant frequency amplitude is more than 50% of the acceleration amplitude even if the peak acceleration is low.

1. Introduction

With the development of China’s railway and urban rail transportation, the speed of the train is rapidly improving the safety and comfort issues have become more prominent. Many domestic and global scholars have carried on a lot of research [1]. In 1922, Timoshenko studied track dynamic by using elastic foundation beam model, which had become a classic method [2]. In 1972, Lyon and Mather established the basic model of wheel-rail dynamic analysis by considering the basic parameters of track and vehicle [3, 4]. In the model, track was described as continuous elastic foundation “Euler Beam” and vehicle is simplified to spring-mass system. Hertz nonlinear spring model was established to simulate the wheel-rail contact.

In 1987, Li gave finite element analysis formulation of the Derby basic model. Xu et al. made a simulation analysis on wheel-rail impact of track joints by using the Timoshenko beam model of continuous elastic foundation [5, 6]. In 1992, Zhai, a professor in Southwest Jiaotong University, combined the vehicle system with the track system as a coupled system for further study. This researcher established a series of vehicle-track interaction coupled model; the major dynamic factors of the vehicle-track structure had been taken into detailed consideration through the whole system [7, 8]. In 2000, Lei and Chen used finite element theory to establish a freedom vibration model of track structure; an analysis on the natural frequencies and vibration characteristics of existing railway lines in different support conditions had been made in the researcher’s paper [9]. Zhang and Dhanasekar examined the effect of braking/traction torque to the longitudinal and lateral dynamics of railway wagons in 2009 [10]; this paper dealt with contribution of defective tracks to vehicle safely. Then in 2013, he presented a computational method for eliminating severe stress concentration at the unsupported railhead ends in rail joints through innovative shape optimization of the contact zone. With a view to minimizing the computational efforts, hybrid genetic algorithm method coupled with parametric finite element has been developed and compared with the traditional genetic algorithm (GA) [11].

At present, there are not many researches on high-frequency nonlinear vibration characteristics in the field of high-speed trains. Over the past 10 years, many scholars had studied the wheel-rail coupling dynamics, which mainly used the method of numerical calculation. Numerical calculation method, such as finite element method, can simulate the coupling situation with a good applicability and powerful analysis ability [1214]. However, the time step must be short when the characteristic was calculated at high-speed condition. The workload could be quite large by using the iteration method. In addition to numerical calculation method, some current researches are carried out by using analytical method [1517].

Track irregularity is an important source of vehicle vibration. It is one of the core issues of high-speed railway to realize the high track smoothness. Therefore, it is very important to control the influence of track irregularity on the dynamic performance of high speed. The short vertical track irregularity mainly causes the high-frequency vibration, whose wavelength is less than 2 m in the small range.

This paper presents a solution for inherent characteristic of wheel system of type. Meanwhile, FE dynamic model of the wheel-rail system using the contact dynamics theory is established. The wheel system’s high-frequency nonlinear vibration response is presented, which is in the influence of vertical track irregularity by utilizing ABAQUS explicit dynamics simulation module with a wheel-rail three-dimensional solid finite element dynamic model. The important influence of short wave track vertical irregularity on vertical vibration of wheel/rail system is discovered. The paper provides a theoretical basis for the study of high-frequency nonlinear vibration characteristics in wheel-rail system.

2. The Linear Inherent Characteristic Calculation of the Vertical Vibration in Wheel System

The structure of type wheel is shown in Figure 1; the wheel tread surface is shown in Figure 2.

In the vertical track irregularity excitation, due to the role of spring damper, high-frequency vibration component of the bogie frame is greatly weakened, which is mainly on the wheel set. Therefore, in the analysis of high-frequency vibration of wheel-rail system, model above incentive points ignores spring structure and only retains wheel set. The wheel and axle section are set up on sketch, and the section is rotated in order to obtain the wheel finite element model.

The ABAQUS finite element wheel model is established by using three-dimensional deformable elastomeric structure type. The global seeds dimension is 0.5 and the unit type is C3D20. The model is divided into 3552 units. The model can satisfy the requirement of calculation accuracy and speed. The finite element model is shown in Figure 3.

Lanczos algorithm is chosen to do the matrix operations through linear perturbation frequency method, and the high order natural frequencies can be effectively solved by using ABAQUS program. In this paper, the finite element model natural frequencies are solved. Natural frequencies calculation results are presented in Table 1.

Through the ABAQUS program, all the natural frequencies and vibration modes from 50 Hz to 5000 Hz are solved, totally 63 orders. Vertical, transverse, and bending vibrations are included. There are 16 vertical vibration natural frequencies, and the frequency distribution is relatively dispersed.

In order to verify the accuracy of the finite element model natural frequency, the conventional spring-mass method can be used to make a comparison. The spring-mass method is commonly used to solve the natural characteristic, which converts the wheel axle into Timoshenko beam with variable cross section and then the wheel is converted into concentrated mass.

The differential equations of the motion represented by the displacement of free vibration Timoshenko beam areIn the equation, is cross-sectional shape factor; is cross-sectional area; is the cross-sectional moment of inertia; is vibration displacement; is density; is the beam’s elastic modulus; is the beam’s shear modulus.

The general dynamic vibration differential equation isIn the equation, is mass matrix; is damping matrix; is stiffness matrix; is inertial element displacement vector; is inertial element velocity vector; is inertial element acceleration vector; is the matrix of system’s external load.

Spring-mass method is more suitable for solving the low band natural frequency, because the model complexity should not be too high. The natural frequency number solved by the spring-mass method is restricted to the model degrees of freedom. In condition of considering the model simplification without affecting the accuracy of solution, a six-degree-of-freedom wheel set model is established, which can solve a total of five natural frequencies.

, are defined. Then the differential equations of system for the free vibration without damping are obtained:According to the conventional theory of simple harmonic vibration method, the equation solution is set up asIn the equation, is natural circular frequency; is amplitude of inertia .

Substitute (4) into (3),In the equation, is the characteristic matrix; calculating natural frequency is converted into computing the algebraic Eigen values for . Meanwhile, the corresponding eigenvectors of the Eigen values are the vibration modes of natural frequencies; the five vertical vibration frequencies are 82.5 Hz, 166.8 Hz, 344.2 Hz, 627.5 Hz, and 820.1 Hz.

Compared with the result of conventional spring-mass method, the result difference of ABAQUS FE model is within 5%. Furthermore, the high order natural frequencies can be effectively solved by using ABAQUS program. It is presented that the ABAQUS FE model is reasonable.

3. Nonlinear Contact Dynamic Finite Element Model of Wheel-Rail System

Usually, in multibody dynamics and conventional spring-mass method, the system’s double nonlinear vertical vibration equation containing nonlinear damping and nonlinear stiffness is established [18].

For example, a two-degree-of-freedom system, Duffing oscillator, is used to define nonlinear stiffness, and then cubic nonlinear elastic force is obtained. Vanderpol oscillator is used to define nonlinear damping, and then cubic nonlinear friction force is obtained:In the equation, is nonlinear stiffness value; is nonlinear damping coefficient.

However, multibody dynamics method is mainly suitable for solving the low-frequency dynamics of integral structure problems such as the whole vehicle model on high-speed railway. The Runge-Kutta algorithm commonly used for solving ordinary differential equations in the multibody dynamics method is time consuming, and it is prone to divergence instability when solving the problem of high-frequency vibration.

The ABAQUS finite element analysis is quite valid to solve contact problems, such as high-frequency vibration and system’s nonlinear characteristics [19]. Based on the former analysis of system’s linear natural characteristics, the ABAQUS deformable contact nonlinear finite element dynamic model is used.

The model is the symmetrical complete one, which takes actual complexity of boundary condition at the symmetry axis into account without distorted rather than one-sided wheel model. In addition, the boundary constraint is relatively simple without taking other regulation of couple moment into consideration. The incorrect result caused by overconstraint is avoided.

When defining load and boundary condition, the load of gravitational field is taking vertical downward applied to the entire model. Commonly, the load of high-speed EMU wheel axle is from 14 t to 16 t. Spring damper is set to the wheel and track system; furthermore, spring damper of track system is connected to the ground. Then, the wheel velocity parameter is applied, as in Figure 4.

The type and quality of grid unit closely influence the accuracy, the smoothness, and spectrum of simulation results. In order to acquire the accurate solution of contact dynamics, the fine linear hexahedral hourglass-stiffness-control-type reducing-integral element C3D8R and hexahedral nonconforming element C3D8I are divided into wheel and track, respectively, in this paper [11]. The grid of wheel is separated quite densely, divided into 113,500 units. Grid seed of the track is laid out relatively sparsely but fine in the region of irregularity ensures the accuracy of wave curve, totally divided into 1860 units. The grid of model is shown in Figure 5. In this model, the same irregularity curves are set up in two tracks.

4. Nonlinear Dynamic Finite Element Model Algorithm and the Optimization of Contact Stiffness Factor Based on Experiment

4.1. Solving Nonlinear Dynamics of Wheel-Rail System Based on ABAQUS/Explicit by Improved Central Difference Method

For nonlinear dynamics, explicit dynamics analysis and implicit dynamics analysis are usually used. For implicit and explicit integral procedures,where is the node mass matrix; is the applied force; is the internal force; is the node acceleration.

Implicit dynamics analysis method is universally used to solve the general nonlinear problems. The method has no inherent limit on the size of time step increment. The solving time is typically dependent on the solution accuracy and convergence [20].

This method itself has the characteristics of unconditional stability. However, the stiffness matrix is very complicated. Particularly when the mesh is meticulous, it may lead to huge cost of iteration, the accuracy has no significant increase, and it has no obvious advantages in solving the high-frequency characteristics.

Explicit analysis method is quite effective on the nonlinear model with an incentive impact in short time, which has the advantage in the fast speed for each incremental step without iterative calculation. The fidelity and accuracy could be higher than the implicit analysis with any time integral parameter, which can solve the nonlinear high-frequency vibration system quite accurately. The explicit module of ABAQUS uses central difference method to give explicit solution of the motion’s equation with parametric improvement. The dynamic condition of incremental step is used to calculate the next one [21].

The dynamic explicit method is used in this paper. Explicit method requires minor step increment, which only depends on the model highest natural frequency and is unrelated to the type and duration of load. At the beginning of increment (), the acceleration isIt is simple to calculate the acceleration, since explicit method always uses diagonal or concentrate mass matrix.

Constant acceleration for the change in calculating velocity is assumed. To determine the midpoint velocity of incremental step , changes in the velocity value are added to the former midpoint velocity :The beginning increment displacement plus the integral of velocity to determine the displacement at the end of increment is as follows: The acceleration with the condition of dynamics balance is provided at the beginning of incremental step. When the acceleration is obtained, the velocity and displacement are calculated “explicitly” on time [19].

In order to ensure the analysis correctness for the elastic body of large displacement, the geometric nonlinear must be considered in the nonlinear contact.

4.2. Contact Parameter Optimization

Because the contact dynamic calculation is required in nonlinear analysis of the finite element model, contact parameter has great influence on the simulation results of the amplitude and frequency spectrum. Therefore, contact theory is the key to nonlinear dynamics. Surface to surface contact is chosen to be the interaction mode. Track surface is set as the first surface while wheel surface is set as the second surface. The mesh of first surface is divided more sparsely than the second one. A rational weighting coefficient of contact could prevent wheel nodes from inserting into the track surface, which leads a more accurate result in the important region. Due to the distinct time-varying characteristic of contact surface, it is necessary to define the slip formula as the finite slip.

Then, the contact property should be defined accurately and scientifically. The contact stiffness is quite important to dynamic calculation. The larger the contact stiffness is, the smaller the penetration is. However, excessive contact stiffness may lead the overall stiffness matrix into morbid state; the convergence difficulty of solution is aroused [22]. Based on Hertz contact theory, the contact stiffness factor is usually set from 0.1 to 1.

The optimal value not only makes the results more reasonable but also makes the displacement and acceleration achieve the best convergence effect. The specific value should be obtained through experiment and simulation of several attempts.

The peak accelerations and the convergence of various operating conditions are obtained by different stiffness factor values. It is shown in Figure 6 and Table 2.

In principle, the appropriate contact stiffness factor is easy to converge. Factors 0.1, 0.3, 0.4, and 0.7 are appropriate. For further verification of the stiffness factor, it is necessary to select the best one according to the actual experimental results.

4.3. Vibration Measurement Experiment of Wheel

When the load applied in different track irregularities and different high velocities, it is difficult to measure the wheel axle vertical acceleration, especially difficult to acquire the high-frequency characteristics within a very short period. Because of the elastic contact in wheel-rail system, the actual wheel sets on irregularity excitation is not a simple harmonic excitation. However, in general, only the equivalent technique is used to verify the amplitude. At present, the vertical displacement excitation of harmonic wave is usually equivalent to the track irregularity in relevant test.

According to the vibration measurement experiment of wheel by China Academy of Railway Sciences, different types of harmonic waves are taken as the track irregularity. The corresponding relationship iswhere , is the frequency of track irregularity; is waveform; is the amplitude of track irregularity; is the velocity of wheel; is the phase of track irregularity, two tracks with the same phase; is time; is the wavelength of track irregularity.

The axle vertical accelerations are tested in the condition of 16 t axial loaded, and then the vertical peak accelerations are presented in Table 3.

Combined with the test results and simulation results under different stiffness factors mentioned above, it can be concluded that the average difference of factor 0.1 is 5.89%, the average difference of factor 0.3 is 8.71%, the average difference of factor 0.5 is 8.25%, and the average difference of factor 0.7 is 2.73%.

It is presented that the contact stiffness factor 0.7 is the optimal value. The results analyzed are close to the measured value of these experiments with a fine convergence.

5. The Vertical Nonlinear Vibration Simulation Results and Analysis

Predecessors have done some related research including wheel-rail system joint spectrum such as Sato formula to input stimulated spectrum of wheel-rail random high-frequency vibration and noise radiation model. The equation is [23]where is spectral density of wheel-rail surface roughness; is roughness wave number; is wheel-rail surface roughness coefficient.

Wheel-rail under load action has elastic deformation, and the formation of an elliptical contact zone occurs. The influence of wheel-rail surface roughness on the wheel-rail system’s vibration is affected by the size of contact area. When the roughness wavelength is less than the contact patch, the excitation effect will be weakened. Considering the effect, filter function is used:where is the roughness correlation coefficient; is the contact zone circle radius; is Bessel function; is the longitudinal distance along the track.

In this paper, the authors set the wavelengths from 0.2 m to 2.0 m. The irregularity amplitudes are generally from 0.2 mm to 1 mm based on the actual tests; 0.3 mm, 0.6 mm, and 0.9 mm are chosen as examples for simulation. Wheel axle is the focus. According to the velocity criteria of China high-speed railway, 200 km/h and 350 km/h are selected to be main velocities.

Figures 714 present the axle vertical acceleration in some typical conditions and the results of Fourier spectrum. The spectrum results of low frequency below 400 Hz are ignored.

5.1. The effect of 2 m Wavelength-0.9 mm Amplitude Irregularity

When the velocity is 200 km/h, the peak acceleration is 60 m/s2; acceleration amplitude is 30 m/s2. It is shown on Fourier spectrum that the main dominant frequency is 1403 Hz while the secondary dominant frequencies are 2382 Hz and 770.3 Hz. Because of the wave filter and restriction on the sampling number, the highest frequency is generally lower than 5000 Hz.

When the train velocity is 350 km/h, the peak acceleration is 150 m/s2; acceleration amplitude is 97 m/s2. It is shown on Fourier spectrum that the main dominant frequency is 1938 Hz, and the secondary dominant frequencies are 2372 Hz and 1403 Hz.

5.2. The Effect of 1 m Wavelength-0.9 mm Amplitude Irregularity

When the velocity is 200 km/h, the peak acceleration is 130 m/s2; acceleration amplitude is 72 m/s2. The main dominant frequency is 1403 Hz; the secondary dominant frequency is 2306 Hz.

When the velocity is 350 km/h, the peak acceleration is 720 m/s2; acceleration amplitude is 340 m/s2. It is shown on Fourier spectrum that the main dominant frequency is 1914 Hz, and the secondary dominant frequencies are 2324 Hz and 3600 Hz.

5.3. The Effect of 0.4 m Wavelength-0.9 mm Amplitude Irregularity

When the train velocity is 200 km/h, the peak acceleration is 400 m/s2; acceleration amplitude is 160 m/s2. It is shown on Fourier spectrum that the main dominant frequency is 2372 Hz, and the secondary dominant frequencies are 1403 Hz and 4344 Hz.

When the train velocity is 350 km/h, the peak acceleration is 950 m/s2; acceleration amplitude is 490 m/s2. It is shown on Fourier spectrum that the main dominant frequencies are 1914 Hz and 2324 Hz, and the secondary dominant frequency is 4283 Hz.

5.4. The Effect of 1 m Wavelength-0.6 mm Amplitude Irregularity

When the train velocity is 200 km/h, the peak acceleration is 91 m/s2; acceleration amplitude is 49 m/s2. It is shown on Fourier spectrum that the main dominant frequency is 2372 Hz, and the secondary dominant frequency is 1403 Hz.

When the train velocity is 350 km/h, the peak acceleration is 310 m/s2; acceleration amplitude is 255 m/s2. It is shown on Fourier spectrum that the main dominant frequencies are 2324 Hz and 1914 Hz, and the secondary dominant frequencies are 3554 Hz and 4283 Hz.

As can be seen from the results list above in different working conditions, changes of track irregularity wavelength will affect Fourier spectrum of the corresponding acceleration curve. The main dominant frequency varies significantly when the velocity is 350 km/h.

Finally, in order to make the analysis comprehensive and persuasive, the analysis includes the conditions of eight irregularity-wavelengths from 0.2 m to 2 m. Combine with the simulation above; Figure 15 and Table 4 are obtained.

It is concluded through the analysis of the simulation results in various working conditions that(1)the main dominant frequency includes 1400 Hz, 1900 Hz, and 2350 Hz;(2)in the majority of operating conditions, the peak acceleration at velocity of 350 km/h is higher than that at 200 km/h significantly, and the main frequency also shifts or changes. Meanwhile, the nonlinear strength of system is enhanced obviously.In the past few years, many experts have made a conclusion about the influence of track irregularity on wheel set. When the speed is constant, the increase of irregularity wavelength will lead to the peak acceleration and nonlinearity decreasing; the increase of irregularity-amplitude will lead to the peak acceleration increasing. However, the decrease of wavelength or amplitude has an opposite result. When the amplitude and wavelength are constant, the decrease of velocity will weaken the vibration.

However, it is not the case. In vibration, it does not have a single dynamic characteristic for a system. Gao et al. had discovered the sensitive wavelength of the long-wavelength track irregularity on the vehicle vibration [24].

In this paper, it is presented from the simulation results that vertical track irregularity (wavelength less than 2.0 m) has a critical wavelength on the wheel vertical vibration because of its inherent characteristics and nonlinear factors. When irregularity wavelength is close to the critical range, the dynamic characteristics of vibration change obviously are summarized as follows.

For the irregularity whose wavelength is more than 0.8 m, whether the velocity is 200 km/h or 350 km/h, the acceleration is reduced by decreasing the irregularity-amplitude or increasing wavelength. The increase of irregularity wavelength will significantly lead to the main frequency amplitude decreasing when the velocity is 350 km/h. The main frequency amplitude changes little when the velocity is 200 km/h.

However, for the irregularity whose wavelength is less than 0.8 m, one has the following.

The increase of irregularity wavelength will significantly lead to the main frequency amplitude decreasing when the velocity is 200 km/h. However, the main frequency amplitude changes little when the velocity is 350 km/h.

The peak acceleration is not sensitive to the vertical track irregularity-amplitude when the velocity is 200 km/h. However, it is sensitive to the wavelength. The decrease of wavelength will significantly increase the acceleration, while the increase of irregularity-amplitude affects little increasing the acceleration.

The peak acceleration is not sensitive to the vertical track irregularity wavelength when the velocity is 350 km/h. However, it is sensitive to the amplitude. The increase of amplitude will significantly increase the acceleration, while the decrease of irregularity wavelength affects little increasing the acceleration.

It is shown in Table 4 that, for example, under the action of 0.2 m wavelength-0.3 mm amplitude irregularity, the vertical peak acceleration at 200 km/h is obviously 2.5 times the value at 350 km/h. Under the action of 0.6 m wavelength-0.3 mm amplitude irregularity, the vertical peak acceleration at 200 km/h is almost equal to the value at 350 km/h. However, under the action of 0.6 m wavelength-0.9 mm amplitude irregularity, the peak acceleration at 350 km/h is 3.3 times the value at 200 km/h. Under the action of 0.8 m wavelength-0.9 mm amplitude irregularity, the peak acceleration at 350 km/h is 6 times the value at 200 km/h. It is explained that the acceleration at 300 km/h is not always absolutely higher than the value at 200 km/h. Sometimes, under the same conditions of vertical irregularities, the vibration is less when the speed is high.

It can be determined that the critical wavelength is about 0.8 m.

In addition, special attention should be paid to some conditions that the main dominant high-frequency amplitude is more than 50% of the acceleration amplitude. From the viewpoint of energy, the abnormal state of spectrum amplitude is definitely not being overlooked. Even if the peak acceleration is low, the destruction will be huge to the internal of wheel-rail system. Furthermore, these working conditions are mostly the “soft” vibration conditions at 200 km/h. It is presented that the wheel vibration response under the action of a vertical track irregularity is not always positively related to the train speed.

6. Conclusions

In this paper, Based on the theory of nonlinear contact dynamics of wheel-rail system, a symmetrical complete three-dimensional finite element model of wheel-rail system is established. The authors discover the important influence of vertical track irregularity on the wheel set high-frequency vibration.(1)The inherent linear characteristic of wheel set is solved. Compared with the result of conventional spring-mass method, the result difference of ABAQUS FE model is within 5%. Furthermore, the high order natural frequencies can be effectively solved by using ABAQUS program. It is presented that the ABAQUS FE model is reasonable.(2)The improved central difference method by ABAQUS explicit nonlinear dynamics module is suitable for calculating the high-frequency vertical vibration of wheel axle. Through the vibration measurement experiment of wheel by China Academy of Railway Sciences, the contact stiffness factor is optimized, defined as 0.7.(3)The nonlinear vibrations and frequencies response characteristics of wheel axle are simulated. It is shown that the main frequencies of vertical vibration acceleration on wheel axle in different conditions are concentrated in the range of 1400 Hz to 4500 Hz, very close to the vertical vibration natural frequency of 7th, 8th, 9th, 10th, 13th, and 15th mode with the obvious characteristics of high-frequency vibration. In addition, the frequencies offset in a certain range is affected by nonlinear characteristics.(4)The critical wavelength of vertical track irregularity on the vertical vibration at wheel axle of wheel-rail system is about 0.8 m. When the velocity is 350 km/h, system’s nonlinear characteristics are higher than the results at 200 km/h. Furthermore, a particular attention should be paid to the Fourier spectrum whose “energy” is quite strong even if the peak acceleration is quite “weak,” especially whose main dominant frequency amplitude is more than 50% of the acceleration amplitude.

The authors believe that the conclusion of this paper is of relatively high value. Referring to the amplitude-frequency dynamics characteristics, the appropriate measures are formulated, which have special guiding significance to the on-site tests and have already been applied to the operational maintenance.

Additional Points

Research area includes dynamics and infrastructure inspection.

Disclosure

Yang Lei is the first author.

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

First, it is an honor for the authors that the paper is based on the project named “Investigation of Long-Term Vibration Response and Deterioration Mechanism of a HS Railway Tunnel Subjected to Train Loading and Aerodynamic Loading under the Cumulative Effect of the Typical Initial Damages.” The “Joint Fund for Basic Research of High-speed Railway” whose Award no. is U1434211 supports the project. “The Joint Fund for Basic Research of High-speed Railway” focuses on the basic theoretical research of the high-speed train and railway engineering technology, which is designed to play the guiding and coordinating role of the National Natural Science Foundation of China (NSFC). The recipient is researcher Qianli Zhang. The authors would like to express the heartfelt thanks to them. This works is supported by Joint Fund for Basic Research of High-speed Railway of National Natural Science Foundation of China, China Railway Corporation (U1434211).