Safety evaluation of a vehicle–bridge interaction system using the pseudo-excitation method

A method for analysing the vehicle–bridge interaction system with enhanced objectivity is proposed in the paper, which considers the time-variant and random characteristics and allows finding the power spectral densities (PSDs) of the system responses directly from the PSD of track irregularity. The pseudo-excitation method is adopted in the proposed framework, where the vehicle is modelled as a rigid body and the bridge is modelled using the finite element method. The vertical and lateral wheel–rail pseudo-excitations are established assuming the wheel and rail have the same displacement and using the simplified Kalker creep theory, respectively. The power spectrum function of vehicle and bridge responses is calculated by history integral. Based on the dynamic responses from the deterministic and random analyses of the interaction system, and the probability density functions for three safety factors (derailment coefficient, wheel unloading rate, and lateral wheel axle force) are obtained, and the probabilities of the safety factors exceeding the given limits are calculated. The proposed method is validated by Monte Carlo simulations using a case study of a high-speed train running over a bridge with five simply supported spans and four piers.


Introduction
The viaduct bridges are widely used in high-speed railway systems, in order to control the foundation settlement and to reduce the workload in railway maintenance. As a critical problem in bridge design, the dynamic effects of live loads must be considered, which include the moving effects of the train and the track irregularity. The vehicle-bridge interaction system is a popular model to optimise the design parameters of the vehicle and the bridge. However, since the parameters used in calculations usually differ from the actual values [1], probabilistic analysis is necessary for the vehicle-bridge interaction system. As a deterministic and time-variant system subject to random excitations, it is important to better understand the random responses of the running safety indices, including the derailment coefficient, the wheel unloading rate, and the lateral wheel axle force, so as to ensure the running safety and to enhance the reliability and cost-effectiveness of bridge design.
Parameters of the vehicle, bridge and track structure can be regarded as certain deterministic factors, while the track irregularity, as one of the system excitations, is a random factor and expressed by known power spectrum formulas in most cases. Traditionally, the track irregularity is input into the interaction system model as a spatial sample, which is numerically simulated, such as the harmony superposition method (HSM) [2,3]. The spatial sample is necessarily random; thus, the dynamic response of the vehicle-bridge interaction system is also a stochastic process and the maximum value of the factors calculated from spatial irregularity samples has inevitable variability. Some attempt for statistical parameters of the dynamic response is carried out by simulating many times with different random phase angles in HSM. The work improved the reliability of calculation but its error is still difficult to estimate.
The problem of random response of the vehicle-bridge interaction system has attracted attention of many researches. Perrin et al. [4] established a nonstationary four-dimensional random field model of track irregularities and obtained the probability density function of responses. Wetzel et al. [5] studied the random dynamic responses of trains under wind load by a subset simulation method. Majka et al. [6] evaluated the safety of the vehicle-bridge interaction system under random track irregularities considering bridge skewness. Wu et al. [7] described the bridge using non-Gaussian uncertainty model and investigated the dynamic criteria by the method of polynomial chaos. Lombaert et al. [8] solved the interaction problem in the frequency domain and studied the effect of the constant moving load under road surface unevenness by decomposing the vehicle-bridge interaction force using the Fourier transform. For other types of excitations present in the interaction system, Kiureghian et al. [9] adopted random structural vibrations produced by multi-point seismic excitation, and Alduse et al. [10] considered the uncertainty of wind speed and direction to study the fatigue damage in long-span bridges.
The classical random vibration theory usually adopts traditional sampling methods, which without exception require a large size of samples and are therefore computationally ineffective. Therefore, other mathematical methods are used to obtain the power spectra of the dynamic responses directly from the input power spectrum of track irregularity. In recent years, the pseudo-excitation method (PEM) [11] has become the most popular of these methods. Zhang et al. [12] improved the train riding comfort condition by optimising the vehicle suspension parameters with the min-max approach. Zhu et al. [13] proposed an approach for predicting the train-induced ground vibrations considering random track unevenness. Li et al. [14] and Zhu et al. [15] proved that the bridge responses follow Gaussian distributions and established a framework to calculate the upper and lower limits of vehicle-bridge system responses from the auto-spectra and the cross-spectra of bridge vibrations.
The multiple-excitation analysis of vehicle-bridge interaction system can be established by considering the spectrum of the wind load and the seismic load. Zhang et al. [16] analysed the nonstationary random responses of interaction systems subjected to lateral horizontal earthquakes by PEM and the precise integration method (PIM) and obtained the time-dependent power spectral density (PSD) functions and standard deviations of the responses. He et al. [17] proposed an efficient analysis framework for the interaction system subjected to wind loads and studied the nonstationary effects on the vibration responses. The methods presented in [10][11][12][13][14][15][16] were all validated by Monte Carlo simulations. By the above-mentioned studies, the statistic characteristics of response in frequency domain can be fully obtained, free of the accidental effect; also the dynamic behaviour in each frequency band of the system movements and forces can be obtained. The work improves the objectivity of bridge dynamic analysis.
The relative position of the vehicle with respect to bridge subsystems is ever-changing when the train is running over the bridge, and thus, the interaction system is time variant. However, PEM requires each subsystem and their interactions to be linear. For solving the problem, Zhai et al. [1] defined different time steps for the subsystems of the vehicle, the track, and the bridge, in order to meet the accuracy requirements. Lei et al. [18] proposed an iterative procedure in time and avoided solving the nonlinear equations when a nonlinear wheel-rail interaction model was adopted. Zhang et al. [19] established an interhistory iteration method by which the two subsystems were solved separately. The iterative process was interposed for faster convergence. Liu et al. [20] developed an easier iteration method, in which the wheel-rail interaction force is defined by the system response from the previous time step. Zhu et al. [13] divided the vehicle-track-bridge system into the vehicle-track subsystem and the bridge subsystem and improved the accuracy by using a much lower interface stiffness between the two subsystems. Mohammdzadeh et al. [21] simulated the wheel-train interaction with SIMPACK software and assembled the two subsystems using the importance sampling and response surface methods. It should be noted that the calculation efficiency was suboptimal in the former research; for example, it is reported by He et al. [17] that a single case of PEM analysis required 900 s of computation time, which is much longer than that of the traditional vehicle-bridgewind interaction analysis.
It is known that a train has high risk of derailment while passing over a bridge if the stiffness of the bridge is insufficient, with which the safety factors defined in the design code TB 10002-2017 (Code on Design of Railway Bridges and Culverts) are concerned. Clearly, the safety factors are random with certain probability density functions, and it is required that the maximum values of the safety factors are below the limits stipulated in the code. A safer and more effective framework for vehicle-bridge interaction system assessment is proposed in this paper, in which the dynamic responses are calculated by PEM, the probability of the safety factors exceeding the given limits can be calculated, and the vehicle-bridge interaction system performance can be evaluated. The Monte Carlo method is used to validate the proposed method using a case study of a bridge with five simply supported spans.

Theoretical analysis method for vehicle-bridge interaction system based on PEM
The load acting on the vehicle-bridge interaction system has two parts: the deterministic one and the random one. The deterministic load is the moving load of the train, for which the system response can be found by the existing methods [19] without considering the track irregularities. The random load is induced by the dynamic effects of the track irregularity, for which PEM can be adopted. PEM is a popular method in random vibration analysis [11], by which the power spectrum of system response can be obtained directly from the power spectrum of excitation. The track irregularity is a zero-mean Gaussian process; thus, the dynamic response due to track irregularity must be another zero-mean Gaussian process, too. In other words, the mean value of the response is determined by the moving load of the train and the standard deviation of the response is determined by the track irregularity. It can be proved that the standard deviation, r i , of response for the ith degree of freedom of the system can be obtained from the auto-power spectrum S ii (x, t) as follows: where x is frequency and t is time. Then, the probabilistic characteristics of the system response can also be obtained. The flowchart of the solution procedure is shown in Fig. 1.
The following assumptions [19] are adopted in the modelling of the vehicle-bridge interaction system: A1: The vehicle and bridge subsystems are linear. A2: The vertical wheel-rail interaction force is derived from the assumption that the wheel and the rail have the same vertical displacement at any time. A3: The lateral wheel-rail interaction force is derived from the simplified Kalker creep theory, i.e. the lateral wheel-rail interaction force is proportional to the wheel-rail relative velocity, and the proportionality coefficient is defined in [19]. A4: The wheel-rail interaction force is distributed to the two neighbouring bridge nodes inversely proportional to the wheel-node distances.
The vehicle-bridge interaction system is illustrated in Fig. 2. The equations of motion for the two interacting subsystems can be expressed as follows: where M, C, and K are the mass, damping, and stiffness matrices, respectively; X and F are the displacement and force vectors; the subscripts 'v' and 'b' refer to the vehicle and the bridge, respectively. The dynamic matrices of the vehicle and the bridge are derived, respectively, from the Safety evaluation of a vehicle-bridge interaction system using the pseudo-excitation method 43 rigid body dynamics method and finite element method (FEM) [19], where the damping matrix of the bridge subsystem C b is formed by the Rayleigh damping, and the proportional factors are calculated by the first two frequencies of the bridge model. As an interaction system, the right-hand side terms, F v and F b , in Eq. (2) are functions of X v and X b of the wheel sets, as expressed in Eq. (3): where R stands for summation over the wheel sets; m, c, k, and f stand for the mass, damping, stiffness and irregularity-induced forces, respectively, of a single wheel set; quantities with subscripts 'vv' and 'bb' are for the vehicle and bridge only, respectively; quantities with subscripts 'vb' and 'bv' are attributed to the interaction between the vehicle and bridge; subscripts 'vi' and 'bi' stand for the vehicle and bridge subsystems related to vehicle i, respectively. Then, the dynamic equilibrium equations of the vehicle-bridge interaction system can be formed by moving the right-hand side terms in Eq. (2) to the left-hand side as in Eq. (4).
For a given wheel set, its lateral (y) displacement is denoted as y 1 ; the relative vertical (z) displacement, and pitch (about axis y) and yaw (about axis x) angles of the vehicle bogie are denoted as z 1 , u 1 , and v 1 , respectively; for the neighbouring nodes I and J, their lateral displacements are denoted as y 2 and y 3, respectively, vertical displacements as z 2 and z 3, respectively, and torsional rotations as u 2 and u 3, respectively. The wheel-rail interaction force is distributed to the two neighbouring bridge nodes inversely proportional to the wheel-node distances, where the weighting factors are denoted as p for node I and q for node J, respectively. For a given wheel set, the equilibrium equations of the interaction system subjected to pseudo-excitation are expressed as follows: Wheel-rail interaction Track irregularity where Eq. (5) describes the vertical pseudo-excitation, Eq. (6) the torsional pseudo-excitation, Eq. (7) the lateral pseudo-excitation, and Eq. (8) the pseudo-inertia excitation, respectively; definitions of the symbols used in these interaction equations are listed in Table 1.
The power spectra of the vertical, lateral, and torsional pseudo-excitations are as follows: where V is the train speed; S dd (x), S ee (x), and S bb (x) are the vertical, lateral, and torsional power spectra of excitations; X is the spatial frequency of track irregularity; S zz (X), S yy (X), and S uu (X) are the power spectra of track irregularity in the vertical, lateral, and torsional directions [14,15]. From Eqs. (3)(4)(5)(6)(7)(8)(9)(10)(11), the additional loads of a single wheel set related to the pseudo-excitations caused by irregularities d, e, and b are respectively as follows: ð12Þ where Dt is the time difference between the concerned wheel set and the first wheel set passing over the same point, i is the imaginary unit and equal to the square root of -1. The pseudo-excitation of irregularities can be obtained by adding the pseudo-excitations for all the wheel sets. Based on the fundamental assumption of PEM, the track irregularities related to vertical direction, d, torsional angle, b, and lateral direction, e, must be analysed independently.
In this way, the power spectra of vehicle-bridge interaction system can be calculated. For a single wheel set, the force of the primary suspension system under pseudo-excitation related to directions y and z, and rotation u (about axis x) are as follows: The vertical acceleration, a w , and angular acceleration about axis x, b w , are as follows: Thus, the interaction forces at the left-and right-hand side wheel-rail contact points are as follows: where F wy1 and F wy2 are the lateral contact forces at the left-and right-hand sides, respectively; F wz1 and F wz2 are the lateral contact forces at the left-and right-hand sides, respectively; G is the static load of wheel set; and g 0 is the gauge.
The maximum values of the derailment coefficient, X DR , the wheel unloading rate, X OL , and the lateral wheel axle force, F AY , are required to be evaluated as safety factors in the designing codes [2]. They are functions of the wheelrail forces as follows: where P 0 is the static wheel load, which is a constant; P(t) and Q(t) are the vertical and lateral wheel-rail forces at a single side of the wheel set, respectively. It can be  Irregularity indicator vector related to torsional inertia pseudoexcitation proved by Eqs. (5-8) that Q(t) and P(t) are both an independent random process. As a random process, the probability of a factor exceeding its limit can be calculated from its probability density function. The mean values of the wheel-rail vertical and lateral forces, l P and l Q , can be obtained using the deterministic loads and the standard deviation can be calculated by Eq. (1). Then, it can be proved by the central limit theorem that the wheel-rail vertical and lateral forces follow Gaussian distributions when the number of samples is large enough.
The wheel unloading rates and the lateral wheel axle forces are functions of the wheel-rail vertical and lateral forces; thus, they will also follow Gaussian distributions with the following mean values and standard deviations: where l OL and r OL are the mean value and standard deviation of the wheel unloading rate, respectively; l FAY and r FAY are the mean value are standard deviation of the lateral wheel axle force, respectively. The derailment coefficient is the ratio of two random variables, P(t) and Q(t), both following Gaussian distributions with known mean values and standard deviations. Thus, the probability density function of the derailment coefficient is as follows: where U Á ð Þ is the cumulative distribution function of the standard Gaussian distribution, and Then, the probability of the concerned factor exceeding its limit can be calculated by integrating the probability density function.

Case study
In this case study, a typical high-speed train runs over a bridge with four column piers and five pre-stressed concrete simply supported spans. A mass density per metre and frequencies are assigned to the span beams instead of a specified cross section. The nodes at the abutments and pier bottoms are fixed in the model. The train comprises eight vehicles, each 25 m in length, and runs at speeds of 250, 300, and 350 km/h. To reduce the effect of sudden load application, the first wheel set of the train is assumed to be 100 m from the bridge end when the simulations starts, i.e. the train is on the bridge between 1.03 and 3.03 s for 350 km/h. The bridge is shown in Fig. 3, and a complete list of the vehicle and bridge parameters are provided in Tables 2 and 3, respectively. The power spectra of track irregularity, S(f), are adopted from TB/T 3352-2014 (PSD of Ballastless Track Irregularities of High-Speed Railway) in the directions of profile (z direction), alignment (y direction), and cross-lever (about x axis) as follows: where S(f) is in mm 2 /m -1 ; f is spatial frequency in m -1 ; A and k are coefficients. The values of coefficients A and k and spatial frequency f are listed in Table 4 and Table 5, respectively. The cross-lever is the height difference of the two rails, which is the product of the torsional irregularity and the gauge.
To validate the proposed method, 100 Monte Carlo simulations were carried out with track irregularity samples generated by HSM. A comparison of the PEM and Monte Carlo simulation results for the train speed of 350 km/h is shown in Figs. 4-11, where the black lines represent the results of PEM, and the red lines represent the mean values or standard deviations of the Monte Carlo simulations. (In graphs throughout the paper, abbreviations Disp., Acc., and AA stand for the displacement, acceleration, and angular acceleration, respectively.) It is found that the two methods agree well for all the concerned factors.
However, some noticeable differences are also found between the PEM and Monte Carlo results, especially for the factors with zero-mean value. This is due to the limited sample numbers used in the Monte Carlo simulations. Larger sample numbers were tried and it was found that the       In order to analyse the distribution of the safety factors, the results produced by the proposed method were compared to those by the Monte Carlo method, as illustrated in Table 6 and Figs. 16-18, where the numbers of samples in    It is clear that the probability of the safety factors exceeding the limits increased with the analysis duration, which may lead to loss of objectivity. The limits concerned in this paper are from the codes for evaluating the safety performance based on the maximum values measured in a single experiment or calculation. Thus, it is not the original intention of the code to ensure that the safety factors are less than the limits for every train over the entire service life of the bridge. The probability density functions of the three safety factors are studied in this paper, based on which reasonable limits can be probabilistically established in future work.
10 l ? 4r ? F -1 (U(4)) ? 10 DR, OL, and FAY stand for the derailment coefficient, the wheel unloading rate, and the lateral wheel axel force, respectively. F(x) is the probability distribution function of the derailment coefficient obtained from its probability density function (Eq. (27)) by integration, and U(Á) is the cumulative distribution function of the standard Gaussian distribution. (1) Despite its time variety and randomness, the vehiclebridge interaction problem can be solved by the method proposed in the paper based on PEM, by which the power spectrum densities of vehicle-bridge interaction system responses can be obtained directly from the power spectrum density of the track irregularity. The proposed method has been validated by Monte Carlo simulations. (2) A typical Chinese high-speed railway bridge has been adopted as a case study, in which the dynamic responses of the vehicle, the bridge, and wheel-rail interaction forces were calculated by the proposed method. It has been found that the dominant frequency of bridge vertical response is the load frequency of train, while the dominant frequencies of bridge lateral responses are the natural frequencies of bridge. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/. Safety evaluation of a vehicle-bridge interaction system using the pseudo-excitation method 55