# Formulation of a Linearized Model for Vehicle-Track Interactions

Authors: Fei Huo; Huanyun Dai
Journal: Shock and Vibration
 (2019)

---

## Abstract

In this paper, a fundamental wheel-rail interaction (WRI) element accompanied by its coupling matrices with other vehicle-track components have been derived taking into consideration the aspect of linearization. The key to the presented formulation is the use of the geometrical relationships of relative motions between degrees of freedom (DOFs) and energy principle. To the WRI element, both of the conditions of wheel-rail contacts and wheel-rail separations are allowed in the numerical computations; besides, the effects of the linear creepage and the gravitational restoring are considered in the description of wheel-rail interactions. By comparing with an advanced three-dimensional nonlinear model, the capability of the linear model in characterizing the response amplitude and frequency characteristic of vehicle-track systems is demonstrated. Moreover, the method for the random vibration analysis of the linear model is presented by treating the creep coefficients as the random sources, through which the safety margin of system response can be predicted well. From the numerical examples, it is, additionally, concluded that the lateral creep coefficient holds significant influence on wheel-rail lateral interactions and track vibrations, especially for the responses at low frequency ranges.


---

## Body

## 1. Introduction

To ensure a comfortable, safe, and efficient transportation, it was critical to scientifically assess the dynamic performance of vehicle-track systems under varying loads; consequently, a vehicle-track interaction system was formulated. The vehicle-track interactions are essentially nonlinear, while for lowering the modelling complexity, increasing the computational efficiency, and broadening the model application, it was rather popular to simplify the vehicle-track systems as a coupled and linearized system.In the last decades, the models aiming at describing the vehicle-track interactions have experienced a rapid development from the earliest two-dimensionally vertical models [1–8] to more advanced three-dimensionally spatial ones. Concerning the system nonlinearity, Sun et al. [9] presented a 3D wagon-track system dynamics model, where a wagon system with 37 degrees of freedom (DOFs), a four-layer track, and the wheel-rail interfacial contacts were comprehensively considered; Zhai et al. [10] presented a theoretical framework and methodologies for characterizing the 3D nonlinear vehicle-track interactions; recently, Xu and Zhai [11, 12] made pioneering and the most advanced work in developing temporal-spatial stochastic models for train-track (bridge) interactions, in which the finite element method (FEM) and random theory were introduced to construct the track systems and to achieve the random vibration analysis, respectively.With the development of nonlinear models, the methodologies are constantly developed to mathematically investigate the vehicle-track (bridge) interactions in the field of linearization modelling. See for instance, Zeng and Guo [13] made pioneering work on building a train-bridge time-dependent model by an application of the principle of a stationary value of total potential energy of dynamic system [14]; Yang and Wu [15] proposed a versatile element that is capable of considering various vehicle-bridge interaction (VBI) effects; however, the condition of wheel-rail separation is neglected in both of the these two models. Some other related researches can be consulted in [16, 17]. Moreover, Zhang et al. [18] regarded the vehicle as a spring-mass-damper system and the track as an infinitely long substructural chain consisting three layers of rail, sleeper, and ballast, and the vehicle and the tracks were coupled by linear springs; based on this developed model, the pseudoexcitation method (PEM) was introduced to achieve the high-efficient frequency characteristic analysis for system indices. Nguyen et al. [19] presented simulations for vehicle-substructure dynamic interactions and wheel movements, where the two compatibilities, i.e., force equilibrium and geometrical relation, are satisfied in the wheel-rail interfaces.Additionally, the train-track dynamic simulations under transient conditions also show remarkable significance to railway maintenance, assessment, rolling fatigue, safety, etc. See for instance, Handoko and Dhanasekar [20] pointed out that the traction and braking forces are seldom considered in the practice of wagon dynamics simulation, although these forces will greatly modify the wheel-rail contact parameters and then the wheelset dynamics; Liu et al. [21] also studied the mechanism of wheelset longitudinal vibration by analysing the process of wheel/rail rolling contact; generally, the track defects can cause profound effects to the dynamics of the railway wagon; Zhang and Dhanasekar [22] noticed this problem and published a model for the dynamics of wagons subject to braking/traction torques on a perfect track by explicitly considering the pitch degree of freedom for wheelsets [23, 24] and extended this model for cases of lateral and vertical track geometry defects and worn railhead and wheel profiles; Grossoni et al. [25] examined the dynamic behaviour at a rail joint using a two-dimensional vehicle-track coupling model, where the influence of the number of layers and the number of elements between two sleepers and the beam model are investigated. Laterally, Zong and Dhanasekar [26] considered the gap between rail joints to account for thermal movement and to maintain electrical insulation for the control of signals and/or broken rail detection circuits; besides, railhead can provide high stresses due to the passage of heavily loaded wheels through a very small contact patch [27], and a multibody dynamic model was developed in [28] to accurately model and analyse the track dynamic behaviour in vicinity of rail discontinuities; recently Zong and Dhanasekar [29] provided a idea of simplifying the design of the IRJs consisting of only two pieces of insulated rails embedded into a concrete sleeper; Ling et al. [30] presented a formulation for a passive road-rail crossing involving stiffened edges of the raised road pavement to minimise the risk of failure of wheel-rail contact using a nonlinear three-dimensional multibody dynamics model; additionally, a series of work on impact derailment due to lateral collisions between heavy road vehicles and passenger trains at level crossings and the associated derailments had also been conducted in [31–34].In the above presentations, though the nonlinear and linear vehicle-track interaction, models have been established by many researchers, while the applicability of the linear model and its comparison to the nonlinear model have rarely been reported. In this paper, a contribution, which is not beneficial to train-track dynamic to severe transient impact/rolling issues, was made to develop a highly practical and efficient model, and fully deriving the equation of motion for the linearized vehicle-track interactions. The coupling mechanism between a railway vehicle and the tracks in relation to wheel-rail interactions and system responses that are of great interest to railway engineers will be revealed.The organization of this paper is listed as follows:(1)
In section2, a general wheel-rail interaction (WRI) element is presented, in which the jumping of the vehicle off the rail, the linear creepage, and gravitational restoring actions are properly considered.(2)
In section3, the coupling matrices for building the physical-mechanical connections between system components and WRI element are elaborately illustrated and finally the unified equation of motion for the vehicle-track system is established.(3)
In section4, numerical examples, including comparisons with the nonlinear model and a practice of random vibration analysis, are conducted to validate the reliability and feasibility of the proposed model.(4)
In section5, some conclusions are drawn accordingly.


## 2. WRI Element

In the whole vehicle-track model presented in Figure1, the wheel-rail interactions play a key role in determining the dynamic trajectories of vehicle/track behaviours. In this part, a wheel-rail interaction (WRI) element will be developed to depict the wheel-rail interaction mechanism from aspect of linearization.Figure 1
Model, dimensions, and parameters of a vehicle. (a) Side view; (b) front view; (c) top view; (d) sign convention of degree of freedoms.
(a)
(b)
(c)
(d)


### 2.1. Displacement Vector

For a wheel-rail contact pair, as shown in Figure2, the displacement vector, X˜wr,k, with order of T˜wr×1, can be expressed by(1)X˜wr,k=xw,kyw,kzw,kψw,kLv,kRv,kLz,kRz,kLh,kRh,kLy,kRy,kLn,kRn,kLx,kRx,kT,where the superscript “T” denotes the transposition of vectors or matrices; the subscript “k” denotes the kth wheel-rail contact element; xw, yw, zw, and ψw are the longitudinal, lateral, vertical, and yaw motions of the wheel; by assuming the track irregularity as virtual displacement of the rail at the wheel-rail contact positions, Lv and Rv are the vertical irregularities of the rail at the left and right side, respectively, Lh and Rh are the lateral irregularities of the rail at the left and right side, respectively; moreover, Lz and Rz, respectively, denote the vertical displacement of the left- and right-side rail at the wheel-rail contact points; in the same manner, Ly and Ry denote the lateral displacement of the left- and right-side rail at the wheel-rail contact points, respectively; Ln and Rn denote the angle of torsion of the left- and right-side rail around the Y axis, respectively; and Lx and Rx denote the longitudinal displacement of the rails with respect to the left and right wheel-rail contact points, respectively.Figure 2
3D random wheel/rail contact model: (a) front view; (b) elliptical contact spot at location of C.
(a)
(b)In normal wheel-rail contacts, it is assumed that the vertical motion of the wheel is nonindependent and being totally determined by the rail vertical displacement, namely,(2)zw,k=12Lv,k+Rv,k+Lz,k+Rz,k,with(3)Lz,k=DR,zl,kNR,z,kT,Rz,k=DR,zr,kNR,z,kT,DR,zl,k=WR1,lθR1,lyWR2,lθR2,ly,DR,zr,k=WR1,rθR1,ryWR2,rθR2,ry,Nr,1=1−3ζw,k/lwr2+2ζw,k/lwr3,Nr,2=ζw,k1−2ζw,k/lwr+ζw,k/lwr2,Nr,3=3ζw,k/lwr2−2ζw,k/lwr3,Nr,4=ζw,kζw,k/lwr2−ζw,k/lwr,NR,z,k=Nr,1Nr,2Nr,3Nr,4,where DR,zl and DR,zr are the nodal displacement vector used to calculate the vertical displacement of the rail beam at the left and right side along the X axis; W and θy are the nodal vertical displacement and angular displacement of the node around the Y axis, in which the subscripts “l” and “r” denote the left- and right-side rails, respectively and the subscripts “R1” and “R2” denote the left and right nodes with respect to the rail beam element, respectively; ζw denotes the distance between the left node of the rail element and the wheel-rail contact position; and NR,z is the Hermite polynomial interpolation used for Euler–Bernoulli beam.When the wheels lose contact with the rails, the vertical motion of the wheel is independent, by assumingzw in Equation (1) as an independent degree of freedom (DOF). In the numerical simulation, the inequalities below are applied to judge whether or not the wheelsets are contacting to the rails [35] (set the left wheel-rail contact side as an example):(4)Lv+Lz−zwminLv+Lz,zw<ξ1,1−ξ2≤Lv+Lzzw≤1,whenzw>0,Lv+Lz>0,1−ξ3≤zwLv+Lz≤1,whenzw<0,Lv+Lz<0,maxLv+Lz,zw≤ξ4,whenzw<0,Lv+Lz>0,where ξi, i=1,2,3,4, is a small positive number depending on the analytical precision.


### 2.2. Wheel-Rail Coupling Matrices

In most of the running conditions, the wheels constantly keep contact with the rails, thus the wheel-rail coupling matrices that consider the wheelset’s vertical displacement as nonindependent will be presented here as an emphasis. The wheel-rail coupling matrices are marked by subscript “wr” with order ofT˜wr×T˜wr.The mass submatrix of WRI element can be written as(5)M˜wr=∑k=1nwrMwr,kΘm,k,Θm,k+Mwr,kℓxw,k,ℓxw,k+Mwr,kℓyw,k,ℓyw,k+Mwr,kℓψw,k,ℓψw,k,=∑k=1nwrmw,z,kNmTNm+mw,y,k+mw,x,k+Iwz,k,with Nm=1/21/21/21/2, where mw,z, mw,y, and mw,x denote the mass of the wheelset against the vertical, lateral, and longitudinal motion of the wheelset; Iwz is the mass moment of inertia of the wheelset around the Z axis; Θm is the number of the freedoms in Equation (1) related to the mass submatrix, Θm=ℓLvRvLzRz, ℓ⋅ denotes the operator characterizing the corresponding number of DOF in the displacement vector, and nwr=4 for a four-axle vehicle.The damping submatrix of WRI element is mainly induced by the virtual work of wheel-rail creep forces, which can be derived as(6)C˜wr=∑k=1nwrCwr,1l,kΘc,1l,k,Θc,1l,k+Cwr,1r,kΘc,1r,k,Θc,1r,k+Cwr,2l,kΘc,2l,k,Θc,2l,k+Cwr,2r,kΘc,2r,k,Θc,2r,k=∑k=1nwrCXV¯Nc1,lTNc1,l+Nc1,rTNc1,r+CYV¯Nc2,lTNc2,l+Nc2,rTNc2,r,with(7)Nc1,l=1llr−1,Nc1,r=1−llr−1,Nc2,l=1−1−1lcp,Nc2,r=1−1−1lcp,Θc,1l,k=ℓxw,kψw,kLx,k,Θc,1r,k=ℓxw,kψw,kRx,k,Θc,2l,k=ℓyw,kLh,kLy,kLn,k,Θc,2r,k=ℓyw,kRh,kRy,kRn,k,where CX and CY denote the longitudinal and lateral creep coefficient, respectively; V¯ is the forward speed of the wheelset nominally; Cwr,1l and Cwr,1r denote the damping matrices derived by the virtual work of wheel-rail longitudinal creep forces; Cwr,2l and Cwr,2r denote the damping matrices derived by the virtual work of wheel-rail lateral creep forces; llr is the lateral half distance the two wheel-rail contact points at the left and right side, and lcp is the vertical distance between the wheel-rail contact point and the centroid of the wheelset.The stiffness submatrix of WRI element is mainly derived from the work of wheel-rail equivalent gravity stiffness, namely,(8)K˜wr=∑k=1nwrKwr,kΘs,k,Θs,k=∑n=1nwrkg,kNs,kTNs,k,with(9)Θs,k=ℓyw,kLh,kRh,kLy,kRy,k,Ns,k=ℓ1−1/2−1/2−1/2−1/2,where kg,k is the wheel-rail equivalent gravity stiffness, which is calculated using the following equation [36]:(10)kg=Wbϕςδ0+ε,with(11)ε=bR−R′ϕ11+Rδ0b,ς′=RR−R′ϕ11+R′δ0b,ϕ1=1−r0δ0b,where ε is the angular parameter of the wheel-rail contact, ς is the change rate of the lateral relative displacement between the wheel centroid and wheel-rail contact point, b is the half of lateral distance between two rolling circles, R and R′ are the radii of curvature of the wheel and the rail at the contact point, respectively, δ0 is the wheel-rail contact angle when the wheelset is placed in the central position, and r0 is the radius of the nominal rolling circle.The load vector of WRI element is originated from the gravity of the wheelset, namely,(12)F˜wr=∑k=1nwrFwr,kΘfl,k,Θfl,k+Fwr,kΘfr,k,Θfr,k=∑k=1nwr12mw,kgNfl+12mw,kgNfr,with (13)Θfl,k=ℓLv,kLz,k,Θfr,k=ℓRv,kRz,k,Nfl=Nfr=11,where g is the gravitational acceleration.Equations (5), (6), and (8) have briefly presented the mass, damping, and stiffness matrices of WRI element, namely, M˜wr, C˜wr, and K˜wr all with order of nwrT˜wr×nwrT˜wr. Moreover, it should be noted that the some of the DOFs in Equation (1), i.e., Lz,k, Rz,k, Ly,k, Ry,k, Ln,k, Rn,k, Lx,k, and Rx,k, which can be called as the interior node, are the interpolation results of the nodal displacements of the rail beam element. Besides Lz,k=DR,zlNR,zT and Rz,k=DR,zrNR,zT illustrated in Equation (2), the other DOFs can be equivalently transformed by the nodal displacement as(14)Ly,k=DR,ylNR,yT,Ry,k=DR,yrNR,yT,with(15)DR,yl=VR1,lθR1,lzVR2,lθR2,lz,DR,yr=VR1,rθR1,rzVR2,rθR2,rz,NR,y=Nr,1−Nr,2Nr,3−Nr,4,where V denotes the lateral displacement of the rail node, θz denotes the angular displacement of the rail node around the Z axis.Moreover,(16)Ln,k=DR,nlNR,nT,Rn,k=DR,nrNR,nT,with(17)DR,nl=θR1,lxθR2,lx,DR,nr=θR1,rxθR2,rx,NR,n=1−ζw,k/lwrζw,k/lwr,where θx denotes the angular displacement of the rail node around the X axis:(18)Lx,k=DR,xlNR,xT,Rx,k=DR,xrNR,xT,with(19)DR,xl=UR1,lxUR2,lx,DR,xr=UR1,rxUR2,rx,NR,x=NR,n,where Ux denotes the longitudinal displacement of the rail node.Based on Equations (2), (14), (16), and (18), the displacement vector, X˜wr,k, presented in Equation (1) can be further extended as(20)X˜wr,k=xw,kyw,kzw,kψw,kTLv,kRv,kWR1,l,kθR1,ly,kWR2,l,kθR2,ly,kWR1,r,kθR1,ry,kWR2,r,kθR2,ry,kTLh,kRh,kVR1,l,kθR1,lz,kVR2,l,kθR2,lz,kVR1,r,kθR1,rz,kVR2,r,kθR2,rz,kTθR1,lx,kθR2,lx,kθR1,rx,kθR2,rx,kTUR1,lx,kUR2,lx,kUR1,rx,kUR2,rx,kT.In the meantime, the dynamic matrices and load vector over WRI element can be conveniently redistributed into this extensive nodal displacement vector; for example, the newly developed mass matrixMwr with order Twr×Twr can be obtained by(21)Mwr=M˜wrΘm,ℓDR,zr+M˜wrΘm,ℓDR,zl+M˜wrℓDR,zr,Θm+M˜wrℓDR,zl,Θm,with(22)M˜wrΘm,ℓDR,zr=M˜wrΘm,ℓRzNR,z,M˜wrΘm,ℓDR,zl=M˜wrΘm,ℓLzNR,z,M˜wrℓDR,zr,Θm=NR,zTM˜wrℓRz,Θm,M˜wrℓDR,zl,Θm=NR,zTM˜wrℓLz,Θm.Using the same strategy, the damping matrix, stiffness matrix, and load vector over the extensive displacement vector, i.e.,Cwr, Kwr, and Fwr, can be derived in the same way as Mwr.


## 2.1. Displacement Vector

For a wheel-rail contact pair, as shown in Figure2, the displacement vector, X˜wr,k, with order of T˜wr×1, can be expressed by(1)X˜wr,k=xw,kyw,kzw,kψw,kLv,kRv,kLz,kRz,kLh,kRh,kLy,kRy,kLn,kRn,kLx,kRx,kT,where the superscript “T” denotes the transposition of vectors or matrices; the subscript “k” denotes the kth wheel-rail contact element; xw, yw, zw, and ψw are the longitudinal, lateral, vertical, and yaw motions of the wheel; by assuming the track irregularity as virtual displacement of the rail at the wheel-rail contact positions, Lv and Rv are the vertical irregularities of the rail at the left and right side, respectively, Lh and Rh are the lateral irregularities of the rail at the left and right side, respectively; moreover, Lz and Rz, respectively, denote the vertical displacement of the left- and right-side rail at the wheel-rail contact points; in the same manner, Ly and Ry denote the lateral displacement of the left- and right-side rail at the wheel-rail contact points, respectively; Ln and Rn denote the angle of torsion of the left- and right-side rail around the Y axis, respectively; and Lx and Rx denote the longitudinal displacement of the rails with respect to the left and right wheel-rail contact points, respectively.Figure 2
3D random wheel/rail contact model: (a) front view; (b) elliptical contact spot at location of C.
(a)
(b)In normal wheel-rail contacts, it is assumed that the vertical motion of the wheel is nonindependent and being totally determined by the rail vertical displacement, namely,(2)zw,k=12Lv,k+Rv,k+Lz,k+Rz,k,with(3)Lz,k=DR,zl,kNR,z,kT,Rz,k=DR,zr,kNR,z,kT,DR,zl,k=WR1,lθR1,lyWR2,lθR2,ly,DR,zr,k=WR1,rθR1,ryWR2,rθR2,ry,Nr,1=1−3ζw,k/lwr2+2ζw,k/lwr3,Nr,2=ζw,k1−2ζw,k/lwr+ζw,k/lwr2,Nr,3=3ζw,k/lwr2−2ζw,k/lwr3,Nr,4=ζw,kζw,k/lwr2−ζw,k/lwr,NR,z,k=Nr,1Nr,2Nr,3Nr,4,where DR,zl and DR,zr are the nodal displacement vector used to calculate the vertical displacement of the rail beam at the left and right side along the X axis; W and θy are the nodal vertical displacement and angular displacement of the node around the Y axis, in which the subscripts “l” and “r” denote the left- and right-side rails, respectively and the subscripts “R1” and “R2” denote the left and right nodes with respect to the rail beam element, respectively; ζw denotes the distance between the left node of the rail element and the wheel-rail contact position; and NR,z is the Hermite polynomial interpolation used for Euler–Bernoulli beam.When the wheels lose contact with the rails, the vertical motion of the wheel is independent, by assumingzw in Equation (1) as an independent degree of freedom (DOF). In the numerical simulation, the inequalities below are applied to judge whether or not the wheelsets are contacting to the rails [35] (set the left wheel-rail contact side as an example):(4)Lv+Lz−zwminLv+Lz,zw<ξ1,1−ξ2≤Lv+Lzzw≤1,whenzw>0,Lv+Lz>0,1−ξ3≤zwLv+Lz≤1,whenzw<0,Lv+Lz<0,maxLv+Lz,zw≤ξ4,whenzw<0,Lv+Lz>0,where ξi, i=1,2,3,4, is a small positive number depending on the analytical precision.


## 2.2. Wheel-Rail Coupling Matrices

In most of the running conditions, the wheels constantly keep contact with the rails, thus the wheel-rail coupling matrices that consider the wheelset’s vertical displacement as nonindependent will be presented here as an emphasis. The wheel-rail coupling matrices are marked by subscript “wr” with order ofT˜wr×T˜wr.The mass submatrix of WRI element can be written as(5)M˜wr=∑k=1nwrMwr,kΘm,k,Θm,k+Mwr,kℓxw,k,ℓxw,k+Mwr,kℓyw,k,ℓyw,k+Mwr,kℓψw,k,ℓψw,k,=∑k=1nwrmw,z,kNmTNm+mw,y,k+mw,x,k+Iwz,k,with Nm=1/21/21/21/2, where mw,z, mw,y, and mw,x denote the mass of the wheelset against the vertical, lateral, and longitudinal motion of the wheelset; Iwz is the mass moment of inertia of the wheelset around the Z axis; Θm is the number of the freedoms in Equation (1) related to the mass submatrix, Θm=ℓLvRvLzRz, ℓ⋅ denotes the operator characterizing the corresponding number of DOF in the displacement vector, and nwr=4 for a four-axle vehicle.The damping submatrix of WRI element is mainly induced by the virtual work of wheel-rail creep forces, which can be derived as(6)C˜wr=∑k=1nwrCwr,1l,kΘc,1l,k,Θc,1l,k+Cwr,1r,kΘc,1r,k,Θc,1r,k+Cwr,2l,kΘc,2l,k,Θc,2l,k+Cwr,2r,kΘc,2r,k,Θc,2r,k=∑k=1nwrCXV¯Nc1,lTNc1,l+Nc1,rTNc1,r+CYV¯Nc2,lTNc2,l+Nc2,rTNc2,r,with(7)Nc1,l=1llr−1,Nc1,r=1−llr−1,Nc2,l=1−1−1lcp,Nc2,r=1−1−1lcp,Θc,1l,k=ℓxw,kψw,kLx,k,Θc,1r,k=ℓxw,kψw,kRx,k,Θc,2l,k=ℓyw,kLh,kLy,kLn,k,Θc,2r,k=ℓyw,kRh,kRy,kRn,k,where CX and CY denote the longitudinal and lateral creep coefficient, respectively; V¯ is the forward speed of the wheelset nominally; Cwr,1l and Cwr,1r denote the damping matrices derived by the virtual work of wheel-rail longitudinal creep forces; Cwr,2l and Cwr,2r denote the damping matrices derived by the virtual work of wheel-rail lateral creep forces; llr is the lateral half distance the two wheel-rail contact points at the left and right side, and lcp is the vertical distance between the wheel-rail contact point and the centroid of the wheelset.The stiffness submatrix of WRI element is mainly derived from the work of wheel-rail equivalent gravity stiffness, namely,(8)K˜wr=∑k=1nwrKwr,kΘs,k,Θs,k=∑n=1nwrkg,kNs,kTNs,k,with(9)Θs,k=ℓyw,kLh,kRh,kLy,kRy,k,Ns,k=ℓ1−1/2−1/2−1/2−1/2,where kg,k is the wheel-rail equivalent gravity stiffness, which is calculated using the following equation [36]:(10)kg=Wbϕςδ0+ε,with(11)ε=bR−R′ϕ11+Rδ0b,ς′=RR−R′ϕ11+R′δ0b,ϕ1=1−r0δ0b,where ε is the angular parameter of the wheel-rail contact, ς is the change rate of the lateral relative displacement between the wheel centroid and wheel-rail contact point, b is the half of lateral distance between two rolling circles, R and R′ are the radii of curvature of the wheel and the rail at the contact point, respectively, δ0 is the wheel-rail contact angle when the wheelset is placed in the central position, and r0 is the radius of the nominal rolling circle.The load vector of WRI element is originated from the gravity of the wheelset, namely,(12)F˜wr=∑k=1nwrFwr,kΘfl,k,Θfl,k+Fwr,kΘfr,k,Θfr,k=∑k=1nwr12mw,kgNfl+12mw,kgNfr,with (13)Θfl,k=ℓLv,kLz,k,Θfr,k=ℓRv,kRz,k,Nfl=Nfr=11,where g is the gravitational acceleration.Equations (5), (6), and (8) have briefly presented the mass, damping, and stiffness matrices of WRI element, namely, M˜wr, C˜wr, and K˜wr all with order of nwrT˜wr×nwrT˜wr. Moreover, it should be noted that the some of the DOFs in Equation (1), i.e., Lz,k, Rz,k, Ly,k, Ry,k, Ln,k, Rn,k, Lx,k, and Rx,k, which can be called as the interior node, are the interpolation results of the nodal displacements of the rail beam element. Besides Lz,k=DR,zlNR,zT and Rz,k=DR,zrNR,zT illustrated in Equation (2), the other DOFs can be equivalently transformed by the nodal displacement as(14)Ly,k=DR,ylNR,yT,Ry,k=DR,yrNR,yT,with(15)DR,yl=VR1,lθR1,lzVR2,lθR2,lz,DR,yr=VR1,rθR1,rzVR2,rθR2,rz,NR,y=Nr,1−Nr,2Nr,3−Nr,4,where V denotes the lateral displacement of the rail node, θz denotes the angular displacement of the rail node around the Z axis.Moreover,(16)Ln,k=DR,nlNR,nT,Rn,k=DR,nrNR,nT,with(17)DR,nl=θR1,lxθR2,lx,DR,nr=θR1,rxθR2,rx,NR,n=1−ζw,k/lwrζw,k/lwr,where θx denotes the angular displacement of the rail node around the X axis:(18)Lx,k=DR,xlNR,xT,Rx,k=DR,xrNR,xT,with(19)DR,xl=UR1,lxUR2,lx,DR,xr=UR1,rxUR2,rx,NR,x=NR,n,where Ux denotes the longitudinal displacement of the rail node.Based on Equations (2), (14), (16), and (18), the displacement vector, X˜wr,k, presented in Equation (1) can be further extended as(20)X˜wr,k=xw,kyw,kzw,kψw,kTLv,kRv,kWR1,l,kθR1,ly,kWR2,l,kθR2,ly,kWR1,r,kθR1,ry,kWR2,r,kθR2,ry,kTLh,kRh,kVR1,l,kθR1,lz,kVR2,l,kθR2,lz,kVR1,r,kθR1,rz,kVR2,r,kθR2,rz,kTθR1,lx,kθR2,lx,kθR1,rx,kθR2,rx,kTUR1,lx,kUR2,lx,kUR1,rx,kUR2,rx,kT.In the meantime, the dynamic matrices and load vector over WRI element can be conveniently redistributed into this extensive nodal displacement vector; for example, the newly developed mass matrixMwr with order Twr×Twr can be obtained by(21)Mwr=M˜wrΘm,ℓDR,zr+M˜wrΘm,ℓDR,zl+M˜wrℓDR,zr,Θm+M˜wrℓDR,zl,Θm,with(22)M˜wrΘm,ℓDR,zr=M˜wrΘm,ℓRzNR,z,M˜wrΘm,ℓDR,zl=M˜wrΘm,ℓLzNR,z,M˜wrℓDR,zr,Θm=NR,zTM˜wrℓRz,Θm,M˜wrℓDR,zl,Θm=NR,zTM˜wrℓLz,Θm.Using the same strategy, the damping matrix, stiffness matrix, and load vector over the extensive displacement vector, i.e.,Cwr, Kwr, and Fwr, can be derived in the same way as Mwr.


## 3. Vehicle-Track Interaction Model

As illustrated in Figure1, a ballastless track system was formulated with a moving vehicle subjected to constant velocity. The railway vehicle was modelled as a four-axle mass-spring-damper system, which consisted of a car body, two bogie frames, four wheelsets, and two stage suspensions. The car body and bogie frames have 6 DOFs each including motions of longitudinal X, transverse Y, bounce Z, roll Φ, yaw ψ, and pitch β with respect to their mass centroids, respectively, while for the wheelsets, only the motions of longitudinal, transverse, vertical, and yaw were considered. The total number of DOFs of vehicle system was therefore 34. The vehicle proceeded with velocity v in the longitudinal direction.In the track system, the rails were modelled as spatial Bernoulli–Euler beam supported by discrete viscoelastic supported with finite length, and the rail pads were assumed as linear spring-damper elements. Moreover, slab tracks were assumed to be thin plate elements connected to subgrade through cement asphalt mortar (CA mortar) layer, which was regarded as continuous plane springs and dashpots.


### 3.1. Equation of Motion for the Vehicle-Track Interaction Systems

The wheel-rail coupling system had been formulated and characterized by dynamic matrices in Section2, based on which one can derive the equation of motion written in submatrix form for the vehicle-track interaction systems by energy principle [13]:(23)Mvv000Mwr000MssX¨vX¨wrX¨s+CvvCv−wr0Cwr−vCwrCwr−s0Cs−wrCssX˙vX˙wrX˙s+KvvKv−wr0Kwr−vKwrKwr−s0Ks−wrKssXvXwrXs=FvFwrFs,where M, C, and K denote the mass, damping, and stiffness submatrices, respectively; X and F denote the displacement and force vectors, respectively; the subscripts “v,” “wr,” and “s” denote the vehicle, wheel-rail interaction element, and slab track, respectively.The wheel-rail coupling matrices, namely,Mwr, Cwr, and Kwr, were core terms in Equation (23). In Equation (23), the elements in the wheel-rail coupling matrices can be partitioned and inserted into different positions based on the “set-in-right-position” rule [37]. The submatrices of the slab track accompanied by their interactions with the rail elements can be found in [38], here not given for brevity.The submatrices of vehicle systems, the coupling submatrices between vehicle systems and wheel-rail interaction element and the force vectors will be elaborated in the following sections.


### 3.2. The “Set-in-Right-Position” Rule

It is assumed that there are totallyn displacement parameters γi, i=1,2,...,n, and the total elastic potential energy is ∏γi. Based on the total potential energy with stationary value, the first-order variation of ∏ is equal to zero, namely,(24)δ∏=∑i=1n∂∏∂γiδγi=0,and accordingly, there is(25)δγi∂∏∂γi=0,i=1,2,...,n.It is obvious that Equation (25) denotes the i-th equilibrium equation of this dynamic system. Herein, δγi denotes the i-th row in the matrix; besides, there are displacement parameters γj, j=1,2,...,n derived from ∂∏/∂γi. In the dynamic matrices, the serial number j denotes the j-th column. Hence, the coefficient multiplied by δγi and γj should be set in the cross position of the matrix at the i-th row and the j-th column.


### 3.3. Matrices for the Vehicle Systems

The matrices of the vehicle were marked with the subscript ‘vv’. The mass matrix of the vehicle can be written as(26)Mvv=McΘc,Θc+MGqΘGq,ΘGq+MGhΘGh,ΘGh+∑k=1nwrMw,kΘw,k,Θw,k,with(27)Θq=ℓxqyqzqψqβqϕq,MqΘq,Θq=diagmqmqmqIq,zIq,yIq,x,Θw,k=ℓxw,kyw,kψw,k,Mw,kΘw,k,Θw,k=diagmw,kmw,kIw,z,k,where x, y, z, ψ, β, and ϕ denote the longitudinal, lateral, vertical, yaw, pitch, and roll motions, respectively; the subscript “q” will be substituted by “c,” “Gq,” and “Gh” denoting the car body, front bogie frame, rear bogie frame, respectively; m denotes the mass; Ix, Iy, and Iz denote the moment of inertia of rigid-bodies around the X, Y, and Z axes, respectively.The stiffness matrix of the vehicle can be written as(28)Kvv=Ksz1Θsz1,Θsz1+Ksz2Θsz2,Θsz2+Ksy1Θsy1,Θsy1+Ksy2Θsy2,Θsy2+Ksx1Θsx1,Θsx1+Ksx2Θsx2,Θsx2,where Ksz1, Ksy1, and Ksx1 denote the stiffness matrices induced by interactions between the car body and the front bogie frame in the vertical, lateral, and longitudinal directions, respectively; Ksz2, Ksy2, and Ksx2 denote the stiffness matrices induced by interactions between the car body and the rear bogie frame in the vertical, lateral, and longitudinal directions, respectively; Nsz, Nsy, and Nsx denote the equivalent shape functions used to characterize the geometrical relationships of the involved DOFs. The detail expressions of the parameters in Equation (28) are illustrated in Appendix A.The damping matrix, i.e.,Cvv, with the same order of Kvv, can be obtained by substituting the stiffness coefficient k in Kvv with damping coefficient c.


### 3.4. Matrices for the Vehicle-WRI Element Interactions

The coupling matrices for vehicle-WRI element interactions were marked with the subscript ‘v-wr’ or ‘wr-v’. Before deriving the matricesKv−wr, Kwr−v, Cv−wr, and Cwr−v, one can deduce the interaction matrices between the vehicle and the rails at first, namely, Kvr and Cvr.K
vr can be written as(29)Kvr=Kvr,z+Kvr,y+Kvr,x,where the subscripts ‘z,’ ‘y,’ and ‘x’ denote the interactions between the bogies and the i-th wheelset in the vertical, lateral, and longitudinal direction, respectively. The details about the formation of these involved stiffness matrices have been presented in Appendix B.The damping matrix, i.e.,Cvr, with the same order of Kvr, can be obtained by substituting the stiffness coefficient k in Kvr with damping coefficient c.The stiffness matricesKv−wr and Kwr−v can be obtained by partitioning Kvr, that is,(30)Kv−wr=KvrΘv,Θwr,Kwr−v=KvrΘwr,Θv,in which(31)Θv=ℓxGpyGpzGpψGpβGpϕGpxGhyGhzGhψGhβGhϕGh,Θwr=ℓxw,iyw,iψw,iLv,iRv,iLz,iRz,i.The damping matricesCv−wr and Cwr−v had the same expression as Kv−wr and Kwr−v, respectively, by merely substituting the stiffness coefficient k with damping coefficient c.


### 3.5. Force Vectors

The force vector of the vehicle can be written as(32)Fv=FvΘv,Θv,with(33)Θv=xcyczcψcβcϕcxGqyGqzGqψGqβGqϕGqxGhyGhzGhψGhβGhϕGh,FvΘv,Θv=mcmcmcIczIcyIcxmGqmGqmGqIGqzIGqyIGqxmGhmGhmGhIGhzIGhyIGhx.Moreover, it can be inferred from Section2 that the WRI element can be regarded as a self-excited system, in which the dynamic control equations can be unified as a linearized homogeneous equation without consideration of the gravitational forces, that is,(34)MwrX¨wr+CwrX˙wr+KwrXwr=0.In Equation (34), there exist DOFs with already known displacement, i.e., Lv, Rv, Lh, and Rh, and thus Equation (34) can be partitioned as(35)MUUMUNMNUMNNX¨UX¨N+CUUCUNCNUCNNX˙UX˙N+KUUKUNKNUKNNXUXN=00,in which the subscript “U” denotes the unknown displacement quantities and “N” denotes the known displacement quantities.Then, Equation (35) can be extended as(36)MUUX¨U+CUUX˙U+KUUXU=−MUNX¨N−CUNX˙N−KUNXN,MNNX¨N+CNNX˙N+KNNXN=−MNUX¨U−CNUX˙U−KNUXU.In the above Equation, the corresponding rows and columns of the known displacement quantities in the matrix equation should be deleted. Obviously, the force vectorFwr can be obtained by(37)Fwr=F˜wr−MwrΘk,ΘnX¨n−CwrΘk,ΘnX˙n−KwrΘk,ΘnXn.All elements in the load vector of the slab trackFs are zero.Till now, all details about the construction of this linear vehicle-track dynamic model had been illustrated. Since the shape functions, i.e.,NR,z in Equation (2) and the differentiation with respect to the local coordinate ξw are time-dependent, the dynamic matrices and force vectors in Equation (23) should be updated in each time-step. Moreover, unlike the most vehicle-track models presented in Refs. [3, 39], no iteration procedures are required to guarantee the response stability of the dynamic model.Equation (23) can be solved by the step-by-step integration method to obtain simultaneously the dynamic responses of vehicle-track systems even using relatively large integration-steps.


## 3.1. Equation of Motion for the Vehicle-Track Interaction Systems

The wheel-rail coupling system had been formulated and characterized by dynamic matrices in Section2, based on which one can derive the equation of motion written in submatrix form for the vehicle-track interaction systems by energy principle [13]:(23)Mvv000Mwr000MssX¨vX¨wrX¨s+CvvCv−wr0Cwr−vCwrCwr−s0Cs−wrCssX˙vX˙wrX˙s+KvvKv−wr0Kwr−vKwrKwr−s0Ks−wrKssXvXwrXs=FvFwrFs,where M, C, and K denote the mass, damping, and stiffness submatrices, respectively; X and F denote the displacement and force vectors, respectively; the subscripts “v,” “wr,” and “s” denote the vehicle, wheel-rail interaction element, and slab track, respectively.The wheel-rail coupling matrices, namely,Mwr, Cwr, and Kwr, were core terms in Equation (23). In Equation (23), the elements in the wheel-rail coupling matrices can be partitioned and inserted into different positions based on the “set-in-right-position” rule [37]. The submatrices of the slab track accompanied by their interactions with the rail elements can be found in [38], here not given for brevity.The submatrices of vehicle systems, the coupling submatrices between vehicle systems and wheel-rail interaction element and the force vectors will be elaborated in the following sections.


## 3.2. The “Set-in-Right-Position” Rule

It is assumed that there are totallyn displacement parameters γi, i=1,2,...,n, and the total elastic potential energy is ∏γi. Based on the total potential energy with stationary value, the first-order variation of ∏ is equal to zero, namely,(24)δ∏=∑i=1n∂∏∂γiδγi=0,and accordingly, there is(25)δγi∂∏∂γi=0,i=1,2,...,n.It is obvious that Equation (25) denotes the i-th equilibrium equation of this dynamic system. Herein, δγi denotes the i-th row in the matrix; besides, there are displacement parameters γj, j=1,2,...,n derived from ∂∏/∂γi. In the dynamic matrices, the serial number j denotes the j-th column. Hence, the coefficient multiplied by δγi and γj should be set in the cross position of the matrix at the i-th row and the j-th column.


## 3.3. Matrices for the Vehicle Systems

The matrices of the vehicle were marked with the subscript ‘vv’. The mass matrix of the vehicle can be written as(26)Mvv=McΘc,Θc+MGqΘGq,ΘGq+MGhΘGh,ΘGh+∑k=1nwrMw,kΘw,k,Θw,k,with(27)Θq=ℓxqyqzqψqβqϕq,MqΘq,Θq=diagmqmqmqIq,zIq,yIq,x,Θw,k=ℓxw,kyw,kψw,k,Mw,kΘw,k,Θw,k=diagmw,kmw,kIw,z,k,where x, y, z, ψ, β, and ϕ denote the longitudinal, lateral, vertical, yaw, pitch, and roll motions, respectively; the subscript “q” will be substituted by “c,” “Gq,” and “Gh” denoting the car body, front bogie frame, rear bogie frame, respectively; m denotes the mass; Ix, Iy, and Iz denote the moment of inertia of rigid-bodies around the X, Y, and Z axes, respectively.The stiffness matrix of the vehicle can be written as(28)Kvv=Ksz1Θsz1,Θsz1+Ksz2Θsz2,Θsz2+Ksy1Θsy1,Θsy1+Ksy2Θsy2,Θsy2+Ksx1Θsx1,Θsx1+Ksx2Θsx2,Θsx2,where Ksz1, Ksy1, and Ksx1 denote the stiffness matrices induced by interactions between the car body and the front bogie frame in the vertical, lateral, and longitudinal directions, respectively; Ksz2, Ksy2, and Ksx2 denote the stiffness matrices induced by interactions between the car body and the rear bogie frame in the vertical, lateral, and longitudinal directions, respectively; Nsz, Nsy, and Nsx denote the equivalent shape functions used to characterize the geometrical relationships of the involved DOFs. The detail expressions of the parameters in Equation (28) are illustrated in Appendix A.The damping matrix, i.e.,Cvv, with the same order of Kvv, can be obtained by substituting the stiffness coefficient k in Kvv with damping coefficient c.


## 3.4. Matrices for the Vehicle-WRI Element Interactions

The coupling matrices for vehicle-WRI element interactions were marked with the subscript ‘v-wr’ or ‘wr-v’. Before deriving the matricesKv−wr, Kwr−v, Cv−wr, and Cwr−v, one can deduce the interaction matrices between the vehicle and the rails at first, namely, Kvr and Cvr.K
vr can be written as(29)Kvr=Kvr,z+Kvr,y+Kvr,x,where the subscripts ‘z,’ ‘y,’ and ‘x’ denote the interactions between the bogies and the i-th wheelset in the vertical, lateral, and longitudinal direction, respectively. The details about the formation of these involved stiffness matrices have been presented in Appendix B.The damping matrix, i.e.,Cvr, with the same order of Kvr, can be obtained by substituting the stiffness coefficient k in Kvr with damping coefficient c.The stiffness matricesKv−wr and Kwr−v can be obtained by partitioning Kvr, that is,(30)Kv−wr=KvrΘv,Θwr,Kwr−v=KvrΘwr,Θv,in which(31)Θv=ℓxGpyGpzGpψGpβGpϕGpxGhyGhzGhψGhβGhϕGh,Θwr=ℓxw,iyw,iψw,iLv,iRv,iLz,iRz,i.The damping matricesCv−wr and Cwr−v had the same expression as Kv−wr and Kwr−v, respectively, by merely substituting the stiffness coefficient k with damping coefficient c.


## 3.5. Force Vectors

The force vector of the vehicle can be written as(32)Fv=FvΘv,Θv,with(33)Θv=xcyczcψcβcϕcxGqyGqzGqψGqβGqϕGqxGhyGhzGhψGhβGhϕGh,FvΘv,Θv=mcmcmcIczIcyIcxmGqmGqmGqIGqzIGqyIGqxmGhmGhmGhIGhzIGhyIGhx.Moreover, it can be inferred from Section2 that the WRI element can be regarded as a self-excited system, in which the dynamic control equations can be unified as a linearized homogeneous equation without consideration of the gravitational forces, that is,(34)MwrX¨wr+CwrX˙wr+KwrXwr=0.In Equation (34), there exist DOFs with already known displacement, i.e., Lv, Rv, Lh, and Rh, and thus Equation (34) can be partitioned as(35)MUUMUNMNUMNNX¨UX¨N+CUUCUNCNUCNNX˙UX˙N+KUUKUNKNUKNNXUXN=00,in which the subscript “U” denotes the unknown displacement quantities and “N” denotes the known displacement quantities.Then, Equation (35) can be extended as(36)MUUX¨U+CUUX˙U+KUUXU=−MUNX¨N−CUNX˙N−KUNXN,MNNX¨N+CNNX˙N+KNNXN=−MNUX¨U−CNUX˙U−KNUXU.In the above Equation, the corresponding rows and columns of the known displacement quantities in the matrix equation should be deleted. Obviously, the force vectorFwr can be obtained by(37)Fwr=F˜wr−MwrΘk,ΘnX¨n−CwrΘk,ΘnX˙n−KwrΘk,ΘnXn.All elements in the load vector of the slab trackFs are zero.Till now, all details about the construction of this linear vehicle-track dynamic model had been illustrated. Since the shape functions, i.e.,NR,z in Equation (2) and the differentiation with respect to the local coordinate ξw are time-dependent, the dynamic matrices and force vectors in Equation (23) should be updated in each time-step. Moreover, unlike the most vehicle-track models presented in Refs. [3, 39], no iteration procedures are required to guarantee the response stability of the dynamic model.Equation (23) can be solved by the step-by-step integration method to obtain simultaneously the dynamic responses of vehicle-track systems even using relatively large integration-steps.


## 4. Numerical Examples

In the numerical examples, it is assumed that a railway vehicle runs with a constant velocity of 350 km/h on the ballastless slab-tracks. For constructing the vehicle-track model, a track segment element is taken between two adjacent rail pads [40]; moreover, a cyclic calculation model proposed by Xu et al. [41] is applied to increase the computational efficiency lowered by the high DOFs of the long track finite element model, and the boundary effect of the tracks can be reduced fairly well.The physical parameters of the vehicle and the tracks can be referred in Tables1 and 2. The track irregularity, as the excitation of the wheel-rail interactions, is obtained by a track irregularity probabilistic model [42] where the measured track irregularities from a high-speed line of China are used as the data sources.Table 1
Main parameters of railway vehicles used in the simulation.
Notation
Parameter
Values
M
c
Car body mass (kg)
34500
M
t
Bogie mass (kg)
3500
M
w
Wheelset mass (kg)
2200
I
cx
Mass moment of inertia of car body about theX axis (kg m2)
124304
I
cy
Mass moment of inertia of car body about theY axis (kg m2)
1560400
I
cz
Mass moment of inertia of car body about theZ axis (kg m2)
1521000
I
tx
Mass moment of inertia of bogie about theX axis (kg m2)
2522
I
ty
Mass moment of inertia of bogie about theY axis (kg m2)
1722
I
tz
Mass moment of inertia of bogie about theZ axis (kg m2)
3100
I
wx
Mass moment of inertia of wheelset about theX axis (kg m2)
730
I
wy
Mass moment of inertia of wheelset about theY axis (kg m2)
84
I
wz
Mass moment of inertia of wheelset about theZ axis (kg m2)
990
K
px
Stiffness coefficient of primary suspension along theX axis (MN/m)
11.68
K
py
Stiffness coefficient of primary suspension along theY axis (MN/m)
7.40
K
pz
Stiffness coefficient of primary suspension along theZ axis (MN/m)
1.173
K
sx
Stiffness coefficient of secondary suspension along theX axis (MN/m)
0.1790
K
sy
Stiffness coefficient of secondary suspension along theY axis (MN/m)
0.1790
K
sz
Stiffness coefficient of secondary suspension along theZ axis (MN/m)
0.61
C
pz
Damping coefficient of primary suspension along theZ axis (kN s/m)
19.60
C
sy
Damping coefficient of secondary suspension along theY axis (kN s/m)
54.80
C
sz
Damping coefficient of secondary suspension along theZ axis (kN s/m)
9.80
L
c
Semilongitudinal distance between bogies (m)
8.75
L
t
Semilongitudinal distance between wheelsets in bogie (m)
1.25
R
0
Wheel radius (m)
0.43Table 2
Main parameters of railway tracks used in the simulation.
Notation
Parameter
Values
E
Elastic modulus of rail (N/m2)
1970400
ρ
Density of rail (kg/m3)
7.16 × 103
I
0
Torsional inertia of rail (m4)
3.731 × 10−5
I
y
Rail second moment of area about theY axis (m4)
3.257 × 10−5
I
z
Rail second moment of area about theZ axis (m4)
5.20 × 10−6
GK
Rail torsional stiffness (N m/rad)
1.9587 × 105
K
pv
Fastener stiffness in vertical direction (N/m)
2.9 × 107
K
ph
Fastener stiffness in lateral direction (N/m)
1.5 × 107
C
hp
Fastener damping in vertical direction (N s/m)
5 × 104
C
h
Fastener damping in lateral direction (N s/m)
5 × 104
E
s
Elastic modulus of rail (N/m2)
3.6 × 1011
μ
Poisson ratio of slab
0.25
M
s
Mass per unit volume of slab (kg/m3)
2400
H
s
Height of slab (m)
0.19
L
s
Length of slab (m)
4.95
W
s
Width of slab (m)
2.40
K
vCA
Equivalent vertical stiffness coefficient of CA layer (N/m)
6.0 × 107
C
vCA
Equivalent vertical damping coefficient of CA layer (N/m)
9.0 × 104
K
hCA
Equivalent vertical stiffness coefficient of CA layer (N/m)
1.03 × 106
C
hCA
Equivalent vertical damping coefficient of CA layer (N/m)
1.0 × 106


### 4.1. Comparison between the Linear Model and the Nonlinear Model

To investigate the efficiency and feasibility of this presented model, a comparative analysis is conducted against the nonlinear model developed by Xu and Zhai [11]. In the nonlinear model, the three-dimensional wheel-rail nonlinear contacts are considered comprehensively.Figures3 and 4 show the comparisons on the lateral and vertical accelerations of the car body. As seen from these two figures, though there exists unavoidable discrepancy between the models on response amplitudes, the tendency of the time-domain responses of the linear model coincides well with that of the nonlinear model; moreover, the dominant frequencies of the car body accelerations drawn from these two models are almost the same. For the vertical acceleration of the car body, the absolute maximum values are 1.6 m/s2 and 1.41 m/s2; with respect to the nonlinear model and linear model, the dominant frequencies are 1.17 Hz and 6.25 Hz. While for the lateral acceleration, the absolute maximum values of the nonlinear and linear model are 0.145 m/s2 and 0.144 m/s2, respectively; besides, the domain frequencies for these two models are 1.72 Hz, 5.46 Hz, and 10.94 Hz and are similar, but there exists a slight difference in spectral densities inevitably.Figure 3
Comparison of the vertical acceleration of the car body between the linear model and the nonlinear model in the (a) time domain and (b) frequency domain.
(a)
(b)Figure 4
Comparison of the lateral acceleration of the car body between the linear model and the nonlinear model in the (a) time domain and (b) frequency domain.
(a)
(b)Instead of calculating the wheel-rail force using the nonlinear Hertzian contact theory [43]:(38)Fv=Kδ3/2,where K is the Hertz contact coefficient and δ is the wheel-rail relative elastic compression, the wheel-rail vertical force Fv in this paper is obtained by introducing D’Alembert’s principle:(39)∑iFi−miaiΔri=0,where i is an integer indicating a variable; Fi is the total applied force on the i-th particle; mi is the mass of the i-th particle; ai is the acceleration, and Δri is the virtual displacement of the i-th particle; Fvconsists of the rolling force Froll, vertical inertial force Finer, and gravitational force Fgrav, namely, Fv=Froll+Finer+Fgrav.As an example, the forces acting on the left side of the wheelsets can be written, respectively, as(40)Froll,i=mcy¨c4−Iczψ¨c4lcHcb+Hbt+Htw+Rw2b0−Icxϕc8bo+mGqy¨Gq4b0Htw+Rw−IGqxϕGq4b0−−1i+1IGqzψ¨Gq2ltHtw+Rw2b0,Froll,i=mcz¨c4−Icyβ¨c4lc+mGqz¨Gq2−−1i+1IGq,yβ¨Gq2lt+mw,iz¨w,i,Fgrav,i=mwg2,where the subscript “i”, i=1,2,3,4, denotes the ith wheelset; the subscript “Gq” represents the front bogie frame, which is only applicable for i = 1 and 2; when i = 3 or 4, the parameters “mGq,” “y¨Gq,” “z¨Gq,” “IGq,x,” “ϕ¨Gq,” “IGq,z,” “ψ¨Gq,” “IGq,y,” and “β¨Gq” should be substituted by “mGh,” “y¨Gh,” “z¨Gh,” “IGh,x,” “ϕ¨Gh,” “IGh,z,” “ψ¨Gh,” “IGh,y,” and “β¨Gh,” respectively, in which the subscript “Gh” represents the rear bogie frame; b0 denotes the half of the lateral distance between the two wheel-rail contact point at the initial condition; Rwis the nominal rolling radius of the wheel; and Htw is the vertical distance between the centroid of bogie frame and the central line of the wheelset.Figure5 illustrates the discrepancy of the wheel-rail vertical force between the nonlinear model and the linear model due to the implementation of different wheel-rail coupling models. As observed from Figure 4, there is a significant correlation between the results of these two models; see for instance, the sections possessing violent wheel-rail interactions are almost the same, but the forces derived by the nonlinear model are locally larger than those of the linear model, especially at high frequency ranges shown in Figure 5(b), because the time interval of integration of the nonlinear model should be greatly smaller than the linear model for guaranteeing the stability of time integration. Moreover, Figure 6 shows the dynamic responses of the rail lateral displacement at both time and frequency domains, from which one can observe that the results of these two models coincide well with each other. The linear model has accurately characterized the dominant frequencies of rail lateral displacement in the whole frequency range of 1–200 Hz though there are differences in spectral values.Figure 5
Comparison of the wheel-rail vertical force between the linear model and the nonlinear model in the (a) time domain and (b) frequency domain.
(a)
(b)Figure 6
Comparison of rail lateral displacement between the linear model and the nonlinear model in the (a) time domain and (b) frequency domain.
(a)
(b)In the numerical calculation, if the railway vehicle moves on the track for one second, the computation time for the linear model is 15 sec, while for the nonlinear model, it is 156 sec in MATLAB®, and thus the computational efficiency has been significantly improved.


### 4.2. Random Vibrations concerning Randomness of Creep Coefficients

Creep coefficient possesses key influence on wheel-rail tangential interactions. Through numerical investigations, it is noticed that the lateral creep coefficientCY in the linearized model will majorly determine the vehicle-track lateral dynamics under irregular excitations. Undoubtedly, the creep coefficients are time-dependent and show randomicity. For obtaining more complete results, there is a necessity to treat the model as an ergodic dynamical system where the most influential parameters, e.g., creep coefficient, should be considered with an ergodicity.In this study, in order to ensure the reliability of dynamic assessment of a linear model, further work on random vibration analysis will be probed into, where the detailed procedures can be illustrated as follows:(1)
To reveal the full amplitude and probability characteristics of the creep coefficients, the track irregularity probabilistic model [42] and the nonlinear vehicle-track coupled model [11] are applied to characterize the probability density function (PDF) of CY.As an example, Figure7 plots the PDF fitting results of lateral creep coefficient, from which it is found out that the τ location-scale distribution can be properly used to reflect the probability property of the creep coefficients, where τ location-scale distribution can be expressed by(41)fx=Γυ+1/2υπΓυ/21+x−μ/σ2/υ−υ+1/2,where Γ⋅ is the gamma function and μ, σ, and υ denote the location, scale, and shape parameter, respectively.(2)
Since the PDF ofCY denoted by pCy is obtained in the above step, the methods for random parameters selection presented in Ref. [5] can be implemented to get the represented random field of CY, that is, ΩCy.(3)
SetH· as the operator of vehicle-track interaction model and loading ΩCy and track irregularity vector TIrre into H·, the random vibration of vehicle-track systems Rindi can be calculated out consequently, that is,(42)Rindi=HΩCyTirre,where the subscript “indi” denotes the dynamic indices of vehicle-track systems.Figure 7
PDF fitting for lateral creep coefficient (GEV, generalized extreme value).Figure8 plots the up and down boundary of power spectral density (PSD), denoted by Pu and Pd, respectively, of rail lateral displacement and wheel-rail lateral force, respectively, within which the confidence probability is 99%, that is,(43)Pru−dω=∫PuxPdxfPindix;ωdx=0.99,where Pr· denotes the probability function; ω is the frequency; fPindix;ω is the PDF of spectral densities against specific frequency ω; and x is the spectral density. The red lines are the PSD results of the nonlinear model.Figure 8
The up- and downboundary PSDs of the linear model and the PSD of the nonlinear model. (a) Rail lateral displacement. (b) Wheel-rail lateral force.
(a)
(b)It can be observed from Figure8 that the boundary PSDs derived by the linear model can properly envelop the results of the nonlinear model, based on the mathematical relationship between spectral density and time-domain amplitude, that is,(44)A=1π∫ωi−1ωiPωdω.It is obvious that time-domain responses of the linear model will inevitably envelop those of the nonlinear model. Therefore, with an implementation of the random vibration analysis, the results derived by the linear model can be provided as a reference to guarantee the safety margin of the vehicle-track dynamic behaviours.However, it should be noted that though the linear model is capable of characterizing the characteristic frequencies of dynamic indices, the chosen of the creep coefficient will exert significant influence on the response amplitude; moreover, the responses of wheel-rail lateral force at low frequency of 10 Hz below will increase violently, which should be treated as one of the deficiencies of the linear model with an emphasis. It is therefore a rather important work to choose the creep coefficient optimally.


## 4.1. Comparison between the Linear Model and the Nonlinear Model

To investigate the efficiency and feasibility of this presented model, a comparative analysis is conducted against the nonlinear model developed by Xu and Zhai [11]. In the nonlinear model, the three-dimensional wheel-rail nonlinear contacts are considered comprehensively.Figures3 and 4 show the comparisons on the lateral and vertical accelerations of the car body. As seen from these two figures, though there exists unavoidable discrepancy between the models on response amplitudes, the tendency of the time-domain responses of the linear model coincides well with that of the nonlinear model; moreover, the dominant frequencies of the car body accelerations drawn from these two models are almost the same. For the vertical acceleration of the car body, the absolute maximum values are 1.6 m/s2 and 1.41 m/s2; with respect to the nonlinear model and linear model, the dominant frequencies are 1.17 Hz and 6.25 Hz. While for the lateral acceleration, the absolute maximum values of the nonlinear and linear model are 0.145 m/s2 and 0.144 m/s2, respectively; besides, the domain frequencies for these two models are 1.72 Hz, 5.46 Hz, and 10.94 Hz and are similar, but there exists a slight difference in spectral densities inevitably.Figure 3
Comparison of the vertical acceleration of the car body between the linear model and the nonlinear model in the (a) time domain and (b) frequency domain.
(a)
(b)Figure 4
Comparison of the lateral acceleration of the car body between the linear model and the nonlinear model in the (a) time domain and (b) frequency domain.
(a)
(b)Instead of calculating the wheel-rail force using the nonlinear Hertzian contact theory [43]:(38)Fv=Kδ3/2,where K is the Hertz contact coefficient and δ is the wheel-rail relative elastic compression, the wheel-rail vertical force Fv in this paper is obtained by introducing D’Alembert’s principle:(39)∑iFi−miaiΔri=0,where i is an integer indicating a variable; Fi is the total applied force on the i-th particle; mi is the mass of the i-th particle; ai is the acceleration, and Δri is the virtual displacement of the i-th particle; Fvconsists of the rolling force Froll, vertical inertial force Finer, and gravitational force Fgrav, namely, Fv=Froll+Finer+Fgrav.As an example, the forces acting on the left side of the wheelsets can be written, respectively, as(40)Froll,i=mcy¨c4−Iczψ¨c4lcHcb+Hbt+Htw+Rw2b0−Icxϕc8bo+mGqy¨Gq4b0Htw+Rw−IGqxϕGq4b0−−1i+1IGqzψ¨Gq2ltHtw+Rw2b0,Froll,i=mcz¨c4−Icyβ¨c4lc+mGqz¨Gq2−−1i+1IGq,yβ¨Gq2lt+mw,iz¨w,i,Fgrav,i=mwg2,where the subscript “i”, i=1,2,3,4, denotes the ith wheelset; the subscript “Gq” represents the front bogie frame, which is only applicable for i = 1 and 2; when i = 3 or 4, the parameters “mGq,” “y¨Gq,” “z¨Gq,” “IGq,x,” “ϕ¨Gq,” “IGq,z,” “ψ¨Gq,” “IGq,y,” and “β¨Gq” should be substituted by “mGh,” “y¨Gh,” “z¨Gh,” “IGh,x,” “ϕ¨Gh,” “IGh,z,” “ψ¨Gh,” “IGh,y,” and “β¨Gh,” respectively, in which the subscript “Gh” represents the rear bogie frame; b0 denotes the half of the lateral distance between the two wheel-rail contact point at the initial condition; Rwis the nominal rolling radius of the wheel; and Htw is the vertical distance between the centroid of bogie frame and the central line of the wheelset.Figure5 illustrates the discrepancy of the wheel-rail vertical force between the nonlinear model and the linear model due to the implementation of different wheel-rail coupling models. As observed from Figure 4, there is a significant correlation between the results of these two models; see for instance, the sections possessing violent wheel-rail interactions are almost the same, but the forces derived by the nonlinear model are locally larger than those of the linear model, especially at high frequency ranges shown in Figure 5(b), because the time interval of integration of the nonlinear model should be greatly smaller than the linear model for guaranteeing the stability of time integration. Moreover, Figure 6 shows the dynamic responses of the rail lateral displacement at both time and frequency domains, from which one can observe that the results of these two models coincide well with each other. The linear model has accurately characterized the dominant frequencies of rail lateral displacement in the whole frequency range of 1–200 Hz though there are differences in spectral values.Figure 5
Comparison of the wheel-rail vertical force between the linear model and the nonlinear model in the (a) time domain and (b) frequency domain.
(a)
(b)Figure 6
Comparison of rail lateral displacement between the linear model and the nonlinear model in the (a) time domain and (b) frequency domain.
(a)
(b)In the numerical calculation, if the railway vehicle moves on the track for one second, the computation time for the linear model is 15 sec, while for the nonlinear model, it is 156 sec in MATLAB®, and thus the computational efficiency has been significantly improved.


## 4.2. Random Vibrations concerning Randomness of Creep Coefficients

Creep coefficient possesses key influence on wheel-rail tangential interactions. Through numerical investigations, it is noticed that the lateral creep coefficientCY in the linearized model will majorly determine the vehicle-track lateral dynamics under irregular excitations. Undoubtedly, the creep coefficients are time-dependent and show randomicity. For obtaining more complete results, there is a necessity to treat the model as an ergodic dynamical system where the most influential parameters, e.g., creep coefficient, should be considered with an ergodicity.In this study, in order to ensure the reliability of dynamic assessment of a linear model, further work on random vibration analysis will be probed into, where the detailed procedures can be illustrated as follows:(1)
To reveal the full amplitude and probability characteristics of the creep coefficients, the track irregularity probabilistic model [42] and the nonlinear vehicle-track coupled model [11] are applied to characterize the probability density function (PDF) of CY.As an example, Figure7 plots the PDF fitting results of lateral creep coefficient, from which it is found out that the τ location-scale distribution can be properly used to reflect the probability property of the creep coefficients, where τ location-scale distribution can be expressed by(41)fx=Γυ+1/2υπΓυ/21+x−μ/σ2/υ−υ+1/2,where Γ⋅ is the gamma function and μ, σ, and υ denote the location, scale, and shape parameter, respectively.(2)
Since the PDF ofCY denoted by pCy is obtained in the above step, the methods for random parameters selection presented in Ref. [5] can be implemented to get the represented random field of CY, that is, ΩCy.(3)
SetH· as the operator of vehicle-track interaction model and loading ΩCy and track irregularity vector TIrre into H·, the random vibration of vehicle-track systems Rindi can be calculated out consequently, that is,(42)Rindi=HΩCyTirre,where the subscript “indi” denotes the dynamic indices of vehicle-track systems.Figure 7
PDF fitting for lateral creep coefficient (GEV, generalized extreme value).Figure8 plots the up and down boundary of power spectral density (PSD), denoted by Pu and Pd, respectively, of rail lateral displacement and wheel-rail lateral force, respectively, within which the confidence probability is 99%, that is,(43)Pru−dω=∫PuxPdxfPindix;ωdx=0.99,where Pr· denotes the probability function; ω is the frequency; fPindix;ω is the PDF of spectral densities against specific frequency ω; and x is the spectral density. The red lines are the PSD results of the nonlinear model.Figure 8
The up- and downboundary PSDs of the linear model and the PSD of the nonlinear model. (a) Rail lateral displacement. (b) Wheel-rail lateral force.
(a)
(b)It can be observed from Figure8 that the boundary PSDs derived by the linear model can properly envelop the results of the nonlinear model, based on the mathematical relationship between spectral density and time-domain amplitude, that is,(44)A=1π∫ωi−1ωiPωdω.It is obvious that time-domain responses of the linear model will inevitably envelop those of the nonlinear model. Therefore, with an implementation of the random vibration analysis, the results derived by the linear model can be provided as a reference to guarantee the safety margin of the vehicle-track dynamic behaviours.However, it should be noted that though the linear model is capable of characterizing the characteristic frequencies of dynamic indices, the chosen of the creep coefficient will exert significant influence on the response amplitude; moreover, the responses of wheel-rail lateral force at low frequency of 10 Hz below will increase violently, which should be treated as one of the deficiencies of the linear model with an emphasis. It is therefore a rather important work to choose the creep coefficient optimally.


## 5. Conclusions

In this study, the modelling framework for vehicle-track interactions is proposed with linearized wheel-rail contacts/creepages. Key steps for modelling the vehicle-track interactions lie in the establishment of the coupling matrices between subsystems of a vehicle and the tracks and the WRI elements. The correctness and the higher efficiency of the proposed linear model have been presented in the numerical examples; besides, the following conclusions can be drawn accordingly:(1)
The linearized model can give approachable results to the nonlinear model regarding the response amplitudes and the main dominant frequencies characterized by these two models are also similar. However, there exist obvious differences in frequency spectrum amplitude.(2)
This presented model possesses higher computational stability due to the linearization. However, this linearization might lead to acceptable deviations in engineering practice comparing to the nonlinear model.(3)
In this dynamic condition, the system responses at low frequency range of 0–16.5 Hz will be significantly enlarged by the increase of the lateral creep coefficient.(4)
This current work shows limitations for dynamics problems related to transients, discrete joints, derailment, etc.



---

*Source: 2015465-2019-01-20.xml*