Research Article  Open Access
FuShan Liu, ZhiPing Zeng, WeiDong Wang, Abdulmumin Ahmed Shuaibu, "Stochastic Analysis of Nonlinear VehicleTrack Coupled Dynamic System and Its Application in Vehicle Operation Safety Evaluation", Shock and Vibration, vol. 2019, Article ID 3236093, 23 pages, 2019. https://doi.org/10.1155/2019/3236093
Stochastic Analysis of Nonlinear VehicleTrack Coupled Dynamic System and Its Application in Vehicle Operation Safety Evaluation
Abstract
This paper presents a vehicle operation safety evaluation model; to this end, a nonlinear vehicletrack coupled dynamic system stochastic analysis model under random irregularity excitations based on probability density evolution method was developed. The nonlinear coupled vehicletrack dynamic system is used to accurately describe the wheelrail contact state. The stochastic functionspectral 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 Zhaimodel 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 wheelrail derailment coefficient increases rapidly at running speeds above 400 km/h. The computational model provides a novel direction for vehicle operation safety evaluation.
1. Introduction
Vehicletrack 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 highspeed 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 traintrack coupled dynamics models or trainbridge coupled dynamics models have been established to analyze the dynamic response of the system. Wu et al. [1] developed a vehiclerailbridge 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 threedimensional space, and at the same time, a threedimensional space wheelrail contact relationship has also been developed. Lou and Zeng [4] have derived the equations of motion in matrix form for the railway freight vehiclebridge interaction system based on energy approach; this modeling method also has been widely used. Dinh et al. [5] have developed a formulation of threedimensional dynamic interactions between a bridge and a highspeed train using wheelrail 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 vehicletrack or vehicletrackbridge coupling dynamics models. The development of vehicletrack 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 vehicletrack coupled dynamic system and have obtained many meaningful research results.
Due to the complexity of vehicletrack 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 vehicletrack 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 vehicletrackcoupled power systems. Lu et al. [11], Zhang et al. [12], and Zeng et al. [13] established random vibration analysis model for linear vehicletrackbridge 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 threedimensional trainbridge systems subjected to lateral horizontal earthquakes [14]. Yu et al. [15] studied the random vibrations of a linearized vehiclerailbridge 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 spacetime analysis method for vehicletrackcoupled 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 vehicletrackcoupled 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 vehicletrack coupled power system or vehiclebridge coupled dynamic system. Rocha et al. presented a probabilistic methodology for the safety assessment of shortspan railway bridges for highspeed 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 shortspan railway bridges is compared [21]. Using the principle of virtual work, utilizing a linearized wheelrail 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 vehiclebridge 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 highspeed railway bridge system.
In the existing random vibration analysis of vehicletrack 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 wheelrail relationship is often reduced to a linear model. The wheelrail 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 wheelrail relationship; this makes the wheelrail 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 wheelrail force, and it will also lead to a difference of the input track irregularity excitation input. Other than this, models with nonlinear wheelrail contact relationship allow for the evaluation of vehicle safety in a more realistic way than linear methods [24]. Numerous studies on vehiclebridge dynamics have indicated that nonlinear interaction has profound effects on the system responses [25]. Therefore, the nonlinear vehicletrack coupled dynamics model is more helpful for vehiclesafety evaluation, and the calculation results are more accurate.
In this paper, a stochastic analysis model of nonlinear vehicletrack coupled dynamic system under random irregularity excitations was established based on probability density evolution method. The nonlinear wheelrail 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, functionspectral 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 highdimensional random variables. The probability density evolution of the left/right wheelrail load reduction, the left/right wheelrail 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 VehicleTrack Coupled Dynamic System
Figure 1 shows the calculation flow of the nonlinear vehicletrack 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 N_{pt} random track irregularities timedomain samples, the dynamic response calculation of the coupled vehicletrack system under random excitations, calculation of the probability density evolution of the safety evaluation index, and the train safety evaluation.
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 vehicletrack coupled dynamic system is expressed aswhere t denotes time, and , , 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 wheelrail 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 vehicletrack coupled dynamic system derives from the random track irregularities. A nonlinear vehicletrack dynamic system ξ (t) is a probabilistic conservative stochastic system that excludes other random sources. According to the probabilistic conservation principle of stochastic systems, probabilistic conservationaugmented systems are constructed as [ξ (t), Θ] withwhere P{·} denotes the probability of random events. Ω_{t} and are the time distribution regions of the system state at times t and t_{0}, respectively, and Ω_{θ} is the space distribution region of the system state of sample θ. The above expression is equivalent to
Expressing the joint probability density function (ξ(t), Θ) as , the above expression becomeswhere D/D_{t} is the derivative of matter [16].
The generalized probability density evolves by the following equation [16]:
In this formula, is the free velocity of ξ (t) under the random sample {Θ = θ}. Its initial probability is .
Solving the above evolution equation of generalized probability density under the initial conditions, we obtain the joint probability density function of sample θ, . Furthermore, the probability density function of P_{ξΘ} is obtained by integrating θ:
For the probabilistic conservative vehicletrack 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 vehicletrack coupled dynamic system equation (1) to obtain 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 vehicletrack 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 lowdeviation sequence selected by number theory, which reduces the required number of random samples and accelerates the convergence of the calculation. However, as the lowdeviation sequence is uniformly distributed, it cannot be selected from a highdimensional 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 functionspectral representation method, which requires just eight random variables. The generation method and validation of the sampling will be discussed in Section 2.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
The assigned probability P_{k} of random sample θ_{k} is given bywhere is the probability density of random variable Θ in Voronoi region V(θ_{k}). The joint probability density in the total partitioned space satisfies
This method obtains the initial probability of each random point sample.
2.2. VehicleTrack Coupled Dynamics Model
Figure 2 shows the vehicletrack coupling dynamic model. The vehicle is modeled by multirigid body dynamics, and the rail is regarded as a finitelength, elasticpoint supported Euler beam. The dynamic vibration equation of the beam is solved by the Ritz method. The wheelrail relationship is modeled by nonlinear Hertz contact theory and wheelrail 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 wheelrail relationship connects the vehicle and rail models, whereas the fastener force connects the rail model and the track slabfoundation structure.
2.2.1. Vehicle Dynamics Model
Figure 3 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).
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 bywhere [M_{V}], [C_{V}], [K_{V}], [G_{V}], [F_{V}] are the vehicle mass, damping, stiffness, additional force, and load matrices, respectively.
The freedom matrix of the vehicle is expressed aswhere the subscripts c, b, and 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 aswhere Y, Z, Ψ, φ, and β represent the freedoms in the longitudinal, transverse, bounce, roll, and yaw, respectively.
The vehiclemass matrix is given bywhere . Here m denotes the mass, and Ix, Iy, and Iz are the moments of inertia around the X, Y, and Z axes, respectively.
The vehiclestiffness matrix is given bywhere
The vehicledamping matrix is similar to the vehiclestiffness matrix. In fact, the vehicledamping 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 matrixwith
In the above expressions, G_{wi} is the matrix of moments acting on the ith wheelset caused by wheelset rolling rotations, G_{tj} is the matrix of moments acting on the jth bioge caused by absorber, G_{c} is the matrix of moments acting on the car body caused by absorber, and F_{xsLj} and F_{xsLj} are the additional forces imposed by the snakeshaped shock absorber. They are calculated as follows:
M_{ri} is the antiroll moment, given by
The vehicleload matrix is expressed aswhere A_{i} is the matrix of conversion coefficients between the wheelrail force and the vehicle freedoms. F_{wri} is the wheelrail force array, withwhere F_{L,Rxi} = F_{nL,Rxi} + F_{cL,Rxi}, F_{L,Ryi} = F_{nL,Ryi} + F_{cL,Ryi}, and F_{L,Rzi} = F_{nL,Rzi} + F_{cL,Rzi}. F_{n} and F_{c} 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 Zaxis 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 trackstructure model and rail forces are shown in Figures 4 and 5, respectively.
The following expressions are shown for the right rail (the left rail is treated analogously). The fastener force matrix is expressed as
The torque of the wheelrail force, which induces rail torsion, is
The differential equations of the vertical, transverse, and torsional vibrations of the rail are, respectively given in the following [3].
Vertical:
Transverse:
Torsional:
The track and base plates are considered as horizontally rigid bodies and are modeled as fournode bending plate elements. They are connected by a fastener element, which is simulated by a linear spring damping element. Slabbase 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 bywhere 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 slabbase is the fastening force. The matrix and modeling process of the track and base plates are detailed elsewhere [13].
2.2.3. WheelRail Contact Relation
The wheelrail relationship links the vehicle system to the track system. In the vehicletrack coupled dynamics, the wheelrail relationship mainly includes the wheelrail normalseeking solution, the wheelrail creep solution, and the wheelrail contact geometry. The wheelrail relationship is shown in Figure 6.
(1) The WheelRail 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 aswhere r_{p}, , and r_{r} 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.
The normal force on the wheel and rail is expressed as
And the maximum pressure between the wheel and rail is given by
In the above expressions, is the equivalent modulus of elasticity, given by , and , where a and b are the lengths of the long and short axes, respectively, of the contact ellipse (with b < a). F_{1}(e) and F_{2}(e) are correction factors that depend only on the ratio of relative curvature radii . The variations of these factors with increasing are plotted in Figure 8.
F _{1}(e) and F_{2}(e) can be directly computed by a formula or directly read from Figure 9. To improve the efficiency of the computation, the curvefitting parameters obtained from Figure 9 are listed in Table 1, F_{1}(e) and F_{2}(e) were then obtained aswith x = ln(B/A).

When 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 wheelrail contact. The normal force between the wheel and rail is then obtained by solving the wheelrail normal contact model and determining the direction and magnitude variations in the normal wheelrail force at different contact points. Converting to the absolute coordinate system, the normal loads on the wheel and rail array are, respectively, expressed as
Here, δ_{cLi} and δ_{cRi} are the leftsided and rightsided wheelrail contact angles, respectively, obtained from the wheelrail contact geometry (details will be given in Section (2)).
(2) WheelRail Tangential Contact Model. The wheelrail 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 bywhere α is the correction coefficient, given as with andF_{x} and F_{y} denote the creep forces in the X and Yaxis directions, respectively, and M_{z} is the rotational creep moment. In linear creep theory, these components are, respectively, expressed aswhere f_{11} and f_{22} are the longitudinal and transverse creep coefficients, respectively, and f_{23} and f_{33} are rotational creep coefficients:
In the above expressions, G_{wr} is the composite shear modulus of the wheel and rail materials, and C_{ij} are the Kalker coefficients [27] (after Kalker, who calculated them). Meymandet al. [28] have fitted the Kalker coefficients under different contact conditions. Here, the C_{ij}s are expressed aswhere is the ratio of the long and shortelliptical axes, = a/b, and 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 wheelrail creep force is expressed aswhere
(3) Geometric Model of WheelRail Contact. Various parameters in the wheelrail contactforce calculation, namely, the contact rail curvature radius ρ_{r}, the wheel tread curvature radius , the wheel rolling radius , the wheelrail contact angle δ_{c}, the wheelrail contact point position C, and the relative compression of the wheelrail contact ε_{z} must be solved by the geometric model of wheelrail contact. The wheelrail contact geometry is a spatially dynamic problem that accounts for the freedoms of the wheelset bounce , wheelset yaw , 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]:where , and l_{x}, l_{y}, and l_{z} are the X, Y, and Z direction cosines, respectively: , , and 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:where, η_{yR,L} and η_{zR,L} are the vertical profiles and alignment irregularity of the left and right rail, respectively.
As shown in Figure 10, the threedimensional curve of the wheel tread (, , ) is projected onto the YOZ plane and shifted upward by a fixed distance ε_{z0}. Interpolating between the discrete coordinates of the wheel tread (, ) and the discrete coordinates of the rail (, ), the model finds the point C that minimizes the wheeltorail distance in the YOZ plane. The contact rail curvature radius ρ_{r}, wheel tread curvature radius , wheel rolling radius , and wheelrail 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
(a)
(b)
(c)
(d)
Finally, the normal force between the wheel and rail is obtained by inserting the wheelrail normal compression into the wheelrail normal contact model.
2.3. Numerical Simulation and Verification of the Random Track Irregularity Model
In general, when simulating the random track irregularity timedomain 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 highdimensional 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 timedomain samples were simulated by the spectral representation randomfunction method, which reduces the required number of random variables in the simulation. In a onedimensional, univariate, stochastic process of the power spectral density function S (ω), η(t) is expressed as [29, 30].where U_{t}(ω) and V_{t}(ω) are the spectral processes of a random process X(t), and their increments dU_{t}(ω) and dV_{t}(ω) satisfy the basic conditions of the spectral representation of a stochastic process:
Discretizing the above formula and introducing a set of standard orthogonal random variables {X_{k}, Y_{k}}, η(t) is approximately discretized aswhere Δω is assumed sufficiently small, and dU_{t}(ω) ≈ ΔU_{t}(ω) = σ_{k}X_{k}, dV_{t}(ω) ≈ ΔV_{t}(ω) = σ_{k}Y_{k}, and σ_{k} = [2S_{X}(ω_{k}) · Δω]^{1/2}. It is easily proved that ΔU_{t}(ω) and ΔV_{t}(ω) 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 {X_{k}, Y_{k}} are obtained from the standard orthogonal variables {} by random mapping (as described in [30]). The standard orthogonal random variables {} can be constructed from orthogonal Legendre polynomials or Hartley orthogonal basis functions. The standard orthogonal variables {} are represented aswhere 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/T33522014 highspeed railway ballastless track irregularity spectrum” implemented in 2015. Three hundred samples of random track irregularities were generated by the spectral representation randomfunction 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, crosslevel, and gage irregularity of the track were, respectively, expressed as follows:
Track vertical profile irregularity:
Alignment irregularity:
Crosslevel irregularity:
Gage irregularity:
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 finitelength sample, but validated by comparing the simulated power spectrum of the track irregularity with that obtained from a sample set.
Figure 10 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, crosslevel, or gage irregularity of the track (Figure 10). 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, crosslevel, 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, crosslevel, and gage irregularity of the track all satisfied the 2degreeoffreedom χ^{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 timedomain samples. The method also retained a large number of frequency components. The commonly used trigonometric method, when simulating track random irregularity timedomain samples, requires enough random phase. When there are only 8 random phases, the simulated track random irregularity timedomain samples will have periodic repeatability and cannot be used in the calculation of the vehiclerail 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 timedomain samples.
2.4. Numerical Solution Method
The trackvehicle coupled dynamic model was solved by an explicitimplicit 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 wheelrail force was then calculated by the wheelrail contact model, using the freedoms of the wheelset and rail at the (n + 1)th moment. After forming the wheelrail load matrix, the fastener model was solved with the rail and slab freedoms, obtaining the wheelrail 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 Table 2. The analysis was performed on the Chinese highspeed 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.

3.1.2. Comparison and Verification of the SingleSample Calculations
The vehicletrack 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 Zhaimodel. Figure 11 compares the present calculation results with those of the Zhaimodel in the case of a single timedomain 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 wheelrail are unequal, the wheel load reduction rates and derailment coefficients differed between the left and right wheelrail contacts; see of Figures 11(a) and 11(b) for comparisons of the left and right wheelrail states, respectively. The maximum differences between the left and right wheelrail load reduction rates and derailment coefficients were 0.35 and 0.19, respectively. Therefore, the left or right wheelrail state alone cannot fully reflect the wheelrail state; both substates must be considered together. Comparing the present calculation model with the Zhaimodel, 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 vehicletrack coupled dynamics model and the classical Zhaimodel differed by less than 10%, confirming the accuracy of the model developed here.
(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 vehicletrack 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. Figure 12 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 wheelrail 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.
(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/s^{2} and 0.5 m/s^{2}, respectively, and the lateral force as 13 kN. Figure 13 plots the equalprobability 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 equalprobability evolution curves, because the random track irregularity timedomain 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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)

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 vehicletrack coupled dynamics.
Figure 14 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%).
(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 vehicletrack coupled power system. The model more accurately captures the force state between the wheel and rail than the linearized vehicletrack coupling dynamic system. The random track irregularities were sampled in the time domain by a stochastic functionspectral 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 Zhaimodel and the MCM. Next, the vehiclesafety 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.
Data Availability
The analysis result data used to support the findings of this study are included within the article. The calculation data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The research work described in this paper was supported by the Highspeed Railway Joint Fund of National Natural Science Foundation of China (grant U1734208), the Major Program of National Natural Science Foundation of China (grant 11790283), the Hunan Provincial Natural Science Foundation of China (grant 2019JJ40384), and the Hunan Provincial Innovation Foundation for Postgraduate (CX2017B058).
References
 Y.S. Wu, Y.B. Yang, and J.D. Yau, “Threedimensional analysis of trainrailbridge interaction problems,” Vehicle System Dynamics, vol. 36, no. 1, pp. 1–35, 2001. View at: Publisher Site  Google Scholar
 G. Chen and W. M. Zhai, “A new wheel/rail spatially dynamic coupling model and its verification,” Vehicle System Dynamics, vol. 41, no. 4, pp. 301–322, 2004. View at: Publisher Site  Google Scholar
 W. M Zhai, VehicleTrack Coupled Dynamics, Science Press, Beijing, China, 4th edition, 2015.
 P. Lou and Q.Y. Zeng, “Formulation of equations of motion for a simply supported bridge under a moving railway freight vehicle,” Shock and Vibration, vol. 14, no. 6, pp. 429–446, 2007. View at: Publisher Site  Google Scholar
 V. N. Dinh, K. D. Kim, and P. Warnitchai, “Dynamic analysis of threedimensional bridgehighspeed train interactions using a wheelrail contact model,” Engineering Structures, vol. 31, no. 12, pp. 3090–3106, 2009. View at: Publisher Site  Google Scholar
 N. Zhang, H. Xia, and W. W. Guo, “Vehicle–bridge interaction analysis under highspeed trains,” Journal of Sound and Vibration, vol. 309, no. 3–5, pp. 407–425, 2008. View at: Publisher Site  Google Scholar
 Y. B. Yang and J. D. Yau, “Vertical and pitching resonance of train cars moving over a series of simple beams,” Journal of Sound and Vibration, vol. 337, pp. 135–149, 2015. View at: Publisher Site  Google Scholar
 G. Chen, W. M. Zhai, and H. Zuo, “Analysis of the random vibration responses characteristics of the vehicle track coupling system,” Journal of Traffic and Transportation Engineering, vol. 1, no. 1, pp. 13–16, 2001. View at: Google Scholar
 H. Xia, G. De Roeck, H. R. Zhang, and N. Zhang, “Dynamic analysis of trainbridge system and its application in steel girder reinforcement,” Computers & Structures, vol. 79, no. 2021, pp. 1851–1860, 2001. View at: Publisher Site  Google Scholar
 F. T. K. Au, Y. S. Cheng, and Y. K. Cheung, “Effects of random road surface roughness and longterm deflection of prestressed concrete girder and cablestayed bridges on impact due to moving vehicles,” Computers & Structures, vol. 79, no. 8, pp. 853–872, 2001. View at: Publisher Site  Google Scholar
 F. Lu, J. H. Lin, D. Kennedy, and F. W. Williams, “An algorithm to study nonstationary random vibrations of vehicle–bridge systems,” Computers & Structures, vol. 87, no. 34, pp. 177–185, 2009. View at: Publisher Site  Google Scholar
 Z. Zhang, Y. Zhang, J. Lin, Y. Zhao, W. P. Howson, and F. W. Williams, “Random vibration of a train traversing a bridge subjected to traveling seismic waves,” Engineering Structures, vol. 33, no. 12, pp. 3546–3558, 2011. View at: Publisher Site  Google Scholar
 Z. P. Zeng, F. S. Liu, P. Lou, Y. G. Zhao, and L. Peng, “Formulation of threedimensional equations of motion for train–slab track–bridge interaction system and its application to random vibration analysis,” Applied Mathematical Modelling, vol. 40, no. 1112, pp. 5891–5929, 2016. View at: Publisher Site  Google Scholar
 Z. C. Zhang, J. H. Lin, Y. H. Zhang, Y. Zhao, W. P. Howson, and F. W. Williams, “Nonstationary random vibration analysis for trainbridge systems subjected to horizontal earthquakes,” Engineering Structures, vol. 32, no. 11, pp. 3571–3582, 2010. View at: Publisher Site  Google Scholar
 Z.W. Yu, J.F. Mao, F.Q. Guo, and W. Guo, “Nonstationary random vibration analysis of a 3D trainbridge system using the probability density evolution method,” Journal of Sound and Vibration, vol. 366, pp. 173–189, 2016. View at: Publisher Site  Google Scholar
 J. Li and J. B. Chen, Stochastic Dynamics of Structures, John Wiley & Sons (Asia) Ptd Ltd, Singapore, 2009.
 L. Xu, W. Zhai, and J. Gao, “A probabilistic model for track random irregularities in vehicle/track coupled dynamics,” Applied Mathematical Modelling, vol. 51, pp. 145–158, 2017. View at: Publisher Site  Google Scholar
 L. Xu and W. Zhai, “A new model for temporalspatial stochastic analysis of vehicletrack coupled systems,” Vehicle System Dynamics, vol. 55, no. 3, pp. 427–448, 2017. View at: Publisher Site  Google Scholar
 J. M. Rocha, A. A. Henriques, and R. Calçada, “Probabilistic safety assessment of a short span highspeed railway bridge,” Engineering Structures, vol. 71, pp. 99–111, 2014. View at: Publisher Site  Google Scholar
 J. M. Rocha, A. A. Henriques, R. Calçada, and A. Rönnquist, “Efficient methodology for the probabilistic safety assessment of highspeed railway bridges,” Engineering Structures, vol. 101, pp. 138–149, 2015. View at: Publisher Site  Google Scholar
 J. M. Rocha, A. A. Henriques, and R. Calçada, “Safety assessment of a short span railway bridge for highspeed traffic using simulation techniques,” Engineering Structures, vol. 40, pp. 141–154, 2012. View at: Publisher Site  Google Scholar
 Z. Jin, S. Pei, X. Li, and S. Qiang, “Probabilistic evaluation approach for nonlinear vehiclebridge dynamic performances,” Journal of Sound and Vibration, vol. 339, pp. 143–156, 2015. View at: Publisher Site  Google Scholar
 T. Cho, M. K. Song, and D. H. Lee, “Reliability analysis for the uncertainties in vehicle and highspeed railway bridge system based on an improved response surface method for nonlinear limit states,” Nonlinear Dynamics, vol. 59, no. 12, pp. 1–17, 2010. View at: Publisher Site  Google Scholar
 P. Antolín, N. Zhang, J. M. Goicolea, H. Xia, M. Á. Astiz, and J. Oliva, “Consideration of nonlinear wheelrail contact forces for dynamic vehiclebridge interaction in highspeed railways,” Journal of Sound and Vibration, vol. 332, no. 5, pp. 1231–1251, 2013. View at: Publisher Site  Google Scholar
 Q. Li, Y. L. Xu, D. J. Wu, and Z. W. Chen, “Computeraided nonlinear vehiclebridge interaction analysis,” Journal of Vibration and Control, vol. 16, no. 12, pp. 1791–1816, 2010. View at: Publisher Site  Google Scholar
 J. Li and J.B. Chen, “The probability density evolution method for dynamic response analysis of nonlinear stochastic structures,” International Journal for Numerical Methods in Engineering, vol. 65, no. 6, pp. 882–903, 2006. View at: Publisher Site  Google Scholar
 K. JJ, On the Rolling Contact of Two Elastic Bodies in the Presence of Dry Friction, T H Delft, Delft, Netherlands, 1967.
 S. Z. Meymand, A. Keylin, and M. Ahmadian, “A survey of wheelrail contact models for rail vehicles,” Vehicle System Dynamics, vol. 54, no. 3, pp. 386–428, 2016. View at: Publisher Site  Google Scholar
 J. Liang, S. R. Chaudhuri, and M. Shinozuka, “Simulation of nonstationary stochastic processes by spectral representation,” Journal of Engineering Mechanics, vol. 133, no. 6, pp. 616–627, 2007. View at: Google Scholar
 Z. Liu, W. Liu, and Y. Peng, “Random function based spectral representation of stationary and nonstationary stochastic processes,” Probabilistic Engineering Mechanics, vol. 45, pp. 115–126, 2016. View at: Publisher Site  Google Scholar
 National Railway Administration, PSD of Ballastless Track Irregularities of HighSpeed Railway.(TB/T 33522014), National Railway Administration, Beijing, China, 2014.
Copyright
Copyright © 2019 FuShan Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.