# Stochastic Analysis of Nonlinear Vehicle-Track Coupled Dynamic System and Its Application in Vehicle Operation Safety Evaluation

Authors: Fu-Shan Liu; Zhi-Ping Zeng; Wei-Dong Wang; Abdulmumin Ahmed Shuaibu
Journal: Shock and Vibration
 (2019)

---

## Abstract

This paper presents a vehicle operation safety evaluation model; to this end, a nonlinear vehicle-track coupled dynamic system stochastic analysis model under random irregularity excitations based on probability density evolution method was developed. The nonlinear coupled vehicle-track dynamic system is used to accurately describe the wheel-rail contact state. The stochastic function-spectral representation is used to simulate the random track irregularity in the time domain for the first time; consequently, the frequency components in the irregularity are preserved and random variables are reduced. In the process of evaluating the safety of train operation, the probability evolution, reliability of evaluation indices for different limit values, and evaluation indices for different probability limits are calculated for more accurate evaluation. The dynamic model and safety evaluation method was verified using the Zhai-model and Monte Carlo method. The results show that, when the probability guarantee is increased, the running safety index of the vehicle increases more rapidly with running speed and the left/right wheel-rail derailment coefficient increases rapidly at running speeds above 400 km/h. The computational model provides a novel direction for vehicle operation safety evaluation.


---

## Body

## 1. Introduction

Vehicle-track coupled dynamic system is a nonlinear and stochastic system in nature as random track irregularity is one of its main random excitations. Meanwhile, the vehicle operation safety under random excitations is also the main transportation engineering concerns. In the past few decades, the popularization of high-speed railway and the continuous demand for high operating speed have brought many engineering challenges. In many countries, evaluation standards for train operation safety have been emplaced. Due to the complexity involved in the safety assessment of train operation, the problems generally need to be addressed through simulation methods. Several different train-track coupled dynamics models or train-bridge coupled dynamics models have been established to analyze the dynamic response of the system. Wu et al. [1] developed a vehicle-rail-bridge interaction model for analyzing the 3D dynamic interaction between the moving trains and railway bridge; in this established model, the wheelset has 3 degrees of freedom, the vertical, lateral, and rolling, and no jumps occur between the wheelsets and rails. After that, Chen and Zhai [2] established two dynamic models, train running on the bridge with ballasted track and the ballastless slab track. In these two models, the jump between the wheelset and the rail can be simulated. In reference [3], Zhai further extended these models to three-dimensional space, and at the same time, a three-dimensional space wheel-rail contact relationship has also been developed. Lou and Zeng [4] have derived the equations of motion in matrix form for the railway freight vehicle-bridge interaction system based on energy approach; this modeling method also has been widely used. Dinh et al. [5] have developed a formulation of three-dimensional dynamic interactions between a bridge and a high-speed train using wheel-rail interfaces; in this model, contact loss is allowed, and the contact is nonlinearly. Other than these, many scholars, such as Zhang et al. [6] and Yang and Yau [7], also have established their own vehicle-track or vehicle-track-bridge coupling dynamics models. The development of vehicle-track coupled dynamics modeling facilitates the evaluation of vehicle operational safety through numerical simulation. However, as the main excitation source in the train operation, the track random irregularity cannot be measured and predicted before the completion of the line construction, it also changes with time and environment. Therefore, it is necessary to consider the train safety operation from a stochastic point of view. Many scholars have tried to use random vibration analysis method to study the dynamic response of vehicle-track coupled dynamic system and have obtained many meaningful research results.Due to the complexity of vehicle-track coupled dynamic system, Monte Carlo method (MCM) is a common analytical method in random analysis. Chen et al. [8] and Xia et al. [9] studied the random vibrations in a nonlinear vehicle-track coupling dynamic system by the MCM and established an early nonlinear analysis model of the system from a stochastic perspective. Subsequently, many efficient random vibration methods were applied to vehicle-track-coupled power systems. Lu et al. [11], Zhang et al. [12], and Zeng et al. [13] established random vibration analysis model for linear vehicle-track-bridge coupling system by pseudoexcitation method (PEM). This method analyzes the random vibrations in linear systems efficiently, especially the power characteristics in the frequency domain. Zhang et al. also have used the PEM and the precise integration method (PIM) to compute the nonstationary random responses of three-dimensional train-bridge systems subjected to lateral horizontal earthquakes [14]. Yu et al. [15] studied the random vibrations of a linearized vehicle-rail-bridge coupling system by generalized probability density evolution theory (originally formulated by Li and Chen [16]), which is more efficient than the MCM. Xu et al. [17, 18] proposed a stochastic space-time analysis method for vehicle-track-coupled power systems. Their method also applies generalized probability density evolution theory and is based on the measured track irregularity data. Their probabilistic model efficiently and accurately analyzes the vehicle-track-coupled power system under random track irregularities.The introduction of stochastic vibration analysis method makes it possible to use the probability theory to evaluate the safety of vehicle-track coupled power system or vehicle-bridge coupled dynamic system. Rocha et al. presented a probabilistic methodology for the safety assessment of short-span railway bridges for high-speed traffic [19] and analysed the sensitivity of the dynamic response due to the variability of the main structural parameters [20], and the efficiency of different probabilistic methodologies for safety assessment of short-span railway bridges is compared [21]. Using the principle of virtual work, utilizing a linearized wheel-rail contact equation to avoid iterative solution at each time step of the integration, and applying rail irregularity as random excitations to the system, the probabilistic evaluation approach for nonlinear vehicle-bridge dynamic performances has been analysed by Jin et al. [22]. By developing an improved response surface method for nonlinear limit states, Cho et al. [23] evaluated the reliability of vehicle and high-speed railway bridge system.In the existing random vibration analysis of vehicle-track coupled dynamics, due to the limitation of stochastic analysis methods or in order to improve the computational efficiency, linear or linearized models are often used. For the linearization of a model, the wheel-rail relationship is often reduced to a linear model. The wheel-rail relationship not only reflects the force state between wheel and rail, but also the excitation of track irregularity also needs to be input into the system through wheel-rail relationship; this makes the wheel-rail relationship one of the core issues in the modeling. The linearized model will not be able to accurately simulate the change of the size, direction, and action point of the wheel-rail force, and it will also lead to a difference of the input track irregularity excitation input. Other than this, models with nonlinear wheel-rail contact relationship allow for the evaluation of vehicle safety in a more realistic way than linear methods [24]. Numerous studies on vehicle-bridge dynamics have indicated that nonlinear interaction has profound effects on the system responses [25]. Therefore, the nonlinear vehicle-track coupled dynamics model is more helpful for vehicle-safety evaluation, and the calculation results are more accurate.In this paper, a stochastic analysis model of nonlinear vehicle-track coupled dynamic system under random irregularity excitations was established based on probability density evolution method. The nonlinear wheel-rail relationship has been considered to facilitate more convenient and accurate evaluation of vehicle operational safety. Generalized probability density evolution theory has been used, which has been widely applied to linear and nonlinear stochastic systems and have proven to obviously improve the efficiency the analysis [15, 16, 26]. In addition, function-spectral representation method has been used to accurately simulate the random track irregularities in the time domain, which retains a large number of frequency components in the samples while reducing the required number of random variables. It also avoids the clustering phenomenon in the point selection of high-dimensional random variables. The probability density evolution of the left/right wheel-rail load reduction, the left/right wheel-rail derailment coefficient, and the vertical/transverse acceleration of the car body were calculated by the generalized probability density evolution theory. Based on the results, a safety evaluation model of train operation is presented and validated by comparison with previous study. The developed computational model provides a new research direction for train safety evaluation.


## 2. Stochastic Analysis Model of Nonlinear Vehicle-Track Coupled Dynamic System

Figure1 shows the calculation flow of the nonlinear vehicle-track coupled dynamic system stochastic analysis model, which is based on generalized probability density evolution theory. The main steps are the establishment of the nonlinear system, the numerical simulation of Npt random track irregularities time-domain samples, the dynamic response calculation of the coupled vehicle-track system under random excitations, calculation of the probability density evolution of the safety evaluation index, and the train safety evaluation.Figure 1
Flow chart of the method for evaluating train operation safety.


### 2.1. Stochastic Analysis Method Based on Generalized Probability Density Evolution Theory

#### 2.1.1. Generalized Probability Density Evolution Theory

Under random irregular excitations, the nonlinear vehicle-track coupled dynamic system is expressed as(1)Mξ¨t+Cξ˙t+Kξt=FΘ,t,where t denotes time, and ξt, ξ˙t, ξ¨t are the displacement, velocity, and acceleration vectors of the coupled system, respectively. The matrices [M], [C], and [K] are the mass, damping, and stiffness matrices, respectively, reflecting the mechanical and geometric nonlinearity in the wheel-rail relationship. F(Θ, t) is the input excitation vector of track random irregularities, and Θ is the set of sampled random variables.Most of the randomness in a vehicle-track coupled dynamic system derives from the random track irregularities. A nonlinear vehicle-track dynamic systemξ (t) is a probabilistic conservative stochastic system that excludes other random sources. According to the probabilistic conservation principle of stochastic systems, probabilistic conservation-augmented systems are constructed as [ξ (t), Θ] with(2)Pξt,Θ∈Ωt×Ωθ=Pξt0,Θ∈Ωt0×Ωθ,where P{·} denotes the probability of random events. Ωt and Ωt0 are the time distribution regions of the system state at times t and t0, respectively, and Ωθ is the space distribution region of the system state of sample θ. The above expression is equivalent to(3)ddtPξt,Θ∈Ωt×Ωθ=0.Expressing the joint probability density function (ξ(t), Θ) as pξθξ,θ,t, the above expression becomes(4)DDt∫Ωt×ΩθpξΘξ,θ,tdξdθ=0,where D/Dt is the derivative of matter [16].The generalized probability density evolves by the following equation [16]:(5)∂pξΘξ,θ,t∂t+ξ˙θ,t∂pξΘξ,θ,t∂ξ=0.In this formula,ξ˙θ,t is the free velocity of ξ (t) under the random sample {Θ = θ}. Its initial probability is PξΘξ,θ,tt=0=δξ−ξ0PΘθ.Solving the above evolution equation of generalized probability density under the initial conditions, we obtain the joint probability density function of sampleθ, pξΘξ,θ,t. Furthermore, the probability density function of PξΘ is obtained by integrating θ:(6)PξΘξ,t=∫ΩΘpξΘξ,θ,tdθ.For the probabilistic conservative vehicle-track coupled dynamic system under the random track irregularity excitation, the system dynamic responseξ satisfies the generalized probability density evolution equation shown in equation (5). Simultaneously, the generalized probability density evolution equation decouples the probability space from the physical space. As shown in Figure 1, first solve the vehicle-track coupled dynamic system equation (1) to obtain ξ˙θ,t under each random sample θ, and then solve the generalized probability density evolution equation for each random sample θ. When all random samples have been calculated, by solving equation (6), the probability density evolution process of the dynamic response can be obtained. Finally, a stochastic analysis of the nonlinear vehicle-track coupled dynamic system is realized.


#### 2.1.2. Selection of Random Variable Samples and the Probability Assignment

In generalized probability density evolution theory, selecting the random variable samples is directly related to the computational efficiency of the random analysis. In this paper, the random variables are described by a low-deviation sequence selected by number theory, which reduces the required number of random samples and accelerates the convergence of the calculation. However, as the low-deviation sequence is uniformly distributed, it cannot be selected from a high-dimensional space without forming void and cluster structures. Therefore, the uniformity will degrade and the calculation results will contain large errors.Random irregularity is the superposition of irregular waves of different wavelengths, different phases, and different amplitudes. To accurately simulate the random track irregularities, more random frequencies and phase components (i.e., more selections) are needed. To avoid clustering in the random sample selection, the random track irregularities were simulated in the time domain by a stochastic function-spectral representation method, which requires just eight random variables. The generation method and validation of the sampling will be discussed in Section2.4. The random variables were sampled as Sobol sequences.The random samples in the probabilistic space ΩΘ must then be divided. Probabilistic spaces are commonly divided into representative Voronoi regions with representative volumes V (θk) satisfying(7)⋃kNptVθk=ΩΘ,Vθi∩Vθj,i≠j.The assigned probabilityPk of random sample θk is given by(8)Pk=∫VθkpΘθkdv,where pΘθk is the probability density of random variable Θ in Voronoi region V(θk). The joint probability density in the total partitioned space satisfies(9)∑k=1NptPk=∑k=1Npt∫VθkpΘθkdv=1.This method obtains the initial probability of each random point sample.


### 2.2. Vehicle-Track Coupled Dynamics Model

Figure2 shows the vehicle-track coupling dynamic model. The vehicle is modeled by multirigid body dynamics, and the rail is regarded as a finite-length, elastic-point supported Euler beam. The dynamic vibration equation of the beam is solved by the Ritz method. The wheel-rail relationship is modeled by nonlinear Hertz contact theory and wheel-rail nonlinear creep theory. The fastener is modeled by a linear spring damping element in the transverse and vertical directions. The track structure is modeled by the finite element method; the slab, mortar layer, and base are modeled as plate elements. The foundation support is simplified as a linear spring damping surface element. The wheel-rail relationship connects the vehicle and rail models, whereas the fastener force connects the rail model and the track slab-foundation structure.Figure 2
Vehicle-track coupled dynamics model.


#### 2.2.1. Vehicle Dynamics Model

Figure3 schematizes the vehicle model composed of one car body, two bogies, and four wheelsets. Each part is treated as a rigid body with freedom in the longitudinal, transverse, bounce, roll, and yaw directions (giving 35 degrees of freedom per vehicle).Figure 3
Vehicle dynamics model.The antisnake shock absorber significantly affects the driving quality when the train is running at a high speed. Treating the force of the antisnake shock absorber on the vehicle and the antirolling moment of the vehicle’s central suspension as an additional force on the vehicle, the dynamic equation of the vehicle is given by(10)MVξ¨V+CVξ˙V+KVξV+GV=FV,where [MV], [CV], [KV], [GV], [FV] are the vehicle mass, damping, stiffness, additional force, and load matrices, respectively.The freedom matrix of the vehicle is expressed as(11)ξV=ξc,ξb1,ξb2,ξw1,ξw2,ξw3,ξw4T,where the subscripts c, b, and w represent the car body, bogie, and wheelset, respectively, and the subscripts i and j represent the ith bogie and the jth wheelset, respectively. The degrees of freedom of each part are expressed as(12)ξc,bi,wj=Yc,bi,wj,Zc,bi,wj,ψc,bi,wj,βc,bi,wj,φc,bi,wj,where Y, Z, Ψ, φ, and β represent the freedoms in the longitudinal, transverse, bounce, roll, and yaw, respectively.The vehicle-mass matrix is given by(13)MV=Mw1,Mw2,Mw3,Mw4,Mt1,Mt2,Mc,where Mwi,bj,c=diagmwi,bj,c,mwi,bj,c,Ixwi,bj,c,Iywi,bj,c,Izwi,bj,c. Here m denotes the mass, and Ix, Iy, and Iz are the moments of inertia around the X, Y, and Z axes, respectively.The vehicle-stiffness matrix is given by(14)KV=Kw1000Kwb100Kw200Kwb200Kw300Kwb30Kw40Kwb40Kb10Kbc1Kb2Kbc2symsKc,where(15)Kwi=2kpykpzkpzdw20kpxdw2,i=1,2,3,4,Kbj=22kpy+ksy0−2kpyHtw+ksyHbt002kpz+ksz0002kpyHtw2+2kpzdwk2+kszdsk2+ksyHbt2002kpzlt2+2kpxHtw2+ksxHbt20syms2kpylt2+2kpxdwk2+ksxdsk2,Kc=4ksy0−ksyHcb00ksz000ksyHcb2+kszdsk200kszlc2+ksxHcb20symsksylc2+ksxdsk2,Kwbi=2−kpy0kpyHtw0−1ikpylt−kpz0−1i+1kpzlt0−kpzdwk20000syms−kpxdwk2,i=1,2,3,4,Kbcj=2−ksy0ksyHcb0−1j+1ksylc−ksz0−1j+1kszlc0−kszdsk2+ksyHcbHbt0−1jksylcHbtksxHcbHbt0syms−ksxdsk2,j=1,2.The vehicle-damping matrix(16)CV=Cw1000Cwb100Cw200Cwb200Cw300Cwb30Cw40Cwb40Cb10Cbc1Cb2Cbc2symsCc,is similar to the vehicle-stiffness matrix. In fact, the vehicle-damping matrix is obtained by replacing the stiffness values in the stiffness matrix with the corresponding damping values.The additional moment on the vehicle is expressed as the matrix(17)G=Gw1,Gw2,Gw3,Gw4,Gt1,Gt2,Gc,with(18)Gwi=0,0,Iywiβ˙wi−Ωφ˙wi,0,Iywiβ˙wi−Ωψ˙wiT,Gtj=0,0,0,FxsLj+FxsRjHbt,FxsLj−FxsRjdsxT,Gc=0,0,Mr1+Mr2,FxsL1+FxsR1+FxsL2+FxsR2Hcb,FxsL1+FxsL2−FxsR1−FxsR2dsx.In the above expressions,Gwi is the matrix of moments acting on the ith wheelset caused by wheelset rolling rotations, Gtj is the matrix of moments acting on the jth bioge caused by absorber, Gc is the matrix of moments acting on the car body caused by absorber, and FxsLj and FxsLj are the additional forces imposed by the snake-shaped shock absorber. They are calculated as follows:(19)FxsL,Rj=FmaxvxctL,Rjv0,vxctL,Rj<v0,FmaxsignvxctL,Rj,vxctL,Rj≥v0.M
ri is the antiroll moment, given by(20)Mri=Krxψc−ψti,i=1,2.The vehicle-load matrix is expressed as(21)F=A1Fwr1,A2Fwr2,A3Fwr3,A4Fwr4,05×1,05×1,05×1,where Ai is the matrix of conversion coefficients between the wheel-rail force and the vehicle freedoms. Fwri is the wheel-rail force array, with(22)Awi=01010000000000−100−10000000−rLia00−rLi−a0000000rLirLiφwi0rRirRiφwi0010010a0a0φwi0−a0−a0φwi0001001,Fwri=FLxi,FLyi,FLzi,FRxi,FRyi,FRzi,McLxi,McLyi,McLzi,McRxi,McRyi,McRzi,where FL,Rxi = FnL,Rxi + FcL,Rxi, FL,Ryi = FnL,Ryi + FcL,Ryi, and FL,Rzi = FnL,Rzi + FcL,Rzi. Fn and Fc are the normal and creep forces, respectively, between the wheel and the rail. The subscripts L and R represent the left and right sides, respectively, and the subscripts x, y, and z denote the upper component in the X-, Y-, and Z-axis directions, respectively.


#### 2.2.2. Dynamic Model of the Track Structure

The track structure is divided into two main parts: the rail and the track plate. The rail is treated as an elastic support for an Euler beam, and the track plate and base plate are modeled by a plate element with high rigidity in the transverse direction, lateral displacement using beam element to simulate. Since the numerical integration solution step length is short and requires a lot of numerical calculation work, the method of the literature [3] is used in the modeling of the rail, and the dynamic response of the rail is obtained by solving the rail vibration differential equation, which can greatly reduce the calculation work. At the same time, the track plate is modeled by the finite element method to facilitate the modeling of various ballastless track structures. The track-structure model and rail forces are shown in Figures 4 and 5, respectively.Figure 4
Structural model of the rail-track plate.Figure 5
Rail force diagram.The following expressions are shown for the right rail (the left rail is treated analogously). The fastener force matrix is expressed as(23)FpviFphiMpi=0kpv00−kpvkph0kphHra−kph0−kphHra0kphHra2+kpvLr2kphHra0Yri,Zri,ψri,Ysi,ZsiT.The torque of the wheel-rail force, which induces rail torsion, is(24)Mcr=FRyHrb−FRzLce.The differential equations of the vertical, transverse, and torsional vibrations of the rail are, respectively given in the following [3].Vertical:(25)EIrz∂4Yrx,t∂x4+mr∂2Yrx,t∂t2=−∑i=1NFpyitδx−xi+∑j=14FRyjtδx−xwj.Transverse:(26)ErIrz∂4Zrx,t∂x4+mr∂2Zrx,t∂t2=−∑i=1NFpvitδx−xi+∑i=14FRzjtδx−xwj.Torsional:(27)ρrIr0∂2ψrx,t∂t4−GrIt∂2ψrx,t∂x2=−∑i=1NMpitδx−xi+∑j=14Mcrjtδx−xwj.The track and base plates are considered as horizontally rigid bodies and are modeled as four-node bending plate elements. They are connected by a fastener element, which is simulated by a linear spring damping element. Slab-base has been modeled by the principle of total potential energy with stationary value in elastic system dynamics [4], the dynamics of the slab and base plates are governed by(28)MssMbbδ¨sδ¨b+CssCsbCbsCbbδ˙sδ˙b+KssKsbKbsKbbδsδb=FsFb,where the subscript s and b denote the track and base plates, respectively. The vibration equation of the track and base plates assumes the principle of invariant potential energy of an elastic system. The foundation structure is modeled by the same method. [F] is the load matrix, and the main load acting on the slab-base is the fastening force. The matrix and modeling process of the track and base plates are detailed elsewhere [13].


#### 2.2.3. Wheel-Rail Contact Relation

The wheel-rail relationship links the vehicle system to the track system. In the vehicle-track coupled dynamics, the wheel-rail relationship mainly includes the wheel-rail normal-seeking solution, the wheel-rail creep solution, and the wheel-rail contact geometry. The wheel-rail relationship is shown in Figure6.Figure 6
Wheel-rail contact model.(1) The Wheel-Rail Normal Contact Model. The wheel tread and rail profile are complex curves composed of many simple curves with different radii. According to nonlinear Hertz contact theory (Figure 7), the shapes of the wheel and rail tread near the contact point are approximately described as(29)z1=12ρwx12+12rwy12,for the wheel tread,z2=−12ρrx22,for the rail tread,where rp, rw, and rr represent the curvature radii of the wheel tread, wheel rolling circle, and rail tread, respectively. Positive and negative curvatures indicate a convex and concave surface, respectively.The pressure distribution on the contact surface is given by Hertz contact theory.(30)p=Pn1−xa2−yb21/2.Figure 7
Simplified schematic of wheel-rail contact.The normal force on the wheel and rail is expressed as(31)FnL,R=43E∗Re1/2εnL,RF2e3/2.And the maximum pressure between the wheel and rail is given by(32)PnL,R=6FnL,RE∗2π3Re21/3F1e−2.In the above expressions,E∗ is the equivalent modulus of elasticity, given by 1/E∗=1−vr2/Er+1−vw2/Ew, and e=1−b2/a21/2, where a and b are the lengths of the long and short axes, respectively, of the contact ellipse (with b < a). F1(e) and F2(e) are correction factors that depend only on the ratio of relative curvature radii R′/R″. The variations of these factors with increasing R′/R″ are plotted in Figure 8.Figure 8
Correction factorsF1(e) and (F)2(e) versus ratio of relative curvature radius.F
1(e) and F2(e) can be directly computed by a formula or directly read from Figure 9. To improve the efficiency of the computation, the curve-fitting parameters obtained from Figure 9 are listed in Table 1, F1(e) and F2(e) were then obtained as(33)Fie=q1x5+q2x4+q3x3+q4x2+q5x1+q6x3+p1x2+q2x1+q3,with x = ln(B/A).Figure 9
Geometric wheel-rail contact model.Table 1
Fitting parameters of theF1(e) and F2(e) curves.
F
1(e)
F
2(e)
Parameters
Error of fitting
Parameters
Error of fitting
q
1
0.001227
SSE: 3.522e-07
q
1
0.002195
SSE: 3.156e-07
q
2
−0.05683
R-square: 1
q
2
−0.07143
R-square: 1
q
3
0.7315
Adjusted R-square: 1
q
3
0.4331
Adjusted R-square: 1
q
4
4.087
RMSE: 5.935e-07
q
4
4.816
RMSE: 5.618e-07
q
5
20.1
q
5
41.2
q
6
36.46
q
6
21.94
p
1
4.599
p
1
5.123
p
2
20.1
p
2
41.2
p
3
36.46
p
3
21.94When calculating the normal contact force between the wheel and rail, the relative compression between the wheel and a rail with random irregularities is obtained by geometrically solving the wheel-rail contact. The normal force between the wheel and rail is then obtained by solving the wheel-rail normal contact model and determining the direction and magnitude variations in the normal wheel-rail force at different contact points. Converting to the absolute coordinate system, the normal loads on the wheel and rail array are, respectively, expressed as(34)FnLi′=FnLxi,FnLyi,FnLziT=FnLisinδcLi+ψwisinφwi,−sinδcLi+ψwicosφwi,cosδcLi+ψwiT,FnRi′=FnRxi,FnRyi,FnRziT=FnRi−sinδcRi−ψwisinφwi,sinδcRi−ψwicosφwi,cosδcRi−ψwiT.Here,δcLi and δcRi are the left-sided and right-sided wheel-rail contact angles, respectively, obtained from the wheel-rail contact geometry (details will be given in Section (2)).(2) Wheel-Rail Tangential Contact Model. The wheel-rail tangential contact was solved by Kalker linear creep theory [27], corrected for nonlinearity. The creep force and creep torque between the wheel and rail after the nonlinear correction are given by(35)Fx′=α⋅Fx,Fy′=α⋅Fy,Mz′=α⋅Mz,where α is the correction coefficient, given as α=FR′/FR with FR=Fx2+Fy2 and(36)FR′=μFnL,RFRμFnL,R−13FRμFnL,R2+127FRμFnL,R3,FR≤3μFnL,R,μFnL,R,FR>3μFnL,R,Fx and Fy denote the creep forces in the X- and Y-axis directions, respectively, and Mz is the rotational creep moment. In linear creep theory, these components are, respectively, expressed as(37)Fx=−f11ξx,Fy=−f22ξy−f23ξsp,Mz=f23ξy−f23ξsp,where f11 and f22 are the longitudinal and transverse creep coefficients, respectively, and f23 and f33 are rotational creep coefficients:(38)f11=GwrabC11,f22=GwrabC22,f23=Gwrab3/2C23,f33=Gwrab2C33.In the above expressions,Gwr is the composite shear modulus of the wheel and rail materials, and Cij are the Kalker coefficients [27] (after Kalker, who calculated them). Meymandet al. [28] have fitted the Kalker coefficients under different contact conditions. Here, the Cijs are expressed as(39)C11=2.39+2.338vwr+1.014g+2.736vwr2−0.0728vwrg−0.0139g2,C22=2.353−0.0212vwr+1.074g+1.163vwrg−0.0319g2,C23=0.2866+0.4467vwr+1.049g+0.038vwrg−0.0177g2+0.005g3+0.1479g2vwr,C33=0.7375−1.177vwr+0.45061g+1.076vwr1g−0.0231g2,where g is the ratio of the long- and short-elliptical axes, g = a/b, and vwr is the equivalent Poisson’s ratio.The calculated creep force and creep torque between the wheel and rail needs to be converted to the absolute coordinate system. The load array of the wheel-rail creep force is expressed as(40)FcwrLi=FcLxi,FcLyi,FcLziT=CwrLi⋅FxLi,FyLi,0T,FcwrRi=FcRxi,FcRyi,FcRziT=CwrRi⋅FxRi,FyRi,0T,McwrLi=McLxi,McLyi,McLziT=CwrLi⋅0,0,MzLiT,McwrRi=McRxi,McRyi,McRziT=CwrRi⋅0,0,MzRiT,where(41)CwrLi=cosφwi−cosδcL+ψwisinφwisinδcL+ψwisinφwisinφwicosδcL+ψwicosφwi−sinδcL+ψwicosφwi0sinδcL+ψwicosδcL+ψwi,CwrRi=cosφwi−cosδcR−ψwisinφwi−sinδcR−ψwisinφwisinφwicosδcR−ψwicosφwisinδcR−ψwicosφwi0−sinδcR−ψwicosδcR−ψwi.(3) Geometric Model of Wheel-Rail Contact. Various parameters in the wheel-rail contact-force calculation, namely, the contact rail curvature radius ρr, the wheel tread curvature radius ρw, the wheel rolling radius rw, the wheel-rail contact angle δc, the wheel-rail contact point position C, and the relative compression of the wheel-rail contact εz must be solved by the geometric model of wheel-rail contact. The wheel-rail contact geometry is a spatially dynamic problem that accounts for the freedoms of the wheelset bounce ψw, wheelset yaw φw, and rail bounce ψr. This model discretizes the rail and wheel tread firstly and solves the tread curvature and tangent angle at the discrete points. The contact point of the wheel and rail is considered to lie on a spatial trace line, which is calculated as follows [3]:(42)xwL,R′=xwlx+lxRwtanδL,R,ywL,R′=xwly∓Rw2−lx2lx2lytanδL,R±lzm+Yw,zwL,R′=xwlz∓Rw2−lx2lx2lztanδL,R∓lzm+Zw,where m=1−lx21+tan2δL,R, and lx, ly, and lz are the X, Y, and Z direction cosines, respectively: lx=−cosψwsinφw, ly=cosψwsinφw, and lz=sinψw·dw is the rolling circular abscissa of the wheel tread.The random track random is regarded as the geometric displacement of the rail. Considering the transverse, vertical, and torsional degrees of freedom of the rail and the track irregularity, the discrete coordinate transformation of the rail is expressed as follows:(43)yrL,R′=yrL,RcosψrL±ψ0+zrL,RsinψrL±ψ0+YrL,R+ηyR,L,zrL,R′=−yrL,RsinψrL±ψ0+zrL,RcosψrL±ψ0+ZrL,R+ηzR,L,where, ηyR,L and ηzR,L are the vertical profiles and alignment irregularity of the left and right rail, respectively.As shown in Figure10, the three-dimensional curve of the wheel tread (xrL,R′, yrL,R′, zrL,R′) is projected onto the Y-O-Z plane and shifted upward by a fixed distance εz0. Interpolating between the discrete coordinates of the wheel tread (yzrL,R′, zyrL,R′) and the discrete coordinates of the rail (yzrL,R′, zyrL,R′), the model finds the point C that minimizes the wheel-to-rail distance in the Y-O-Z plane. The contact rail curvature radius ρr, wheel tread curvature radius ρw, wheel rolling radius rw, and wheel-rail contact angle δc are then obtained by interpolating between the coordinates of the C points. The normal compression of the wheel and rail is given by(44)εnL,R=εzL,RcosδcL,R±ψw.Figure 10
Comparison of track irregularity samples. (a) Vertical profiles of track irregularities, (b) track alignment irregularities, (c) track cross-level irregularities, and (d) track gage irregularities.
(a)
(b)
(c)
(d)Finally, the normal force between the wheel and rail is obtained by inserting the wheel-rail normal compression into the wheel-rail normal contact model.


### 2.3. Numerical Simulation and Verification of the Random Track Irregularity Model

In general, when simulating the random track irregularity time-domain samples, they often need dozens or hundreds of random variables. If the number theory method is used to select the random variables, it will inevitably produce certain aggregations, especially when high-dimensional selection is performed. To accurately reflect the entire random variable space, it is necessary to distribute the random variable points in the random variable space as much as possible. The generation of the clustering phenomenon will increase the number of random samples required for the calculation to some extent. Therefore, the random track irregularities time-domain samples were simulated by the spectral representation random-function method, which reduces the required number of random variables in the simulation. In a one-dimensional, univariate, stochastic process of the power spectral density functionS (ω), η(t) is expressed as [29, 30].(45)ηt=∫0∞cosωtdUtω+sinωtdVtω,where Ut(ω) and Vt(ω) are the spectral processes of a random process X(t), and their increments dUt(ω) and dVt(ω) satisfy the basic conditions of the spectral representation of a stochastic process:(46)EdUtω=EdVtω=0,ω≥0,EdUt2ω=EdVt2ω=2SXt,ωdω,ω≥0,EdUtωdUtω′=EdVtωdVtω′=0,ω,ω′≥0,ω′≠ω,EdUtωdVtω′=0,ω,ω′≥0Discretizing the above formula and introducing a set of standard orthogonal random variables {Xk, Yk}, η(t) is approximately discretized as(47)ηkt≈∑k=0∞cosωktΔUωk+sinωktΔVωk≈∑k=1NωcosωktσkXk+sinωktσkYk,where Δω is assumed sufficiently small, and dUt(ω) ≈ ΔUt(ω) = σkXk, dVt(ω) ≈ ΔVt(ω) = σkYk, and σk = [2SX(ωk) · Δω]1/2. It is easily proved that ΔUt(ω) and ΔVt(ω) satisfy the above basic conditions. Nω is the number of discrete frequencies. When the frequency components are uniformly dispersed, then ωk = ωmin + k · Δω, with Δω = (ωu − ωl)/N · ωu and ωl as the upper and lower limits, respectively, of the discrete frequency points.The orthogonal random variables {Xk, Yk} are obtained from the standard orthogonal variables {X¯n,Y¯n} by random mapping (as described in [30]). The standard orthogonal random variables {X¯n,Y¯n} can be constructed from orthogonal Legendre polynomials or Hartley orthogonal basis functions. The standard orthogonal variables {X¯n,Y¯n} are represented as(48)X¯n=casnΘ1,Y¯n=casnΘ2,n=1,2,…,N,where cas(x) denotes the Hartley orthogonal basis function, cas(x) = cos(x) + sin(x). Θ1 and Θ2 are independent random variables uniformly distributed in the interval [0, 2π].As an example, we consider the “TB/T3352-2014 high-speed railway ballastless track irregularity spectrum” implemented in 2015. Three hundred samples of random track irregularities were generated by the spectral representation random-function method. When generating the power spectrum of the track irregularities, the unit length and sampling interval were set to 1024 m and 0.25 m, respectively, with 4096 points in a single unit [31]. Therefore, when simulating the samples in the time domain, the discretized number of frequencies Nω was 4096. The vertical profile, alignment, cross-level, and gage irregularity of the track were, respectively, expressed as follows:Track vertical profile irregularity:(49)ηVit=2SVωkΔω⋅∑k=1NωcosωktX1k+sinωktX2k.Alignment irregularity:(50)ηAit=2SAωkΔω⋅∑k=1NωcosωktX3k+sinωktX4k.Cross-level irregularity:(51)ηCit=2SCωkΔω⋅∑k=1NωcosωktX5k+sinωktX6k.Gage irregularity:(52)ηGit=2SGωkΔω⋅∑k=1NωcosωktX7k+sinωktX8k.It should be noted that the track irregularity spectrum is averaged over many measured data. As the random error between a single sample and the track irregularity spectrum is 100%, the accuracy of the simulation method cannot be validated by comparing the simulated power spectrum of the track irregularity with that obtained from a single or finite-length sample, but validated by comparing the simulated power spectrum of the track irregularity with that obtained from a sample set.Figure10 compares the means and standard deviations of the track irregularity in the time domain and the track irregularity spectrum in [31]. The power spectrum of each frequency point in the sample set was tested by χ2 distribution with a significance level of α = 0.05 and two degrees of freedom with the Kolmogorov–Smirnov method, verifying the accuracy of the random track irregularities simulated by the present method.The power spectral density and track irregularity power spectrum derived from a single sample did not accurately estimate the vertical profile, alignment, cross-level, or gage irregularity of the track (Figure10). However, the mean power spectrum calculated from the sample set of track irregularities favorably agreed with the target spectrum. The maximum difference between the mean vertical profile, alignment, cross-level, and gage irregularity and those obtained from the target spectrum were 0.18 mm, 0.10 mm, 0.07 mm, and 0.08 mm, respectively, with maximum standard deviations of 4.68%, 5.77%, 6.97%, and 6.99%, respectively. In contrast, the maximum deviations in the power spectrum of the sample set were 14.54%, 10.23%, 10.07%, and 9.87%, respectively. In the power spectral densities at different frequencies of the sample set, the vertical profile, alignment, cross-level, and gage irregularity of the track all satisfied the 2-degree-of-freedom χ2 distribution (with 98.00%, 96.33%, 98.00%, and 98.33% of the frequency points satisfying the distribution, respectively). Each index of the sample set well agreed with the target spectrum, and only eight random variables were required to extract the indices from the time-domain samples. The method also retained a large number of frequency components. The commonly used trigonometric method, when simulating track random irregularity time-domain samples, requires enough random phase. When there are only 8 random phases, the simulated track random irregularity time-domain samples will have periodic repeatability and cannot be used in the calculation of the vehicle-rail coupled dynamic system. It can be seen that the numerical simulation method in this paper can greatly reduce the number of random variables required to simulate track random irregularity time-domain samples.


### 2.4. Numerical Solution Method

The track-vehicle coupled dynamic model was solved by an explicit-implicit integration method, with time integral step less than 0.0001 seconds, and the mode number of rail is taken as half the number of rail support points. The displacement and speed of the train and rail at the (n + 1)th moment were computed from the displacement, speed, and acceleration in the nth and (n − 1)th moments. The wheel-rail force was then calculated by the wheel-rail contact model, using the freedoms of the wheelset and rail at the (n + 1)th moment. After forming the wheel-rail load matrix, the fastener model was solved with the rail and slab freedoms, obtaining the wheel-rail load matrix. The vehicle and rail acceleration at the (n + 1)th moment was solved by the vehicle and rail motion equation, and the dynamic response of the slab and base at the (n + 1)th moment under the fastener force load was obtained by the Newmark-β integral method. Finally, the generalized probability density was evolved by the finite difference discretization of the modified Lax–Wendroff scheme [16].


## 2.1. Stochastic Analysis Method Based on Generalized Probability Density Evolution Theory

### 2.1.1. Generalized Probability Density Evolution Theory

Under random irregular excitations, the nonlinear vehicle-track coupled dynamic system is expressed as(1)Mξ¨t+Cξ˙t+Kξt=FΘ,t,where t denotes time, and ξt, ξ˙t, ξ¨t are the displacement, velocity, and acceleration vectors of the coupled system, respectively. The matrices [M], [C], and [K] are the mass, damping, and stiffness matrices, respectively, reflecting the mechanical and geometric nonlinearity in the wheel-rail relationship. F(Θ, t) is the input excitation vector of track random irregularities, and Θ is the set of sampled random variables.Most of the randomness in a vehicle-track coupled dynamic system derives from the random track irregularities. A nonlinear vehicle-track dynamic systemξ (t) is a probabilistic conservative stochastic system that excludes other random sources. According to the probabilistic conservation principle of stochastic systems, probabilistic conservation-augmented systems are constructed as [ξ (t), Θ] with(2)Pξt,Θ∈Ωt×Ωθ=Pξt0,Θ∈Ωt0×Ωθ,where P{·} denotes the probability of random events. Ωt and Ωt0 are the time distribution regions of the system state at times t and t0, respectively, and Ωθ is the space distribution region of the system state of sample θ. The above expression is equivalent to(3)ddtPξt,Θ∈Ωt×Ωθ=0.Expressing the joint probability density function (ξ(t), Θ) as pξθξ,θ,t, the above expression becomes(4)DDt∫Ωt×ΩθpξΘξ,θ,tdξdθ=0,where D/Dt is the derivative of matter [16].The generalized probability density evolves by the following equation [16]:(5)∂pξΘξ,θ,t∂t+ξ˙θ,t∂pξΘξ,θ,t∂ξ=0.In this formula,ξ˙θ,t is the free velocity of ξ (t) under the random sample {Θ = θ}. Its initial probability is PξΘξ,θ,tt=0=δξ−ξ0PΘθ.Solving the above evolution equation of generalized probability density under the initial conditions, we obtain the joint probability density function of sampleθ, pξΘξ,θ,t. Furthermore, the probability density function of PξΘ is obtained by integrating θ:(6)PξΘξ,t=∫ΩΘpξΘξ,θ,tdθ.For the probabilistic conservative vehicle-track coupled dynamic system under the random track irregularity excitation, the system dynamic responseξ satisfies the generalized probability density evolution equation shown in equation (5). Simultaneously, the generalized probability density evolution equation decouples the probability space from the physical space. As shown in Figure 1, first solve the vehicle-track coupled dynamic system equation (1) to obtain ξ˙θ,t under each random sample θ, and then solve the generalized probability density evolution equation for each random sample θ. When all random samples have been calculated, by solving equation (6), the probability density evolution process of the dynamic response can be obtained. Finally, a stochastic analysis of the nonlinear vehicle-track coupled dynamic system is realized.


### 2.1.2. Selection of Random Variable Samples and the Probability Assignment

In generalized probability density evolution theory, selecting the random variable samples is directly related to the computational efficiency of the random analysis. In this paper, the random variables are described by a low-deviation sequence selected by number theory, which reduces the required number of random samples and accelerates the convergence of the calculation. However, as the low-deviation sequence is uniformly distributed, it cannot be selected from a high-dimensional space without forming void and cluster structures. Therefore, the uniformity will degrade and the calculation results will contain large errors.Random irregularity is the superposition of irregular waves of different wavelengths, different phases, and different amplitudes. To accurately simulate the random track irregularities, more random frequencies and phase components (i.e., more selections) are needed. To avoid clustering in the random sample selection, the random track irregularities were simulated in the time domain by a stochastic function-spectral representation method, which requires just eight random variables. The generation method and validation of the sampling will be discussed in Section2.4. The random variables were sampled as Sobol sequences.The random samples in the probabilistic space ΩΘ must then be divided. Probabilistic spaces are commonly divided into representative Voronoi regions with representative volumes V (θk) satisfying(7)⋃kNptVθk=ΩΘ,Vθi∩Vθj,i≠j.The assigned probabilityPk of random sample θk is given by(8)Pk=∫VθkpΘθkdv,where pΘθk is the probability density of random variable Θ in Voronoi region V(θk). The joint probability density in the total partitioned space satisfies(9)∑k=1NptPk=∑k=1Npt∫VθkpΘθkdv=1.This method obtains the initial probability of each random point sample.


## 2.1.1. Generalized Probability Density Evolution Theory

Under random irregular excitations, the nonlinear vehicle-track coupled dynamic system is expressed as(1)Mξ¨t+Cξ˙t+Kξt=FΘ,t,where t denotes time, and ξt, ξ˙t, ξ¨t are the displacement, velocity, and acceleration vectors of the coupled system, respectively. The matrices [M], [C], and [K] are the mass, damping, and stiffness matrices, respectively, reflecting the mechanical and geometric nonlinearity in the wheel-rail relationship. F(Θ, t) is the input excitation vector of track random irregularities, and Θ is the set of sampled random variables.Most of the randomness in a vehicle-track coupled dynamic system derives from the random track irregularities. A nonlinear vehicle-track dynamic systemξ (t) is a probabilistic conservative stochastic system that excludes other random sources. According to the probabilistic conservation principle of stochastic systems, probabilistic conservation-augmented systems are constructed as [ξ (t), Θ] with(2)Pξt,Θ∈Ωt×Ωθ=Pξt0,Θ∈Ωt0×Ωθ,where P{·} denotes the probability of random events. Ωt and Ωt0 are the time distribution regions of the system state at times t and t0, respectively, and Ωθ is the space distribution region of the system state of sample θ. The above expression is equivalent to(3)ddtPξt,Θ∈Ωt×Ωθ=0.Expressing the joint probability density function (ξ(t), Θ) as pξθξ,θ,t, the above expression becomes(4)DDt∫Ωt×ΩθpξΘξ,θ,tdξdθ=0,where D/Dt is the derivative of matter [16].The generalized probability density evolves by the following equation [16]:(5)∂pξΘξ,θ,t∂t+ξ˙θ,t∂pξΘξ,θ,t∂ξ=0.In this formula,ξ˙θ,t is the free velocity of ξ (t) under the random sample {Θ = θ}. Its initial probability is PξΘξ,θ,tt=0=δξ−ξ0PΘθ.Solving the above evolution equation of generalized probability density under the initial conditions, we obtain the joint probability density function of sampleθ, pξΘξ,θ,t. Furthermore, the probability density function of PξΘ is obtained by integrating θ:(6)PξΘξ,t=∫ΩΘpξΘξ,θ,tdθ.For the probabilistic conservative vehicle-track coupled dynamic system under the random track irregularity excitation, the system dynamic responseξ satisfies the generalized probability density evolution equation shown in equation (5). Simultaneously, the generalized probability density evolution equation decouples the probability space from the physical space. As shown in Figure 1, first solve the vehicle-track coupled dynamic system equation (1) to obtain ξ˙θ,t under each random sample θ, and then solve the generalized probability density evolution equation for each random sample θ. When all random samples have been calculated, by solving equation (6), the probability density evolution process of the dynamic response can be obtained. Finally, a stochastic analysis of the nonlinear vehicle-track coupled dynamic system is realized.


## 2.1.2. Selection of Random Variable Samples and the Probability Assignment

In generalized probability density evolution theory, selecting the random variable samples is directly related to the computational efficiency of the random analysis. In this paper, the random variables are described by a low-deviation sequence selected by number theory, which reduces the required number of random samples and accelerates the convergence of the calculation. However, as the low-deviation sequence is uniformly distributed, it cannot be selected from a high-dimensional space without forming void and cluster structures. Therefore, the uniformity will degrade and the calculation results will contain large errors.Random irregularity is the superposition of irregular waves of different wavelengths, different phases, and different amplitudes. To accurately simulate the random track irregularities, more random frequencies and phase components (i.e., more selections) are needed. To avoid clustering in the random sample selection, the random track irregularities were simulated in the time domain by a stochastic function-spectral representation method, which requires just eight random variables. The generation method and validation of the sampling will be discussed in Section2.4. The random variables were sampled as Sobol sequences.The random samples in the probabilistic space ΩΘ must then be divided. Probabilistic spaces are commonly divided into representative Voronoi regions with representative volumes V (θk) satisfying(7)⋃kNptVθk=ΩΘ,Vθi∩Vθj,i≠j.The assigned probabilityPk of random sample θk is given by(8)Pk=∫VθkpΘθkdv,where pΘθk is the probability density of random variable Θ in Voronoi region V(θk). The joint probability density in the total partitioned space satisfies(9)∑k=1NptPk=∑k=1Npt∫VθkpΘθkdv=1.This method obtains the initial probability of each random point sample.


## 2.2. Vehicle-Track Coupled Dynamics Model

Figure2 shows the vehicle-track coupling dynamic model. The vehicle is modeled by multirigid body dynamics, and the rail is regarded as a finite-length, elastic-point supported Euler beam. The dynamic vibration equation of the beam is solved by the Ritz method. The wheel-rail relationship is modeled by nonlinear Hertz contact theory and wheel-rail nonlinear creep theory. The fastener is modeled by a linear spring damping element in the transverse and vertical directions. The track structure is modeled by the finite element method; the slab, mortar layer, and base are modeled as plate elements. The foundation support is simplified as a linear spring damping surface element. The wheel-rail relationship connects the vehicle and rail models, whereas the fastener force connects the rail model and the track slab-foundation structure.Figure 2
Vehicle-track coupled dynamics model.


### 2.2.1. Vehicle Dynamics Model

Figure3 schematizes the vehicle model composed of one car body, two bogies, and four wheelsets. Each part is treated as a rigid body with freedom in the longitudinal, transverse, bounce, roll, and yaw directions (giving 35 degrees of freedom per vehicle).Figure 3
Vehicle dynamics model.The antisnake shock absorber significantly affects the driving quality when the train is running at a high speed. Treating the force of the antisnake shock absorber on the vehicle and the antirolling moment of the vehicle’s central suspension as an additional force on the vehicle, the dynamic equation of the vehicle is given by(10)MVξ¨V+CVξ˙V+KVξV+GV=FV,where [MV], [CV], [KV], [GV], [FV] are the vehicle mass, damping, stiffness, additional force, and load matrices, respectively.The freedom matrix of the vehicle is expressed as(11)ξV=ξc,ξb1,ξb2,ξw1,ξw2,ξw3,ξw4T,where the subscripts c, b, and w represent the car body, bogie, and wheelset, respectively, and the subscripts i and j represent the ith bogie and the jth wheelset, respectively. The degrees of freedom of each part are expressed as(12)ξc,bi,wj=Yc,bi,wj,Zc,bi,wj,ψc,bi,wj,βc,bi,wj,φc,bi,wj,where Y, Z, Ψ, φ, and β represent the freedoms in the longitudinal, transverse, bounce, roll, and yaw, respectively.The vehicle-mass matrix is given by(13)MV=Mw1,Mw2,Mw3,Mw4,Mt1,Mt2,Mc,where Mwi,bj,c=diagmwi,bj,c,mwi,bj,c,Ixwi,bj,c,Iywi,bj,c,Izwi,bj,c. Here m denotes the mass, and Ix, Iy, and Iz are the moments of inertia around the X, Y, and Z axes, respectively.The vehicle-stiffness matrix is given by(14)KV=Kw1000Kwb100Kw200Kwb200Kw300Kwb30Kw40Kwb40Kb10Kbc1Kb2Kbc2symsKc,where(15)Kwi=2kpykpzkpzdw20kpxdw2,i=1,2,3,4,Kbj=22kpy+ksy0−2kpyHtw+ksyHbt002kpz+ksz0002kpyHtw2+2kpzdwk2+kszdsk2+ksyHbt2002kpzlt2+2kpxHtw2+ksxHbt20syms2kpylt2+2kpxdwk2+ksxdsk2,Kc=4ksy0−ksyHcb00ksz000ksyHcb2+kszdsk200kszlc2+ksxHcb20symsksylc2+ksxdsk2,Kwbi=2−kpy0kpyHtw0−1ikpylt−kpz0−1i+1kpzlt0−kpzdwk20000syms−kpxdwk2,i=1,2,3,4,Kbcj=2−ksy0ksyHcb0−1j+1ksylc−ksz0−1j+1kszlc0−kszdsk2+ksyHcbHbt0−1jksylcHbtksxHcbHbt0syms−ksxdsk2,j=1,2.The vehicle-damping matrix(16)CV=Cw1000Cwb100Cw200Cwb200Cw300Cwb30Cw40Cwb40Cb10Cbc1Cb2Cbc2symsCc,is similar to the vehicle-stiffness matrix. In fact, the vehicle-damping matrix is obtained by replacing the stiffness values in the stiffness matrix with the corresponding damping values.The additional moment on the vehicle is expressed as the matrix(17)G=Gw1,Gw2,Gw3,Gw4,Gt1,Gt2,Gc,with(18)Gwi=0,0,Iywiβ˙wi−Ωφ˙wi,0,Iywiβ˙wi−Ωψ˙wiT,Gtj=0,0,0,FxsLj+FxsRjHbt,FxsLj−FxsRjdsxT,Gc=0,0,Mr1+Mr2,FxsL1+FxsR1+FxsL2+FxsR2Hcb,FxsL1+FxsL2−FxsR1−FxsR2dsx.In the above expressions,Gwi is the matrix of moments acting on the ith wheelset caused by wheelset rolling rotations, Gtj is the matrix of moments acting on the jth bioge caused by absorber, Gc is the matrix of moments acting on the car body caused by absorber, and FxsLj and FxsLj are the additional forces imposed by the snake-shaped shock absorber. They are calculated as follows:(19)FxsL,Rj=FmaxvxctL,Rjv0,vxctL,Rj<v0,FmaxsignvxctL,Rj,vxctL,Rj≥v0.M
ri is the antiroll moment, given by(20)Mri=Krxψc−ψti,i=1,2.The vehicle-load matrix is expressed as(21)F=A1Fwr1,A2Fwr2,A3Fwr3,A4Fwr4,05×1,05×1,05×1,where Ai is the matrix of conversion coefficients between the wheel-rail force and the vehicle freedoms. Fwri is the wheel-rail force array, with(22)Awi=01010000000000−100−10000000−rLia00−rLi−a0000000rLirLiφwi0rRirRiφwi0010010a0a0φwi0−a0−a0φwi0001001,Fwri=FLxi,FLyi,FLzi,FRxi,FRyi,FRzi,McLxi,McLyi,McLzi,McRxi,McRyi,McRzi,where FL,Rxi = FnL,Rxi + FcL,Rxi, FL,Ryi = FnL,Ryi + FcL,Ryi, and FL,Rzi = FnL,Rzi + FcL,Rzi. Fn and Fc are the normal and creep forces, respectively, between the wheel and the rail. The subscripts L and R represent the left and right sides, respectively, and the subscripts x, y, and z denote the upper component in the X-, Y-, and Z-axis directions, respectively.


### 2.2.2. Dynamic Model of the Track Structure

The track structure is divided into two main parts: the rail and the track plate. The rail is treated as an elastic support for an Euler beam, and the track plate and base plate are modeled by a plate element with high rigidity in the transverse direction, lateral displacement using beam element to simulate. Since the numerical integration solution step length is short and requires a lot of numerical calculation work, the method of the literature [3] is used in the modeling of the rail, and the dynamic response of the rail is obtained by solving the rail vibration differential equation, which can greatly reduce the calculation work. At the same time, the track plate is modeled by the finite element method to facilitate the modeling of various ballastless track structures. The track-structure model and rail forces are shown in Figures 4 and 5, respectively.Figure 4
Structural model of the rail-track plate.Figure 5
Rail force diagram.The following expressions are shown for the right rail (the left rail is treated analogously). The fastener force matrix is expressed as(23)FpviFphiMpi=0kpv00−kpvkph0kphHra−kph0−kphHra0kphHra2+kpvLr2kphHra0Yri,Zri,ψri,Ysi,ZsiT.The torque of the wheel-rail force, which induces rail torsion, is(24)Mcr=FRyHrb−FRzLce.The differential equations of the vertical, transverse, and torsional vibrations of the rail are, respectively given in the following [3].Vertical:(25)EIrz∂4Yrx,t∂x4+mr∂2Yrx,t∂t2=−∑i=1NFpyitδx−xi+∑j=14FRyjtδx−xwj.Transverse:(26)ErIrz∂4Zrx,t∂x4+mr∂2Zrx,t∂t2=−∑i=1NFpvitδx−xi+∑i=14FRzjtδx−xwj.Torsional:(27)ρrIr0∂2ψrx,t∂t4−GrIt∂2ψrx,t∂x2=−∑i=1NMpitδx−xi+∑j=14Mcrjtδx−xwj.The track and base plates are considered as horizontally rigid bodies and are modeled as four-node bending plate elements. They are connected by a fastener element, which is simulated by a linear spring damping element. Slab-base has been modeled by the principle of total potential energy with stationary value in elastic system dynamics [4], the dynamics of the slab and base plates are governed by(28)MssMbbδ¨sδ¨b+CssCsbCbsCbbδ˙sδ˙b+KssKsbKbsKbbδsδb=FsFb,where the subscript s and b denote the track and base plates, respectively. The vibration equation of the track and base plates assumes the principle of invariant potential energy of an elastic system. The foundation structure is modeled by the same method. [F] is the load matrix, and the main load acting on the slab-base is the fastening force. The matrix and modeling process of the track and base plates are detailed elsewhere [13].


### 2.2.3. Wheel-Rail Contact Relation

The wheel-rail relationship links the vehicle system to the track system. In the vehicle-track coupled dynamics, the wheel-rail relationship mainly includes the wheel-rail normal-seeking solution, the wheel-rail creep solution, and the wheel-rail contact geometry. The wheel-rail relationship is shown in Figure6.Figure 6
Wheel-rail contact model.(1) The Wheel-Rail Normal Contact Model. The wheel tread and rail profile are complex curves composed of many simple curves with different radii. According to nonlinear Hertz contact theory (Figure 7), the shapes of the wheel and rail tread near the contact point are approximately described as(29)z1=12ρwx12+12rwy12,for the wheel tread,z2=−12ρrx22,for the rail tread,where rp, rw, and rr represent the curvature radii of the wheel tread, wheel rolling circle, and rail tread, respectively. Positive and negative curvatures indicate a convex and concave surface, respectively.The pressure distribution on the contact surface is given by Hertz contact theory.(30)p=Pn1−xa2−yb21/2.Figure 7
Simplified schematic of wheel-rail contact.The normal force on the wheel and rail is expressed as(31)FnL,R=43E∗Re1/2εnL,RF2e3/2.And the maximum pressure between the wheel and rail is given by(32)PnL,R=6FnL,RE∗2π3Re21/3F1e−2.In the above expressions,E∗ is the equivalent modulus of elasticity, given by 1/E∗=1−vr2/Er+1−vw2/Ew, and e=1−b2/a21/2, where a and b are the lengths of the long and short axes, respectively, of the contact ellipse (with b < a). F1(e) and F2(e) are correction factors that depend only on the ratio of relative curvature radii R′/R″. The variations of these factors with increasing R′/R″ are plotted in Figure 8.Figure 8
Correction factorsF1(e) and (F)2(e) versus ratio of relative curvature radius.F
1(e) and F2(e) can be directly computed by a formula or directly read from Figure 9. To improve the efficiency of the computation, the curve-fitting parameters obtained from Figure 9 are listed in Table 1, F1(e) and F2(e) were then obtained as(33)Fie=q1x5+q2x4+q3x3+q4x2+q5x1+q6x3+p1x2+q2x1+q3,with x = ln(B/A).Figure 9
Geometric wheel-rail contact model.Table 1
Fitting parameters of theF1(e) and F2(e) curves.
F
1(e)
F
2(e)
Parameters
Error of fitting
Parameters
Error of fitting
q
1
0.001227
SSE: 3.522e-07
q
1
0.002195
SSE: 3.156e-07
q
2
−0.05683
R-square: 1
q
2
−0.07143
R-square: 1
q
3
0.7315
Adjusted R-square: 1
q
3
0.4331
Adjusted R-square: 1
q
4
4.087
RMSE: 5.935e-07
q
4
4.816
RMSE: 5.618e-07
q
5
20.1
q
5
41.2
q
6
36.46
q
6
21.94
p
1
4.599
p
1
5.123
p
2
20.1
p
2
41.2
p
3
36.46
p
3
21.94When calculating the normal contact force between the wheel and rail, the relative compression between the wheel and a rail with random irregularities is obtained by geometrically solving the wheel-rail contact. The normal force between the wheel and rail is then obtained by solving the wheel-rail normal contact model and determining the direction and magnitude variations in the normal wheel-rail force at different contact points. Converting to the absolute coordinate system, the normal loads on the wheel and rail array are, respectively, expressed as(34)FnLi′=FnLxi,FnLyi,FnLziT=FnLisinδcLi+ψwisinφwi,−sinδcLi+ψwicosφwi,cosδcLi+ψwiT,FnRi′=FnRxi,FnRyi,FnRziT=FnRi−sinδcRi−ψwisinφwi,sinδcRi−ψwicosφwi,cosδcRi−ψwiT.Here,δcLi and δcRi are the left-sided and right-sided wheel-rail contact angles, respectively, obtained from the wheel-rail contact geometry (details will be given in Section (2)).(2) Wheel-Rail Tangential Contact Model. The wheel-rail tangential contact was solved by Kalker linear creep theory [27], corrected for nonlinearity. The creep force and creep torque between the wheel and rail after the nonlinear correction are given by(35)Fx′=α⋅Fx,Fy′=α⋅Fy,Mz′=α⋅Mz,where α is the correction coefficient, given as α=FR′/FR with FR=Fx2+Fy2 and(36)FR′=μFnL,RFRμFnL,R−13FRμFnL,R2+127FRμFnL,R3,FR≤3μFnL,R,μFnL,R,FR>3μFnL,R,Fx and Fy denote the creep forces in the X- and Y-axis directions, respectively, and Mz is the rotational creep moment. In linear creep theory, these components are, respectively, expressed as(37)Fx=−f11ξx,Fy=−f22ξy−f23ξsp,Mz=f23ξy−f23ξsp,where f11 and f22 are the longitudinal and transverse creep coefficients, respectively, and f23 and f33 are rotational creep coefficients:(38)f11=GwrabC11,f22=GwrabC22,f23=Gwrab3/2C23,f33=Gwrab2C33.In the above expressions,Gwr is the composite shear modulus of the wheel and rail materials, and Cij are the Kalker coefficients [27] (after Kalker, who calculated them). Meymandet al. [28] have fitted the Kalker coefficients under different contact conditions. Here, the Cijs are expressed as(39)C11=2.39+2.338vwr+1.014g+2.736vwr2−0.0728vwrg−0.0139g2,C22=2.353−0.0212vwr+1.074g+1.163vwrg−0.0319g2,C23=0.2866+0.4467vwr+1.049g+0.038vwrg−0.0177g2+0.005g3+0.1479g2vwr,C33=0.7375−1.177vwr+0.45061g+1.076vwr1g−0.0231g2,where g is the ratio of the long- and short-elliptical axes, g = a/b, and vwr is the equivalent Poisson’s ratio.The calculated creep force and creep torque between the wheel and rail needs to be converted to the absolute coordinate system. The load array of the wheel-rail creep force is expressed as(40)FcwrLi=FcLxi,FcLyi,FcLziT=CwrLi⋅FxLi,FyLi,0T,FcwrRi=FcRxi,FcRyi,FcRziT=CwrRi⋅FxRi,FyRi,0T,McwrLi=McLxi,McLyi,McLziT=CwrLi⋅0,0,MzLiT,McwrRi=McRxi,McRyi,McRziT=CwrRi⋅0,0,MzRiT,where(41)CwrLi=cosφwi−cosδcL+ψwisinφwisinδcL+ψwisinφwisinφwicosδcL+ψwicosφwi−sinδcL+ψwicosφwi0sinδcL+ψwicosδcL+ψwi,CwrRi=cosφwi−cosδcR−ψwisinφwi−sinδcR−ψwisinφwisinφwicosδcR−ψwicosφwisinδcR−ψwicosφwi0−sinδcR−ψwicosδcR−ψwi.(3) Geometric Model of Wheel-Rail Contact. Various parameters in the wheel-rail contact-force calculation, namely, the contact rail curvature radius ρr, the wheel tread curvature radius ρw, the wheel rolling radius rw, the wheel-rail contact angle δc, the wheel-rail contact point position C, and the relative compression of the wheel-rail contact εz must be solved by the geometric model of wheel-rail contact. The wheel-rail contact geometry is a spatially dynamic problem that accounts for the freedoms of the wheelset bounce ψw, wheelset yaw φw, and rail bounce ψr. This model discretizes the rail and wheel tread firstly and solves the tread curvature and tangent angle at the discrete points. The contact point of the wheel and rail is considered to lie on a spatial trace line, which is calculated as follows [3]:(42)xwL,R′=xwlx+lxRwtanδL,R,ywL,R′=xwly∓Rw2−lx2lx2lytanδL,R±lzm+Yw,zwL,R′=xwlz∓Rw2−lx2lx2lztanδL,R∓lzm+Zw,where m=1−lx21+tan2δL,R, and lx, ly, and lz are the X, Y, and Z direction cosines, respectively: lx=−cosψwsinφw, ly=cosψwsinφw, and lz=sinψw·dw is the rolling circular abscissa of the wheel tread.The random track random is regarded as the geometric displacement of the rail. Considering the transverse, vertical, and torsional degrees of freedom of the rail and the track irregularity, the discrete coordinate transformation of the rail is expressed as follows:(43)yrL,R′=yrL,RcosψrL±ψ0+zrL,RsinψrL±ψ0+YrL,R+ηyR,L,zrL,R′=−yrL,RsinψrL±ψ0+zrL,RcosψrL±ψ0+ZrL,R+ηzR,L,where, ηyR,L and ηzR,L are the vertical profiles and alignment irregularity of the left and right rail, respectively.As shown in Figure10, the three-dimensional curve of the wheel tread (xrL,R′, yrL,R′, zrL,R′) is projected onto the Y-O-Z plane and shifted upward by a fixed distance εz0. Interpolating between the discrete coordinates of the wheel tread (yzrL,R′, zyrL,R′) and the discrete coordinates of the rail (yzrL,R′, zyrL,R′), the model finds the point C that minimizes the wheel-to-rail distance in the Y-O-Z plane. The contact rail curvature radius ρr, wheel tread curvature radius ρw, wheel rolling radius rw, and wheel-rail contact angle δc are then obtained by interpolating between the coordinates of the C points. The normal compression of the wheel and rail is given by(44)εnL,R=εzL,RcosδcL,R±ψw.Figure 10
Comparison of track irregularity samples. (a) Vertical profiles of track irregularities, (b) track alignment irregularities, (c) track cross-level irregularities, and (d) track gage irregularities.
(a)
(b)
(c)
(d)Finally, the normal force between the wheel and rail is obtained by inserting the wheel-rail normal compression into the wheel-rail normal contact model.


## 2.2.1. Vehicle Dynamics Model

Figure3 schematizes the vehicle model composed of one car body, two bogies, and four wheelsets. Each part is treated as a rigid body with freedom in the longitudinal, transverse, bounce, roll, and yaw directions (giving 35 degrees of freedom per vehicle).Figure 3
Vehicle dynamics model.The antisnake shock absorber significantly affects the driving quality when the train is running at a high speed. Treating the force of the antisnake shock absorber on the vehicle and the antirolling moment of the vehicle’s central suspension as an additional force on the vehicle, the dynamic equation of the vehicle is given by(10)MVξ¨V+CVξ˙V+KVξV+GV=FV,where [MV], [CV], [KV], [GV], [FV] are the vehicle mass, damping, stiffness, additional force, and load matrices, respectively.The freedom matrix of the vehicle is expressed as(11)ξV=ξc,ξb1,ξb2,ξw1,ξw2,ξw3,ξw4T,where the subscripts c, b, and w represent the car body, bogie, and wheelset, respectively, and the subscripts i and j represent the ith bogie and the jth wheelset, respectively. The degrees of freedom of each part are expressed as(12)ξc,bi,wj=Yc,bi,wj,Zc,bi,wj,ψc,bi,wj,βc,bi,wj,φc,bi,wj,where Y, Z, Ψ, φ, and β represent the freedoms in the longitudinal, transverse, bounce, roll, and yaw, respectively.The vehicle-mass matrix is given by(13)MV=Mw1,Mw2,Mw3,Mw4,Mt1,Mt2,Mc,where Mwi,bj,c=diagmwi,bj,c,mwi,bj,c,Ixwi,bj,c,Iywi,bj,c,Izwi,bj,c. Here m denotes the mass, and Ix, Iy, and Iz are the moments of inertia around the X, Y, and Z axes, respectively.The vehicle-stiffness matrix is given by(14)KV=Kw1000Kwb100Kw200Kwb200Kw300Kwb30Kw40Kwb40Kb10Kbc1Kb2Kbc2symsKc,where(15)Kwi=2kpykpzkpzdw20kpxdw2,i=1,2,3,4,Kbj=22kpy+ksy0−2kpyHtw+ksyHbt002kpz+ksz0002kpyHtw2+2kpzdwk2+kszdsk2+ksyHbt2002kpzlt2+2kpxHtw2+ksxHbt20syms2kpylt2+2kpxdwk2+ksxdsk2,Kc=4ksy0−ksyHcb00ksz000ksyHcb2+kszdsk200kszlc2+ksxHcb20symsksylc2+ksxdsk2,Kwbi=2−kpy0kpyHtw0−1ikpylt−kpz0−1i+1kpzlt0−kpzdwk20000syms−kpxdwk2,i=1,2,3,4,Kbcj=2−ksy0ksyHcb0−1j+1ksylc−ksz0−1j+1kszlc0−kszdsk2+ksyHcbHbt0−1jksylcHbtksxHcbHbt0syms−ksxdsk2,j=1,2.The vehicle-damping matrix(16)CV=Cw1000Cwb100Cw200Cwb200Cw300Cwb30Cw40Cwb40Cb10Cbc1Cb2Cbc2symsCc,is similar to the vehicle-stiffness matrix. In fact, the vehicle-damping matrix is obtained by replacing the stiffness values in the stiffness matrix with the corresponding damping values.The additional moment on the vehicle is expressed as the matrix(17)G=Gw1,Gw2,Gw3,Gw4,Gt1,Gt2,Gc,with(18)Gwi=0,0,Iywiβ˙wi−Ωφ˙wi,0,Iywiβ˙wi−Ωψ˙wiT,Gtj=0,0,0,FxsLj+FxsRjHbt,FxsLj−FxsRjdsxT,Gc=0,0,Mr1+Mr2,FxsL1+FxsR1+FxsL2+FxsR2Hcb,FxsL1+FxsL2−FxsR1−FxsR2dsx.In the above expressions,Gwi is the matrix of moments acting on the ith wheelset caused by wheelset rolling rotations, Gtj is the matrix of moments acting on the jth bioge caused by absorber, Gc is the matrix of moments acting on the car body caused by absorber, and FxsLj and FxsLj are the additional forces imposed by the snake-shaped shock absorber. They are calculated as follows:(19)FxsL,Rj=FmaxvxctL,Rjv0,vxctL,Rj<v0,FmaxsignvxctL,Rj,vxctL,Rj≥v0.M
ri is the antiroll moment, given by(20)Mri=Krxψc−ψti,i=1,2.The vehicle-load matrix is expressed as(21)F=A1Fwr1,A2Fwr2,A3Fwr3,A4Fwr4,05×1,05×1,05×1,where Ai is the matrix of conversion coefficients between the wheel-rail force and the vehicle freedoms. Fwri is the wheel-rail force array, with(22)Awi=01010000000000−100−10000000−rLia00−rLi−a0000000rLirLiφwi0rRirRiφwi0010010a0a0φwi0−a0−a0φwi0001001,Fwri=FLxi,FLyi,FLzi,FRxi,FRyi,FRzi,McLxi,McLyi,McLzi,McRxi,McRyi,McRzi,where FL,Rxi = FnL,Rxi + FcL,Rxi, FL,Ryi = FnL,Ryi + FcL,Ryi, and FL,Rzi = FnL,Rzi + FcL,Rzi. Fn and Fc are the normal and creep forces, respectively, between the wheel and the rail. The subscripts L and R represent the left and right sides, respectively, and the subscripts x, y, and z denote the upper component in the X-, Y-, and Z-axis directions, respectively.


## 2.2.2. Dynamic Model of the Track Structure

The track structure is divided into two main parts: the rail and the track plate. The rail is treated as an elastic support for an Euler beam, and the track plate and base plate are modeled by a plate element with high rigidity in the transverse direction, lateral displacement using beam element to simulate. Since the numerical integration solution step length is short and requires a lot of numerical calculation work, the method of the literature [3] is used in the modeling of the rail, and the dynamic response of the rail is obtained by solving the rail vibration differential equation, which can greatly reduce the calculation work. At the same time, the track plate is modeled by the finite element method to facilitate the modeling of various ballastless track structures. The track-structure model and rail forces are shown in Figures 4 and 5, respectively.Figure 4
Structural model of the rail-track plate.Figure 5
Rail force diagram.The following expressions are shown for the right rail (the left rail is treated analogously). The fastener force matrix is expressed as(23)FpviFphiMpi=0kpv00−kpvkph0kphHra−kph0−kphHra0kphHra2+kpvLr2kphHra0Yri,Zri,ψri,Ysi,ZsiT.The torque of the wheel-rail force, which induces rail torsion, is(24)Mcr=FRyHrb−FRzLce.The differential equations of the vertical, transverse, and torsional vibrations of the rail are, respectively given in the following [3].Vertical:(25)EIrz∂4Yrx,t∂x4+mr∂2Yrx,t∂t2=−∑i=1NFpyitδx−xi+∑j=14FRyjtδx−xwj.Transverse:(26)ErIrz∂4Zrx,t∂x4+mr∂2Zrx,t∂t2=−∑i=1NFpvitδx−xi+∑i=14FRzjtδx−xwj.Torsional:(27)ρrIr0∂2ψrx,t∂t4−GrIt∂2ψrx,t∂x2=−∑i=1NMpitδx−xi+∑j=14Mcrjtδx−xwj.The track and base plates are considered as horizontally rigid bodies and are modeled as four-node bending plate elements. They are connected by a fastener element, which is simulated by a linear spring damping element. Slab-base has been modeled by the principle of total potential energy with stationary value in elastic system dynamics [4], the dynamics of the slab and base plates are governed by(28)MssMbbδ¨sδ¨b+CssCsbCbsCbbδ˙sδ˙b+KssKsbKbsKbbδsδb=FsFb,where the subscript s and b denote the track and base plates, respectively. The vibration equation of the track and base plates assumes the principle of invariant potential energy of an elastic system. The foundation structure is modeled by the same method. [F] is the load matrix, and the main load acting on the slab-base is the fastening force. The matrix and modeling process of the track and base plates are detailed elsewhere [13].


## 2.2.3. Wheel-Rail Contact Relation

The wheel-rail relationship links the vehicle system to the track system. In the vehicle-track coupled dynamics, the wheel-rail relationship mainly includes the wheel-rail normal-seeking solution, the wheel-rail creep solution, and the wheel-rail contact geometry. The wheel-rail relationship is shown in Figure6.Figure 6
Wheel-rail contact model.(1) The Wheel-Rail Normal Contact Model. The wheel tread and rail profile are complex curves composed of many simple curves with different radii. According to nonlinear Hertz contact theory (Figure 7), the shapes of the wheel and rail tread near the contact point are approximately described as(29)z1=12ρwx12+12rwy12,for the wheel tread,z2=−12ρrx22,for the rail tread,where rp, rw, and rr represent the curvature radii of the wheel tread, wheel rolling circle, and rail tread, respectively. Positive and negative curvatures indicate a convex and concave surface, respectively.The pressure distribution on the contact surface is given by Hertz contact theory.(30)p=Pn1−xa2−yb21/2.Figure 7
Simplified schematic of wheel-rail contact.The normal force on the wheel and rail is expressed as(31)FnL,R=43E∗Re1/2εnL,RF2e3/2.And the maximum pressure between the wheel and rail is given by(32)PnL,R=6FnL,RE∗2π3Re21/3F1e−2.In the above expressions,E∗ is the equivalent modulus of elasticity, given by 1/E∗=1−vr2/Er+1−vw2/Ew, and e=1−b2/a21/2, where a and b are the lengths of the long and short axes, respectively, of the contact ellipse (with b < a). F1(e) and F2(e) are correction factors that depend only on the ratio of relative curvature radii R′/R″. The variations of these factors with increasing R′/R″ are plotted in Figure 8.Figure 8
Correction factorsF1(e) and (F)2(e) versus ratio of relative curvature radius.F
1(e) and F2(e) can be directly computed by a formula or directly read from Figure 9. To improve the efficiency of the computation, the curve-fitting parameters obtained from Figure 9 are listed in Table 1, F1(e) and F2(e) were then obtained as(33)Fie=q1x5+q2x4+q3x3+q4x2+q5x1+q6x3+p1x2+q2x1+q3,with x = ln(B/A).Figure 9
Geometric wheel-rail contact model.Table 1
Fitting parameters of theF1(e) and F2(e) curves.
F
1(e)
F
2(e)
Parameters
Error of fitting
Parameters
Error of fitting
q
1
0.001227
SSE: 3.522e-07
q
1
0.002195
SSE: 3.156e-07
q
2
−0.05683
R-square: 1
q
2
−0.07143
R-square: 1
q
3
0.7315
Adjusted R-square: 1
q
3
0.4331
Adjusted R-square: 1
q
4
4.087
RMSE: 5.935e-07
q
4
4.816
RMSE: 5.618e-07
q
5
20.1
q
5
41.2
q
6
36.46
q
6
21.94
p
1
4.599
p
1
5.123
p
2
20.1
p
2
41.2
p
3
36.46
p
3
21.94When calculating the normal contact force between the wheel and rail, the relative compression between the wheel and a rail with random irregularities is obtained by geometrically solving the wheel-rail contact. The normal force between the wheel and rail is then obtained by solving the wheel-rail normal contact model and determining the direction and magnitude variations in the normal wheel-rail force at different contact points. Converting to the absolute coordinate system, the normal loads on the wheel and rail array are, respectively, expressed as(34)FnLi′=FnLxi,FnLyi,FnLziT=FnLisinδcLi+ψwisinφwi,−sinδcLi+ψwicosφwi,cosδcLi+ψwiT,FnRi′=FnRxi,FnRyi,FnRziT=FnRi−sinδcRi−ψwisinφwi,sinδcRi−ψwicosφwi,cosδcRi−ψwiT.Here,δcLi and δcRi are the left-sided and right-sided wheel-rail contact angles, respectively, obtained from the wheel-rail contact geometry (details will be given in Section (2)).(2) Wheel-Rail Tangential Contact Model. The wheel-rail tangential contact was solved by Kalker linear creep theory [27], corrected for nonlinearity. The creep force and creep torque between the wheel and rail after the nonlinear correction are given by(35)Fx′=α⋅Fx,Fy′=α⋅Fy,Mz′=α⋅Mz,where α is the correction coefficient, given as α=FR′/FR with FR=Fx2+Fy2 and(36)FR′=μFnL,RFRμFnL,R−13FRμFnL,R2+127FRμFnL,R3,FR≤3μFnL,R,μFnL,R,FR>3μFnL,R,Fx and Fy denote the creep forces in the X- and Y-axis directions, respectively, and Mz is the rotational creep moment. In linear creep theory, these components are, respectively, expressed as(37)Fx=−f11ξx,Fy=−f22ξy−f23ξsp,Mz=f23ξy−f23ξsp,where f11 and f22 are the longitudinal and transverse creep coefficients, respectively, and f23 and f33 are rotational creep coefficients:(38)f11=GwrabC11,f22=GwrabC22,f23=Gwrab3/2C23,f33=Gwrab2C33.In the above expressions,Gwr is the composite shear modulus of the wheel and rail materials, and Cij are the Kalker coefficients [27] (after Kalker, who calculated them). Meymandet al. [28] have fitted the Kalker coefficients under different contact conditions. Here, the Cijs are expressed as(39)C11=2.39+2.338vwr+1.014g+2.736vwr2−0.0728vwrg−0.0139g2,C22=2.353−0.0212vwr+1.074g+1.163vwrg−0.0319g2,C23=0.2866+0.4467vwr+1.049g+0.038vwrg−0.0177g2+0.005g3+0.1479g2vwr,C33=0.7375−1.177vwr+0.45061g+1.076vwr1g−0.0231g2,where g is the ratio of the long- and short-elliptical axes, g = a/b, and vwr is the equivalent Poisson’s ratio.The calculated creep force and creep torque between the wheel and rail needs to be converted to the absolute coordinate system. The load array of the wheel-rail creep force is expressed as(40)FcwrLi=FcLxi,FcLyi,FcLziT=CwrLi⋅FxLi,FyLi,0T,FcwrRi=FcRxi,FcRyi,FcRziT=CwrRi⋅FxRi,FyRi,0T,McwrLi=McLxi,McLyi,McLziT=CwrLi⋅0,0,MzLiT,McwrRi=McRxi,McRyi,McRziT=CwrRi⋅0,0,MzRiT,where(41)CwrLi=cosφwi−cosδcL+ψwisinφwisinδcL+ψwisinφwisinφwicosδcL+ψwicosφwi−sinδcL+ψwicosφwi0sinδcL+ψwicosδcL+ψwi,CwrRi=cosφwi−cosδcR−ψwisinφwi−sinδcR−ψwisinφwisinφwicosδcR−ψwicosφwisinδcR−ψwicosφwi0−sinδcR−ψwicosδcR−ψwi.(3) Geometric Model of Wheel-Rail Contact. Various parameters in the wheel-rail contact-force calculation, namely, the contact rail curvature radius ρr, the wheel tread curvature radius ρw, the wheel rolling radius rw, the wheel-rail contact angle δc, the wheel-rail contact point position C, and the relative compression of the wheel-rail contact εz must be solved by the geometric model of wheel-rail contact. The wheel-rail contact geometry is a spatially dynamic problem that accounts for the freedoms of the wheelset bounce ψw, wheelset yaw φw, and rail bounce ψr. This model discretizes the rail and wheel tread firstly and solves the tread curvature and tangent angle at the discrete points. The contact point of the wheel and rail is considered to lie on a spatial trace line, which is calculated as follows [3]:(42)xwL,R′=xwlx+lxRwtanδL,R,ywL,R′=xwly∓Rw2−lx2lx2lytanδL,R±lzm+Yw,zwL,R′=xwlz∓Rw2−lx2lx2lztanδL,R∓lzm+Zw,where m=1−lx21+tan2δL,R, and lx, ly, and lz are the X, Y, and Z direction cosines, respectively: lx=−cosψwsinφw, ly=cosψwsinφw, and lz=sinψw·dw is the rolling circular abscissa of the wheel tread.The random track random is regarded as the geometric displacement of the rail. Considering the transverse, vertical, and torsional degrees of freedom of the rail and the track irregularity, the discrete coordinate transformation of the rail is expressed as follows:(43)yrL,R′=yrL,RcosψrL±ψ0+zrL,RsinψrL±ψ0+YrL,R+ηyR,L,zrL,R′=−yrL,RsinψrL±ψ0+zrL,RcosψrL±ψ0+ZrL,R+ηzR,L,where, ηyR,L and ηzR,L are the vertical profiles and alignment irregularity of the left and right rail, respectively.As shown in Figure10, the three-dimensional curve of the wheel tread (xrL,R′, yrL,R′, zrL,R′) is projected onto the Y-O-Z plane and shifted upward by a fixed distance εz0. Interpolating between the discrete coordinates of the wheel tread (yzrL,R′, zyrL,R′) and the discrete coordinates of the rail (yzrL,R′, zyrL,R′), the model finds the point C that minimizes the wheel-to-rail distance in the Y-O-Z plane. The contact rail curvature radius ρr, wheel tread curvature radius ρw, wheel rolling radius rw, and wheel-rail contact angle δc are then obtained by interpolating between the coordinates of the C points. The normal compression of the wheel and rail is given by(44)εnL,R=εzL,RcosδcL,R±ψw.Figure 10
Comparison of track irregularity samples. (a) Vertical profiles of track irregularities, (b) track alignment irregularities, (c) track cross-level irregularities, and (d) track gage irregularities.
(a)
(b)
(c)
(d)Finally, the normal force between the wheel and rail is obtained by inserting the wheel-rail normal compression into the wheel-rail normal contact model.


## 2.3. Numerical Simulation and Verification of the Random Track Irregularity Model

In general, when simulating the random track irregularity time-domain samples, they often need dozens or hundreds of random variables. If the number theory method is used to select the random variables, it will inevitably produce certain aggregations, especially when high-dimensional selection is performed. To accurately reflect the entire random variable space, it is necessary to distribute the random variable points in the random variable space as much as possible. The generation of the clustering phenomenon will increase the number of random samples required for the calculation to some extent. Therefore, the random track irregularities time-domain samples were simulated by the spectral representation random-function method, which reduces the required number of random variables in the simulation. In a one-dimensional, univariate, stochastic process of the power spectral density functionS (ω), η(t) is expressed as [29, 30].(45)ηt=∫0∞cosωtdUtω+sinωtdVtω,where Ut(ω) and Vt(ω) are the spectral processes of a random process X(t), and their increments dUt(ω) and dVt(ω) satisfy the basic conditions of the spectral representation of a stochastic process:(46)EdUtω=EdVtω=0,ω≥0,EdUt2ω=EdVt2ω=2SXt,ωdω,ω≥0,EdUtωdUtω′=EdVtωdVtω′=0,ω,ω′≥0,ω′≠ω,EdUtωdVtω′=0,ω,ω′≥0Discretizing the above formula and introducing a set of standard orthogonal random variables {Xk, Yk}, η(t) is approximately discretized as(47)ηkt≈∑k=0∞cosωktΔUωk+sinωktΔVωk≈∑k=1NωcosωktσkXk+sinωktσkYk,where Δω is assumed sufficiently small, and dUt(ω) ≈ ΔUt(ω) = σkXk, dVt(ω) ≈ ΔVt(ω) = σkYk, and σk = [2SX(ωk) · Δω]1/2. It is easily proved that ΔUt(ω) and ΔVt(ω) satisfy the above basic conditions. Nω is the number of discrete frequencies. When the frequency components are uniformly dispersed, then ωk = ωmin + k · Δω, with Δω = (ωu − ωl)/N · ωu and ωl as the upper and lower limits, respectively, of the discrete frequency points.The orthogonal random variables {Xk, Yk} are obtained from the standard orthogonal variables {X¯n,Y¯n} by random mapping (as described in [30]). The standard orthogonal random variables {X¯n,Y¯n} can be constructed from orthogonal Legendre polynomials or Hartley orthogonal basis functions. The standard orthogonal variables {X¯n,Y¯n} are represented as(48)X¯n=casnΘ1,Y¯n=casnΘ2,n=1,2,…,N,where cas(x) denotes the Hartley orthogonal basis function, cas(x) = cos(x) + sin(x). Θ1 and Θ2 are independent random variables uniformly distributed in the interval [0, 2π].As an example, we consider the “TB/T3352-2014 high-speed railway ballastless track irregularity spectrum” implemented in 2015. Three hundred samples of random track irregularities were generated by the spectral representation random-function method. When generating the power spectrum of the track irregularities, the unit length and sampling interval were set to 1024 m and 0.25 m, respectively, with 4096 points in a single unit [31]. Therefore, when simulating the samples in the time domain, the discretized number of frequencies Nω was 4096. The vertical profile, alignment, cross-level, and gage irregularity of the track were, respectively, expressed as follows:Track vertical profile irregularity:(49)ηVit=2SVωkΔω⋅∑k=1NωcosωktX1k+sinωktX2k.Alignment irregularity:(50)ηAit=2SAωkΔω⋅∑k=1NωcosωktX3k+sinωktX4k.Cross-level irregularity:(51)ηCit=2SCωkΔω⋅∑k=1NωcosωktX5k+sinωktX6k.Gage irregularity:(52)ηGit=2SGωkΔω⋅∑k=1NωcosωktX7k+sinωktX8k.It should be noted that the track irregularity spectrum is averaged over many measured data. As the random error between a single sample and the track irregularity spectrum is 100%, the accuracy of the simulation method cannot be validated by comparing the simulated power spectrum of the track irregularity with that obtained from a single or finite-length sample, but validated by comparing the simulated power spectrum of the track irregularity with that obtained from a sample set.Figure10 compares the means and standard deviations of the track irregularity in the time domain and the track irregularity spectrum in [31]. The power spectrum of each frequency point in the sample set was tested by χ2 distribution with a significance level of α = 0.05 and two degrees of freedom with the Kolmogorov–Smirnov method, verifying the accuracy of the random track irregularities simulated by the present method.The power spectral density and track irregularity power spectrum derived from a single sample did not accurately estimate the vertical profile, alignment, cross-level, or gage irregularity of the track (Figure10). However, the mean power spectrum calculated from the sample set of track irregularities favorably agreed with the target spectrum. The maximum difference between the mean vertical profile, alignment, cross-level, and gage irregularity and those obtained from the target spectrum were 0.18 mm, 0.10 mm, 0.07 mm, and 0.08 mm, respectively, with maximum standard deviations of 4.68%, 5.77%, 6.97%, and 6.99%, respectively. In contrast, the maximum deviations in the power spectrum of the sample set were 14.54%, 10.23%, 10.07%, and 9.87%, respectively. In the power spectral densities at different frequencies of the sample set, the vertical profile, alignment, cross-level, and gage irregularity of the track all satisfied the 2-degree-of-freedom χ2 distribution (with 98.00%, 96.33%, 98.00%, and 98.33% of the frequency points satisfying the distribution, respectively). Each index of the sample set well agreed with the target spectrum, and only eight random variables were required to extract the indices from the time-domain samples. The method also retained a large number of frequency components. The commonly used trigonometric method, when simulating track random irregularity time-domain samples, requires enough random phase. When there are only 8 random phases, the simulated track random irregularity time-domain samples will have periodic repeatability and cannot be used in the calculation of the vehicle-rail coupled dynamic system. It can be seen that the numerical simulation method in this paper can greatly reduce the number of random variables required to simulate track random irregularity time-domain samples.


## 2.4. Numerical Solution Method

The track-vehicle coupled dynamic model was solved by an explicit-implicit integration method, with time integral step less than 0.0001 seconds, and the mode number of rail is taken as half the number of rail support points. The displacement and speed of the train and rail at the (n + 1)th moment were computed from the displacement, speed, and acceleration in the nth and (n − 1)th moments. The wheel-rail force was then calculated by the wheel-rail contact model, using the freedoms of the wheelset and rail at the (n + 1)th moment. After forming the wheel-rail load matrix, the fastener model was solved with the rail and slab freedoms, obtaining the wheel-rail load matrix. The vehicle and rail acceleration at the (n + 1)th moment was solved by the vehicle and rail motion equation, and the dynamic response of the slab and base at the (n + 1)th moment under the fastener force load was obtained by the Newmark-β integral method. Finally, the generalized probability density was evolved by the finite difference discretization of the modified Lax–Wendroff scheme [16].


## 3. Numerical Calculation and Analysis

### 3.1. Model Verification

#### 3.1.1. Calculation Conditions

The vehicle parameters are shown in Table2. The analysis was performed on the Chinese high-speed ballastless railway irregularity spectrum, with an irregularity wavelength range of 2–120 m. The vehicle speed was 350 km/h. The track irregularities in the time domain were sampled by the random irregularity numerical simulation method described in Section 2.3. The evaluation indices were the left/right wheel load reduction rates, the left/right wheel derailment coefficients, lateral force of the wheel axis, and the lateral/vertical accelerations of the car body. The calculation results of single and multiple samples were compared and verified.Table 2
Vehicle dynamics parameters.
Item
Notation
Unit
Value
Car body
Mass
m
c
kg
33.2 × 103
Mass moment of inertia aboutx-axis
Ix
c
kg·m2
1.07568 × 105
Mass moment of inertia abouty-axis
Iy
c
Kg·m2
1.6268 × 106
Mass moment of inertia aboutz-axis
Iz
c
kg·m2
1.4027 × 106
Bogie
Mass of bogie
m
t
kg
2.6 × 103
Mass moment of inertia aboutx-axis
Ix
t
kg·m2
2106
Mass moment of inertia abouty-axis
Iy
t
kg·m2
1424
Mass moment of inertia aboutz-axis
Iz
t
kg·m2
2600
Wheelset
Mass of wheelset
m
w
kg
1.97 × 103
Mass moment of inertia aboutx-axis
I
x
w
kg·m2
623
Mass moment of inertia abouty-axis
I
y
w
kg·m2
78
Mass moment of inertia aboutz-axis
I
z
w
kg·m2
623
Primary suspension system
Longitudinal spring stiffness
k
px
N/m
9.8 × 106
Lateral spring stiffness
k
py
N/m
9.8 × 106
Vertical spring stiffness
k
pz
N/m
1.176 × 106
Longitudinal damping coefficient
c
px
N·s/m
0
Lateral damping coefficient
c
py
N·s/m
0
Vertical damping coefficient
c
pz
N·s/m
0.196 × 103
Second suspension system
Longitudinal spring stiffness
k
sx
N/m
1.784 × 105
Lateral spring stiffness
k
sy
N/m
1.931 × 105
Vertical spring stiffness
k
sz
N/m
1.931 × 105
Longitudinal damping coefficient
c
sx
N·s/m
0
Lateral damping coefficient
c
sy
N·s/m
5.88 × 104
Vertical damping coefficient
c
sz
N·s/m
9.8 × 104
Half the transverse distance between vertical primary suspension systems
d
wk
m
1.0
Half the transverse distance between vertical primary suspension system absorber
d
wc
m
1.0
Half the transverse distance between vertical second suspension systems
d
sk
m
1.25
Half the transverse distance between vertical second suspension systems absorber
d
sc
m
1.25
Vertical distance between center of gravity of car body and lateral secondary suspension system
H
cb
m
0.644
Vertical distance between lateral secondary suspension system and center of gravity of bogie
H
bt
m
0.29
Vertical distance between center of gravity of bogie and lateral primary suspension system
H
tw
m
0.08
Half the longitudinal distance between center of rear bogie and center of front bogie
L
c
m
8.75
Half the bogie axle base
L
t
m
1.25


#### 3.1.2. Comparison and Verification of the Single-Sample Calculations

The vehicle-track coupled dynamic model proposed by Zhai and Chen [2, 3] has been widely used, so the dynamic model is verified by comparing the model in this paper with Zhai-model. Figure 11 compares the present calculation results with those of the Zhai-model in the case of a single time-domain sample from the track irregularity data. The calculation parameters and random track irregularity sample were the same in both methods. As the forces on the left and right wheel-rail are unequal, the wheel load reduction rates and derailment coefficients differed between the left and right wheel-rail contacts; see of Figures 11(a) and 11(b) for comparisons of the left and right wheel-rail states, respectively. The maximum differences between the left and right wheel-rail load reduction rates and derailment coefficients were 0.35 and 0.19, respectively. Therefore, the left or right wheel-rail state alone cannot fully reflect the wheel-rail state; both substates must be considered together. Comparing the present calculation model with the Zhai-model, the wheel load reduction rates and derailment coefficients differed by (at most) 8.74% and 7.62%, respectively. Meanwhile, the maximum differences in the vertical and transverse accelerations of the car body were 7.21% and 8.99%, respectively, and the maximum difference in the transverse force of the wheel axle was 6.78%. In summary, the results of the vehicle-track coupled dynamics model and the classical Zhai-model differed by less than 10%, confirming the accuracy of the model developed here.Figure 11
Comparison of calculation results for a single sample. (a) Reduction rate of wheel load, (b) derailment coefficient of wheel load, (c) acceleration of car body, and (d) lateral wheelset force.
(a)
(b)
(c)
(d)


#### 3.1.3. Comparison and Verification of Multisample Calculation Results

Random variables were selected by the MCM, and the random track irregularity was solved by the trigonometric series method. Five thousand random track irregularity samples were selected for computing the dynamic response of the vehicle-track coupled power system, and the statistical characteristics of each evaluation index were calculated. The probabilistic evolution processes of the indices were then calculated by the developed model, varying the number of samples as 100, 500, and 1000. Figure12 plots the cumulative probabilities of each index for the different numbers of samples and (for comparison) the cumulative probability in the MCM method with 5000 samples. The sectional cumulative probabilities of the left/right wheel load reduction rates, the left/right wheel derailment coefficients, the vertical/transverse accelerations of the car body, and the lateral force of the wheel axle well agreed with the statistical results of many MCM calculations. However, the developed model greatly reduced the number of required irregularity samples. Figures 12(a) and 12(c) show the sectional cumulative probabilities of the wheel load reduction rates and derailment coefficients, respectively, in the four cases at a fixed distance on the left side (the equivalent plots on the right side are shown in Figures 12(b) and 12(d), respectively). In the single sample analysis, the load reduction rates and derailment coefficients largely differed between the left and right sides, but after selecting 100, 500, or 1000 samples, the cumulative probability distributions of the left and right wheel-rail indices were consistent. The wheel load reduction rates of the developed model with 100, 500, and 1000 samples differed by (at most) 0.068, 0.049, and 0.033, respectively, from those of the MCM with 5000 samples. The corresponding differences between the derailment coefficients computed by the developed method and the MCM were 0.075, 0.040, and 0.029, respectively. Clearly, the developed computational model with multisampling achieves the same computational accuracy as the largely sampled MCM.Figure 12
Comparison of cumulative probabilities of the major vehicle-safety evaluation indices. Reduction rate of left wheel load (a) and right wheel load (b). Derailment coefficient of left wheel load (c) and right wheel load (d). Vertical acceleration of car body (e). Lateral acceleration of car body (f). Lateral wheelset force (g).
(a)
(b)
(c)
(d)
(e)
(f)
(g)


### 3.2. Reliability Analysis of Vehicle Operation Safety

The ballastless track structure in China offers a comfortable riding experience. When a single vehicle is traveling at 350 km/h, its safety evaluation index (calculated by the developed method) meets the requirements of common specifications. As limiting values in the safety evaluation, the wheel load reduction rate was taken as 0.3, the derailment coefficient as 0.15, the lateral and vertical accelerations of the car body as 1.0 m/s2 and 0.5 m/s2, respectively, and the lateral force as 13 kN. Figure 13 plots the equal-probability and reliability curves of the safety evaluation indices of the running vehicle. The developed model easily evolved the probability of each evaluation index, and all indices spatiotemporally varied with the running mileage. Although the running mileage influenced the reliability of the different indices, it did not affect their equal-probability evolution curves, because the random track irregularity time-domain samples obtained from the track irregularity spectrum is a stationary random process. Table 3 lists the maximum and minimum probabilities of each safety evaluation index under the calculation conditions. The probability values under different limits were quickly calculated, confirming the convenience of the developed model in train operation safety evaluation.Figure 13
Equal-probability and reliability curves of the train safety evaluation indices. Reduction rates of left wheel load (a) and right wheel load (b). Derailment coefficients of left wheel load (c) and right wheel load (d). Vertical acceleration of car body (e). Lateral acceleration of car body (f). Lateral wheelset force (g).
(a)
(b)
(c)
(d)
(e)
(f)
(g)Table 3
Probabilities of limit values of the vehicle-safety evaluation indices.
Evaluating indicator
Reduction rate
Derailment coefficient
Vertical acceleration of car body
Lateral acceleration of car body
Lateral wheelset force
Left
Right
Left
Right
Min
0.9787
0.9857
0.9874
0.9821
0.9885
0.9853
0.9725
Max
0.9989
0.9988
0.9989
0.9984
0.9997
0.9987
0.9948


### 3.3. Train Safety Evaluation Indices Limit Value with Different Probabilities

Given a vehicle and track structure with random track irregularities, the calculation results of a single sample cannot obtain the limits of the evaluation indices. In the presence of random irregularities, the limit values can be obtained only with a certain probability. For this purpose, the probabilistic index was introduced, providing a new way of computing the limit values of the indices of vehicle-track coupled dynamics.Figure14 plots the evaluation indices variation curve with vehicle speed computed by the developed model with different index probabilities (1, 2, and 3 times the standard deviation of a normal distribution). At 68.3% probability, the evaluation indices were relatively insensitive to vehicle speed, but at 99.7% probability, they significantly increased with speed. The derailment coefficient was barely affected by vehicle speeds below 400 km/h but rapidly increased at higher speeds (at 450 km/h, the derailment coefficient value at 99.7% probability increased by 123.8%).Figure 14
Speed dependences of the probabilistic train safety evaluation indices. Reduction rates of left wheel load (a) and right wheel load (b). Derailment coefficients of left wheel load (c) and right wheel load (d). Vertical acceleration of car body (e). Lateral acceleration of car body (f). Lateral wheelset force (g).
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 3.1. Model Verification

### 3.1.1. Calculation Conditions

The vehicle parameters are shown in Table2. The analysis was performed on the Chinese high-speed ballastless railway irregularity spectrum, with an irregularity wavelength range of 2–120 m. The vehicle speed was 350 km/h. The track irregularities in the time domain were sampled by the random irregularity numerical simulation method described in Section 2.3. The evaluation indices were the left/right wheel load reduction rates, the left/right wheel derailment coefficients, lateral force of the wheel axis, and the lateral/vertical accelerations of the car body. The calculation results of single and multiple samples were compared and verified.Table 2
Vehicle dynamics parameters.
Item
Notation
Unit
Value
Car body
Mass
m
c
kg
33.2 × 103
Mass moment of inertia aboutx-axis
Ix
c
kg·m2
1.07568 × 105
Mass moment of inertia abouty-axis
Iy
c
Kg·m2
1.6268 × 106
Mass moment of inertia aboutz-axis
Iz
c
kg·m2
1.4027 × 106
Bogie
Mass of bogie
m
t
kg
2.6 × 103
Mass moment of inertia aboutx-axis
Ix
t
kg·m2
2106
Mass moment of inertia abouty-axis
Iy
t
kg·m2
1424
Mass moment of inertia aboutz-axis
Iz
t
kg·m2
2600
Wheelset
Mass of wheelset
m
w
kg
1.97 × 103
Mass moment of inertia aboutx-axis
I
x
w
kg·m2
623
Mass moment of inertia abouty-axis
I
y
w
kg·m2
78
Mass moment of inertia aboutz-axis
I
z
w
kg·m2
623
Primary suspension system
Longitudinal spring stiffness
k
px
N/m
9.8 × 106
Lateral spring stiffness
k
py
N/m
9.8 × 106
Vertical spring stiffness
k
pz
N/m
1.176 × 106
Longitudinal damping coefficient
c
px
N·s/m
0
Lateral damping coefficient
c
py
N·s/m
0
Vertical damping coefficient
c
pz
N·s/m
0.196 × 103
Second suspension system
Longitudinal spring stiffness
k
sx
N/m
1.784 × 105
Lateral spring stiffness
k
sy
N/m
1.931 × 105
Vertical spring stiffness
k
sz
N/m
1.931 × 105
Longitudinal damping coefficient
c
sx
N·s/m
0
Lateral damping coefficient
c
sy
N·s/m
5.88 × 104
Vertical damping coefficient
c
sz
N·s/m
9.8 × 104
Half the transverse distance between vertical primary suspension systems
d
wk
m
1.0
Half the transverse distance between vertical primary suspension system absorber
d
wc
m
1.0
Half the transverse distance between vertical second suspension systems
d
sk
m
1.25
Half the transverse distance between vertical second suspension systems absorber
d
sc
m
1.25
Vertical distance between center of gravity of car body and lateral secondary suspension system
H
cb
m
0.644
Vertical distance between lateral secondary suspension system and center of gravity of bogie
H
bt
m
0.29
Vertical distance between center of gravity of bogie and lateral primary suspension system
H
tw
m
0.08
Half the longitudinal distance between center of rear bogie and center of front bogie
L
c
m
8.75
Half the bogie axle base
L
t
m
1.25


### 3.1.2. Comparison and Verification of the Single-Sample Calculations

The vehicle-track coupled dynamic model proposed by Zhai and Chen [2, 3] has been widely used, so the dynamic model is verified by comparing the model in this paper with Zhai-model. Figure 11 compares the present calculation results with those of the Zhai-model in the case of a single time-domain sample from the track irregularity data. The calculation parameters and random track irregularity sample were the same in both methods. As the forces on the left and right wheel-rail are unequal, the wheel load reduction rates and derailment coefficients differed between the left and right wheel-rail contacts; see of Figures 11(a) and 11(b) for comparisons of the left and right wheel-rail states, respectively. The maximum differences between the left and right wheel-rail load reduction rates and derailment coefficients were 0.35 and 0.19, respectively. Therefore, the left or right wheel-rail state alone cannot fully reflect the wheel-rail state; both substates must be considered together. Comparing the present calculation model with the Zhai-model, the wheel load reduction rates and derailment coefficients differed by (at most) 8.74% and 7.62%, respectively. Meanwhile, the maximum differences in the vertical and transverse accelerations of the car body were 7.21% and 8.99%, respectively, and the maximum difference in the transverse force of the wheel axle was 6.78%. In summary, the results of the vehicle-track coupled dynamics model and the classical Zhai-model differed by less than 10%, confirming the accuracy of the model developed here.Figure 11
Comparison of calculation results for a single sample. (a) Reduction rate of wheel load, (b) derailment coefficient of wheel load, (c) acceleration of car body, and (d) lateral wheelset force.
(a)
(b)
(c)
(d)


### 3.1.3. Comparison and Verification of Multisample Calculation Results

Random variables were selected by the MCM, and the random track irregularity was solved by the trigonometric series method. Five thousand random track irregularity samples were selected for computing the dynamic response of the vehicle-track coupled power system, and the statistical characteristics of each evaluation index were calculated. The probabilistic evolution processes of the indices were then calculated by the developed model, varying the number of samples as 100, 500, and 1000. Figure12 plots the cumulative probabilities of each index for the different numbers of samples and (for comparison) the cumulative probability in the MCM method with 5000 samples. The sectional cumulative probabilities of the left/right wheel load reduction rates, the left/right wheel derailment coefficients, the vertical/transverse accelerations of the car body, and the lateral force of the wheel axle well agreed with the statistical results of many MCM calculations. However, the developed model greatly reduced the number of required irregularity samples. Figures 12(a) and 12(c) show the sectional cumulative probabilities of the wheel load reduction rates and derailment coefficients, respectively, in the four cases at a fixed distance on the left side (the equivalent plots on the right side are shown in Figures 12(b) and 12(d), respectively). In the single sample analysis, the load reduction rates and derailment coefficients largely differed between the left and right sides, but after selecting 100, 500, or 1000 samples, the cumulative probability distributions of the left and right wheel-rail indices were consistent. The wheel load reduction rates of the developed model with 100, 500, and 1000 samples differed by (at most) 0.068, 0.049, and 0.033, respectively, from those of the MCM with 5000 samples. The corresponding differences between the derailment coefficients computed by the developed method and the MCM were 0.075, 0.040, and 0.029, respectively. Clearly, the developed computational model with multisampling achieves the same computational accuracy as the largely sampled MCM.Figure 12
Comparison of cumulative probabilities of the major vehicle-safety evaluation indices. Reduction rate of left wheel load (a) and right wheel load (b). Derailment coefficient of left wheel load (c) and right wheel load (d). Vertical acceleration of car body (e). Lateral acceleration of car body (f). Lateral wheelset force (g).
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 3.1.1. Calculation Conditions

The vehicle parameters are shown in Table2. The analysis was performed on the Chinese high-speed ballastless railway irregularity spectrum, with an irregularity wavelength range of 2–120 m. The vehicle speed was 350 km/h. The track irregularities in the time domain were sampled by the random irregularity numerical simulation method described in Section 2.3. The evaluation indices were the left/right wheel load reduction rates, the left/right wheel derailment coefficients, lateral force of the wheel axis, and the lateral/vertical accelerations of the car body. The calculation results of single and multiple samples were compared and verified.Table 2
Vehicle dynamics parameters.
Item
Notation
Unit
Value
Car body
Mass
m
c
kg
33.2 × 103
Mass moment of inertia aboutx-axis
Ix
c
kg·m2
1.07568 × 105
Mass moment of inertia abouty-axis
Iy
c
Kg·m2
1.6268 × 106
Mass moment of inertia aboutz-axis
Iz
c
kg·m2
1.4027 × 106
Bogie
Mass of bogie
m
t
kg
2.6 × 103
Mass moment of inertia aboutx-axis
Ix
t
kg·m2
2106
Mass moment of inertia abouty-axis
Iy
t
kg·m2
1424
Mass moment of inertia aboutz-axis
Iz
t
kg·m2
2600
Wheelset
Mass of wheelset
m
w
kg
1.97 × 103
Mass moment of inertia aboutx-axis
I
x
w
kg·m2
623
Mass moment of inertia abouty-axis
I
y
w
kg·m2
78
Mass moment of inertia aboutz-axis
I
z
w
kg·m2
623
Primary suspension system
Longitudinal spring stiffness
k
px
N/m
9.8 × 106
Lateral spring stiffness
k
py
N/m
9.8 × 106
Vertical spring stiffness
k
pz
N/m
1.176 × 106
Longitudinal damping coefficient
c
px
N·s/m
0
Lateral damping coefficient
c
py
N·s/m
0
Vertical damping coefficient
c
pz
N·s/m
0.196 × 103
Second suspension system
Longitudinal spring stiffness
k
sx
N/m
1.784 × 105
Lateral spring stiffness
k
sy
N/m
1.931 × 105
Vertical spring stiffness
k
sz
N/m
1.931 × 105
Longitudinal damping coefficient
c
sx
N·s/m
0
Lateral damping coefficient
c
sy
N·s/m
5.88 × 104
Vertical damping coefficient
c
sz
N·s/m
9.8 × 104
Half the transverse distance between vertical primary suspension systems
d
wk
m
1.0
Half the transverse distance between vertical primary suspension system absorber
d
wc
m
1.0
Half the transverse distance between vertical second suspension systems
d
sk
m
1.25
Half the transverse distance between vertical second suspension systems absorber
d
sc
m
1.25
Vertical distance between center of gravity of car body and lateral secondary suspension system
H
cb
m
0.644
Vertical distance between lateral secondary suspension system and center of gravity of bogie
H
bt
m
0.29
Vertical distance between center of gravity of bogie and lateral primary suspension system
H
tw
m
0.08
Half the longitudinal distance between center of rear bogie and center of front bogie
L
c
m
8.75
Half the bogie axle base
L
t
m
1.25


## 3.1.2. Comparison and Verification of the Single-Sample Calculations

The vehicle-track coupled dynamic model proposed by Zhai and Chen [2, 3] has been widely used, so the dynamic model is verified by comparing the model in this paper with Zhai-model. Figure 11 compares the present calculation results with those of the Zhai-model in the case of a single time-domain sample from the track irregularity data. The calculation parameters and random track irregularity sample were the same in both methods. As the forces on the left and right wheel-rail are unequal, the wheel load reduction rates and derailment coefficients differed between the left and right wheel-rail contacts; see of Figures 11(a) and 11(b) for comparisons of the left and right wheel-rail states, respectively. The maximum differences between the left and right wheel-rail load reduction rates and derailment coefficients were 0.35 and 0.19, respectively. Therefore, the left or right wheel-rail state alone cannot fully reflect the wheel-rail state; both substates must be considered together. Comparing the present calculation model with the Zhai-model, the wheel load reduction rates and derailment coefficients differed by (at most) 8.74% and 7.62%, respectively. Meanwhile, the maximum differences in the vertical and transverse accelerations of the car body were 7.21% and 8.99%, respectively, and the maximum difference in the transverse force of the wheel axle was 6.78%. In summary, the results of the vehicle-track coupled dynamics model and the classical Zhai-model differed by less than 10%, confirming the accuracy of the model developed here.Figure 11
Comparison of calculation results for a single sample. (a) Reduction rate of wheel load, (b) derailment coefficient of wheel load, (c) acceleration of car body, and (d) lateral wheelset force.
(a)
(b)
(c)
(d)


## 3.1.3. Comparison and Verification of Multisample Calculation Results

Random variables were selected by the MCM, and the random track irregularity was solved by the trigonometric series method. Five thousand random track irregularity samples were selected for computing the dynamic response of the vehicle-track coupled power system, and the statistical characteristics of each evaluation index were calculated. The probabilistic evolution processes of the indices were then calculated by the developed model, varying the number of samples as 100, 500, and 1000. Figure12 plots the cumulative probabilities of each index for the different numbers of samples and (for comparison) the cumulative probability in the MCM method with 5000 samples. The sectional cumulative probabilities of the left/right wheel load reduction rates, the left/right wheel derailment coefficients, the vertical/transverse accelerations of the car body, and the lateral force of the wheel axle well agreed with the statistical results of many MCM calculations. However, the developed model greatly reduced the number of required irregularity samples. Figures 12(a) and 12(c) show the sectional cumulative probabilities of the wheel load reduction rates and derailment coefficients, respectively, in the four cases at a fixed distance on the left side (the equivalent plots on the right side are shown in Figures 12(b) and 12(d), respectively). In the single sample analysis, the load reduction rates and derailment coefficients largely differed between the left and right sides, but after selecting 100, 500, or 1000 samples, the cumulative probability distributions of the left and right wheel-rail indices were consistent. The wheel load reduction rates of the developed model with 100, 500, and 1000 samples differed by (at most) 0.068, 0.049, and 0.033, respectively, from those of the MCM with 5000 samples. The corresponding differences between the derailment coefficients computed by the developed method and the MCM were 0.075, 0.040, and 0.029, respectively. Clearly, the developed computational model with multisampling achieves the same computational accuracy as the largely sampled MCM.Figure 12
Comparison of cumulative probabilities of the major vehicle-safety evaluation indices. Reduction rate of left wheel load (a) and right wheel load (b). Derailment coefficient of left wheel load (c) and right wheel load (d). Vertical acceleration of car body (e). Lateral acceleration of car body (f). Lateral wheelset force (g).
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 3.2. Reliability Analysis of Vehicle Operation Safety

The ballastless track structure in China offers a comfortable riding experience. When a single vehicle is traveling at 350 km/h, its safety evaluation index (calculated by the developed method) meets the requirements of common specifications. As limiting values in the safety evaluation, the wheel load reduction rate was taken as 0.3, the derailment coefficient as 0.15, the lateral and vertical accelerations of the car body as 1.0 m/s2 and 0.5 m/s2, respectively, and the lateral force as 13 kN. Figure 13 plots the equal-probability and reliability curves of the safety evaluation indices of the running vehicle. The developed model easily evolved the probability of each evaluation index, and all indices spatiotemporally varied with the running mileage. Although the running mileage influenced the reliability of the different indices, it did not affect their equal-probability evolution curves, because the random track irregularity time-domain samples obtained from the track irregularity spectrum is a stationary random process. Table 3 lists the maximum and minimum probabilities of each safety evaluation index under the calculation conditions. The probability values under different limits were quickly calculated, confirming the convenience of the developed model in train operation safety evaluation.Figure 13
Equal-probability and reliability curves of the train safety evaluation indices. Reduction rates of left wheel load (a) and right wheel load (b). Derailment coefficients of left wheel load (c) and right wheel load (d). Vertical acceleration of car body (e). Lateral acceleration of car body (f). Lateral wheelset force (g).
(a)
(b)
(c)
(d)
(e)
(f)
(g)Table 3
Probabilities of limit values of the vehicle-safety evaluation indices.
Evaluating indicator
Reduction rate
Derailment coefficient
Vertical acceleration of car body
Lateral acceleration of car body
Lateral wheelset force
Left
Right
Left
Right
Min
0.9787
0.9857
0.9874
0.9821
0.9885
0.9853
0.9725
Max
0.9989
0.9988
0.9989
0.9984
0.9997
0.9987
0.9948


## 3.3. Train Safety Evaluation Indices Limit Value with Different Probabilities

Given a vehicle and track structure with random track irregularities, the calculation results of a single sample cannot obtain the limits of the evaluation indices. In the presence of random irregularities, the limit values can be obtained only with a certain probability. For this purpose, the probabilistic index was introduced, providing a new way of computing the limit values of the indices of vehicle-track coupled dynamics.Figure14 plots the evaluation indices variation curve with vehicle speed computed by the developed model with different index probabilities (1, 2, and 3 times the standard deviation of a normal distribution). At 68.3% probability, the evaluation indices were relatively insensitive to vehicle speed, but at 99.7% probability, they significantly increased with speed. The derailment coefficient was barely affected by vehicle speeds below 400 km/h but rapidly increased at higher speeds (at 450 km/h, the derailment coefficient value at 99.7% probability increased by 123.8%).Figure 14
Speed dependences of the probabilistic train safety evaluation indices. Reduction rates of left wheel load (a) and right wheel load (b). Derailment coefficients of left wheel load (c) and right wheel load (d). Vertical acceleration of car body (e). Lateral acceleration of car body (f). Lateral wheelset force (g).
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 4. Conclusions

This paper presents an evaluation model of safe railway operation based on probability theory and a nonlinear vehicle-track coupled power system. The model more accurately captures the force state between the wheel and rail than the linearized vehicle-track coupling dynamic system. The random track irregularities were sampled in the time domain by a stochastic function-spectral representation method, which reduces the required number of samplings. The developed simulation method was verified in comparisons with an established method. The safety of the running vehicles was probabilistically determined by the evaluation indices (derailment coefficient, acceleration of the car body, and the lateral force of the wheel axle). The model results were compared with those of the Zhai-model and the MCM. Next, the vehicle-safety evaluation indices were computed for different probabilities of their limit values. At the lower probability (68.3%), the limiting evaluation indices were little affected by the vehicle speed, but at the highest probability (99.7%), they were significantly increased at higher vehicle speeds. The derailment coefficient deserves special attention when driving at high speed. The developed computational model provides a new research direction for vehicle operation safety evaluation. In future work, the model will be applied to safe railway operation under different conditions, such as different structures, different track or foundation diseases, and different track irregularities.



---

*Source: 3236093-2019-07-02.xml*