# An Alternative Solution of Train-Track Dynamic Interaction

Authors: Zhaowei Chen; Hui Fang
Journal: Shock and Vibration
 (2019)

---

## Abstract

Train-track interaction (TTI) is a classic research topic in railway engineering, which consists of three main parts, namely, train model, track model, and wheel-rail interaction. To improve the computational accuracy and broaden the application range, an alternative calculation method to investigate TTI based on secondary development technology of the commercial software ANSYS through APDL language is introduced in this article. Primarily, the train-track interaction theory is briefly presented. On this basis, TTI is programmed and implemented on the computing platform of ANSYS by fully taking the nonlinear wheel-rail interaction into consideration. In this calculation method, the train model, which is established based on multibody dynamics theory and solved by an advanced explicit integration method, is programmed into ANSYS through APDL language, while the track part is simulated according to finite element theory. Then, the proposed calculation method is validated with field test results to verify the validity. Finally, a numerical demonstration is conducted employing the present method. Results show that the introduced method is effective and able to investigate TTI. Different complicated track systems can be accurately simulated employing this method. Moreover, this method is also adoptable to explore train-bridge interaction and train-track-bridge interaction.


---

## Body

## 1. Introduction

Train-track interaction (TTI) has been a traditional research topic for quite a long time [1]. In recent years, as the development of high-speed railways and metro systems, a growing attention is paid to this topic. Early research usually did not couple the train system and track system together, and the dynamic responses of the train and the track were investigated separately [2–5]. As a matter of fact, running trains cause track vibrations, which in turn aggravate the train vibrations. Therefore, the train system and track system essentially form an integrated dynamic system by wheel-rail interaction relationship. On this basis, various train-track coupled models have been established and validated in recent years, which are also widely adopted in railway engineering [6–8]. Until now, the train-track interaction theory has already been mature and perfect, and the key issue on TTI is to make the calculation more effective and more accurate.So far, various calculation methods have been proposed to investigate TTI. These methods can be generally classified into the following categories:Method A: dynamic equations of TTI are programmed on different compiled platforms employing different programming languages such as FORTRAN, MATLAB, and C, which are then solved with numerical integration methods. This is the most commonly used method. With this method, the train systems are usually established based on multibody dynamics [9] and the tracks are commonly simplified as Euler beam [10] or Timoshenko beam [11] and solved with the Ritz method [12] or mode superposition method [10]. Adopting this calculation method, the computational efficiency is high and the computational accuracy is enough for engineering [8–10]. However, the disadvantage of this method is that the trains and tracks should be normal and commonly used. While various kinds of trains and tracks are adopted in high-speed railways and metro systems, for each new kind of structure, the dynamic equations and the programs should be rewritten, which costs a lot. Moreover, the structural stress, the nonlinear material property, and the damage behavior of the train-track systems are hardly to be considered through this calculation method.Method B: train systems and track systems are all established in commercial finite element software such as ANSYS and ABAQUS [13, 14]. This calculation method can reach a high computational accuracy; however, in the authors’ opinion, the computational efficiency is extremely low due to the implicit integration algorithm in these finite element software. Another nonnegligible disadvantage is that the track irregularity, which is the most important random excitation for the railway system, is hardly to be accurately considered in this calculation method although the track irregularities can be included either as external time-variant forces or by including the geometrical irregularity of the rail elements in combination with a general contact formulation. Furthermore, the accurate nonlinear wheel-rail contact relationship is also hardly to be simulated.Method C: different commercial software are combined to solve TTI. With this method, train systems are usually modeled by multibody software such as SIMPACK and UM, while the track structures are simulated using finite element software. Then, the modal properties of tracks are exported from finite element software into multibody software to make the two different systems coupled [15, 16]. This method can also reach a high computational efficiency due to the algorithm of multibody dynamics and mode superposition method or modal synthesis method. However, the nonlinear characteristics of the tracks are ignored in the modal properties, which may cause calculation errors for nonlinear structures and systems.Method D: motion equations of trains are programmed with programming languages such as FORTRAN, while the track structures are modeled in finite element software. Then, mass matrix and stiffness matrix of tracks are exported to external files, which are further read in by the self-programed procedure to form the whole dynamic equations of the train-track systems, and finally these equations can be calculated through different integration methods [17]. However, employing this calculation method, the nonlinear material and contact relationship are also hard to be simulated, and the structural stress is different to be acquired.Through the above methods, several generalized conclusions are reached by the authors:(i)
The train systems should be modeled based on multibody dynamics theory to improve the calculation speed(ii)
To conveniently establish various kinds of track structures and to accurately analyze the dynamic behaviors of tracks considering their nonlinear characteristics, finite element software is perfect tools to reach this goal(iii)
Track irregularity is hardly to be considered in finite element software(iv)
Comprehensively considering the computational efficiency and the computational accuracy, the train equations and track systems should be, respectively, solved by the explicit integration method and implicit integration method.Based on the above investigations, to improve the computational accuracy and broaden the application range, an alternative method to solve TTI is proposed in this present work based on development of ANSYS. First, TTI is programmed and implemented on the computing platform of ANSYS in Section2, and then the method is further validated in the next section. After that, the advantage and disadvantage of the proposed calculation method are discussed. Adopting the present calculation method, a numerical case is conducted in Section 5, and some important conclusions are reached in the last section.


## 2. Solution of Train-Track Interaction Based on Secondary Development of ANSYS

Solutions of TTI adopting different technological means have been presented in previous works [6–10, 18]. In this present work, an alternative method for determining TTI based on secondary development of ANSYS with APDL (short for “ANSYS parametric design language”) is clearly introduced.ANSYS is a world-famous finite element software, which contains various kinds of elements and materials. The software has been widely used in railway engineering. Furthermore, APDL broadens the function and application range of ANSYS, resulting in convenient modelling and simulation. Adopting this language, the train system and track irregularity can be easily considered on the computing platform of ANSYS. Hence, the software ANSYS and APDL are adopted to investigate TTI.Taking a 2D vehicle-slab track system as an example, Figure1 illustrates the train-track interaction model consisting of three main parts, namely, the train subsystem, the track subsystem, and the wheel-rail interaction relationship. To further describe the theory of TTI in detail, the adopted notations in the following texts are listed in Table 1.Figure 1
Typical train-track interaction model (slab ballastless track).Table 1
Main notations adopted in this work.
Notation
Physical meaning
Unit
Z
c/Zt/Zw
Vertical displacements of car body, frame, and wheelset
m
β
c/βt/βw
Pitch angles of car body, frame, and wheelset
°
M
c/Mt/Mw
Masses of car body, frame, and wheelset
kg
C
pz/Csz
Damping of primary suspension and secondary suspension
N·s/m
K
pz/Ksz
Stiffness of primary suspension and secondary suspension
N/m
l
t
Wheel base
m
l
c
Length between bogie centers
m
G
Contact constant of wheel and rail
m/N2/3
δZ(t)
Elastic compression deformation
m
Z
w(t)
Wheel displacement
m
Z
r(t)
Rail displacement
m
Z
0(t)
Track irregularity
m
z
r(x, t)/yr(x, t)/ϕr(x, t)
Vertical, lateral and torsional displacements of the rail
m
m
r
Rail mass per unit length
kg/m
ρ
r
Rail density
kg/m3
x
i
/
x
w
j
Locations ofith fastener and jth wheelset
m
N
s/Nw
Numbers of fasteners and wheelsets
—
δ(x)
Dirac delta function
—


### 2.1. 2D Train Subsystem

The train system contains a series of vehicles, which are evenly placed on the track. Each vehicle is modeled as a mass-spring-damper system consisting of a car body, two bogie frames, four wheelsets, and two-stage suspensions, as shown in Figure1. Each vehicle has 10 DOFs, including the vertical and pitch motions of the carbody and two bogie frames, and the vertical motion of each wheelset. The dynamic equations for all the seven bodies are given as follows.Vertical motion of car body:(1)McZ¨c+2CszZ˙c+2KszZc−CszZ˙t1−KszZt1−CszZ˙t2−KszZt2=Mcg.Pitch motion of car body:(2)Jcβ¨c+2Cszlc2β˙c+2Kszlc2βc+CszlcZ˙t1−CszlcZ˙t2+KszlcZt1−KszlcZt2=0.Vertical motion of front frame:(3)MtZ¨t1+2Cpz+CszZ˙t1+2Kpz+KszZt1−CszZ˙c−KszZc−CpzZ˙w1−CpzZ˙w2−KpzZw1−KpzZw2+Cszlcβ˙c+Kszlcβc=Mtg.Pitch motion of front frame:(4)Jtβ¨t1+2Cpzlt2β˙t1+2Kpzlt2βt1+CpzltZ˙w1−CpzltZ˙w2+KpzlcZw1−KpzlcZw2=0.Vertical motion of rear frame:(5)MtZ¨t2+2Cpz+CszZ˙t2+2Kpz+KszZt2−CszZ˙c−KszZc−CpzZ˙w3−CpzZ˙w4−KpzZw3−KpzZw4−Cszlcβ˙c+Kszlcβc=Mtg.Pitch motion of rear frame:(6)Jtβ¨t2+2Cpzlt2β˙t2+2Kpzlt2βt2+CpzltZ˙w3−CpzltZ˙w4+KpzlcZw3−KpzlcZw4=0.Vertical motion of 1st wheelset:(7)MwZ¨w1+CpzZ˙w1+KpzZw1−CpzZ˙t1−KpzZt1+Cpzltβ˙t1+Kpzltβt1+2p1t−Mwg=F01t.Vertical motion of 2nd wheelset:(8)MwZ¨w2+CpzZ˙w2+KpzZw2−CpzZ˙t1−KpzZt1+Cpzltβ˙t1−Kpzltβt1+2p2t−Mwg=F02t.Vertical motion of 3rd wheelset:(9)MwZ¨w3+CpzZ˙w3+KpzZw3−CpzZ˙t2−KpzZt2+Cpzltβ˙t2−Kpzltβt2+2p3t−Mwg=F03t.Vertical motion of 4th wheelset:(10)MwZ¨w4+CpzZ˙w4+KpzZw4−CpzZ˙t2−KpzZt2+Cpzltβ˙t2−Kpzltβt2+2p4t−Mwg=F04t.Furthermore, the dynamic equations of the train system can be given in matrix form as follows:(11)MZ¨+CZ˙+KZ=F.The detailed expressions of the mass, the stiffness, and the damping matrices of the train system can be referred to literature [8]. Further, this equation can be simplified as(12)Z¨+M−1CZ˙+M−1KZ=M−1F.These matrices of the train system cannot be calculated by ANSYS automatically, which should be written into the software manually. The relevant program code is given as follows.! Definition of the mass, the stiffness, and the damping matrices of the train subsystem∗dim,C,array,10,10∗dim,K,array,10,10∗dim,F,array,10∗dim,p,array,4! Damping matrixC(1,1)=2∗Csz/Mc$C(1,3)=-Csz/Mc$C(1,5)=-Csz/McC(2,2)=2∗Csz∗lc∗∗2/Jc$C(2,3)=Csz∗lc/Jc$C(2,5)=-Csz∗lc/JcC(3,1)=-Csz/Mt$C(3,2)=Csz∗lc/Mt$C(3,3)=(Csz+2∗Cpz)/Mt$C(3,7)=-Cpz/Mt$C(3,8)=-Cpz/MtC(4,4)=2∗Cpz∗lt∗∗2/Jt$C(4,7)=Cpz∗lt/Jt$C(4,8)=-Cpz∗lt/JtC(5,1)=-Csz/Mt$C(5,2)=-Csz∗lc/Mt$C(5,5)=(Csz+2∗Cpz)/Mt$C(5,9)=-Cpz/Mt$C(5,10)=-Cpz/MtC(6,6)=2∗Cpz∗lt∗∗2/Jt$C(6,9)=Cpz∗lt/Jt$C(6,10)=-Cpz∗lt/JtC(7,3)=-Cpz/Mw$C(7,4)=Cpz∗lt/Mw$C(7,7)=Cpz/MwC(8,3)=-Cpz/Mw$C(8,4)=-Cpz∗lt/Mw$C(8,8)=Cpz/MwC(9,5)=-Cpz/Mw$C(9,6)=Cpz∗lt/Mw$C(9,9)=Cpz/MwC(10,5)=-Cpz/Mw$C(10,6)=-Cpz∗lt/Mw$C(10,10)=Cpz/Mw! Stiffness matrixK(1,1)=2∗Ksz/Mc$K(1,3)=-Ksz/Mc$K(1,5)=-Ksz/McK(2,2)=2∗Ksz∗lc∗∗2/Jc$K(2,3)=Ksz∗lc/Jc$K(2,5)=-Ksz∗lc/JcK(3,1)=-Ksz/Mt$K(3,2)=Ksz∗lc/Mt$K(3,3)=(Ksz+2∗Kpz)/Mt$K(3,7)=-Kpz/Mt$K(3,8)=-Kpz/MtK(4,4)=2∗Kpz∗lt∗∗2/Jt$K(4,7)=Kpz∗lt/Jt$K(4,8)=-Kpz∗lt/JtK(5,1)=-Ksz/Mt$K(5,2)=-Ksz∗lc/Mt$K(5,5)=(Ksz+2∗Kpz)/Mt$K(5,9)=-Kpz/Mt$K(5,10)=-Kpz/MtK(6,6)=2∗Kpz∗lt∗∗2/Jt$K(6,9)=Kpz∗lt/Jt$K(6,10)=-Kpz∗lt/JtK(7,3)=-Kpz/Mw$K(7,4)=Kpz∗lt/Mw$K(7,7)=Kpz/MwK(8,3)=-Kpz/Mw$K(8,4)=-Kpz∗lt/Mw$K(8,8)=Kpz/MwK(9,5)=-Kpz/Mw$K(9,6)=Kpz∗lt/Mw$K(9,9)=Kpz/MwK(10,5)=-Kpz/Mw$K(10,6)=-Kpz∗lt/Mw$K(10,10)=Kpz/Mw! Loading vectorF(1)=Mc∗g/Mc$F(2)=0$F(3)=Mt∗g/Mt$F(4)=0$F(5)=Mt∗g/Mt$F(6)=0F(7)=(Mw∗g-p(1))/Mw$F(8)=(Mw∗g-p(2))/Mw$F(9)=(Mw∗g-p(3))/Mw$F(10)=(Mw∗g-p(4))/MwIn the above code,C, K, and F denote damping matrix, stiffness matrix, and load vector, respectively; Csz and Cpz are the damping in the secondary suspension and primary suspension, respectively; Ksz and Kpz are the stiffness in the secondary suspension and primary suspension, respectively; lc and lt denote the half of the distance between two bogies on one vehicle and the half of the distance between two wheelsets on one bogie; Mc, Mt, and Mw represent the mass of the carbody, the bogie frame, and the wheelset, respectively.


### 2.2. 2D Track Subsystem

The track system consists of many parts with different mechanical characteristics. ANSYS provides different kinds of elements and materials, which are enough to accurately simulate different components in track structures. According to the modelling principle in research [19], the rails can be modeled with Timoshenko-beam element BEAM188, the fastener systems are usually simulated by spring-damper element COMBIN14, and the slabs are established adopting 3D solid element SOLID185. The code for establishing rail and fasteners is given by! Definition of real constantsr,1,40e6,40e3! for fastenerr,2,0.007745,0.524e-5,0.3217e-4,0.176,0.1115! for rail! Modeling of railk,1,0,0,0k,2,0,0,Lengthk,3,0,1,0l,1,2latt,1,2,1,3,esize, L_raillmesh,all! Modeling of fasteners∗do,i,1,Length/L_space+1k,10+i,0,-0.1,(i-1)∗L_space∗enddotype,2real,1l, kp(0,0,0), kp(0,-0.1,0)lsel,s,loc,z,0lesize,all,1lmesh,alllgen, Length/L_space+1,all,,L_spaceIn the above code,Length, L_rail, and L_space are the total length of the rail, the length of each element on rail, and the fastener space.


### 2.3. 2D Train-Track Interaction

The train system and the track system are coupled by wheel-rail interaction relationship, which can be described by the nonlinear Hertz contact theory [20, 21]. Therefore, the wheel-rail force is determined by(13)pt=1GδZt3/2.Furthermore,δZt is written as(14)δZt=Zwt−Zrt−Z0t.It should be underlined that the wheel-rail contact relationship is programmed into ANSYS rather than employing the default contact in this software.The rail is modeled by the finite element method; thus, all the wheels will not be at the node locations exactly in each integration step. In most cases, the wheels are located between the two adjacent nodes. Therefore, it is necessary to allocate the wheel-rail force to the two adjacent nodes, as shown in Figure2. The wheel-rail force Pw is applied on any point between the ith node and the (i + 1)th node on the beam element.Figure 2
Allocation of wheel-rail force on the beam element of rail.As shown in Figure2, the wheel-rail force applying on the beam element can be allocated to two forces and two bending moments on the adjacent two nodes, as expressed in the following:(15)pwi=pwb2l+2al3,Mwi=pwab2l2,pwi+1=pwa2l+2bl3,Mwi+1=pwa2bl2,where pwi and pwi+1 are the forces allocated on the ith node and the (i + 1)th node, respectively; Mwi and Mwi+1 are the bending moments allocated on the ith node and the (i + 1)th node, respectively; and l is the length of the beam.Also, in the finite element model, the vibration displacements of all the nodes are calculated. In order to obtain the rail displacement at the wheelset location between the adjacent two nodes, the following equation is employed based on the interpolation theory:(16)Zr=S1Zi+S2Ri+S3Zi+1+S4Ri+1,where Zr is the rail displacement at the wheel location; Zi and Zi+1 are the displacements of the ith node and the (i + 1)th node, respectively; Ri and Ri+1 are the rotation angles of the ith node and the (i + 1)th node, respectively; and S1, S2, S3, and S4 are the coefficients, which can be written as Hermite shape functions as follows:(17)S1=1−3a2l2+2a3l3,S2=a−2a2l+a3l2,S3=3a2l2−2a3l3,S4=−a2l+a3l2.As a result, the rail displacement at any point can be obtained combining equations (16) and (17).The allocation of wheel-rail force in APDL language is written as∗do,i,1,4! Determination of rail nodes applied by wheel-rail force∗if,nint(xw(i)/L_rail),GE,xw(i)/L_rail, thenN_front(i)=node_left(nint(xw(i)/L_rail))N_back(i)=node_left(nint(xw(i)/L_rail)-1)L_front(i)=nint(xw(i)/L_rail)∗L_rail-xw(i)L_back(i)=L_rail-L_front(i)∗elseN_front(i)=node_left(nint(xw(i)/L_rail)+1)N_back(i)=node_left(nint(xw(i)/L_rail))L_back(i)=xw(i)-nint(xw(i)/L_rail)∗L_railL_front(i)=L_rail-L_back(i)∗endifallsel,all,all! Allocation of wheel-rail forcep_front(i)=(L_back(i)+3∗L_front(i))∗p(i)∗L_back(i)∗∗2/L_rail∗∗3p_back(i)=(3∗L_back(i)+L_front(i))∗p(i)∗L_front(i)∗∗2/L_rail∗∗3M_front(i)=p(i)∗L_back(i)∗∗2∗L_front(i)/L_rail∗∗2M_back(i)=p(i)∗L_back(i)∗L_front(i)∗∗2/L_rail∗∗2f,N_front(i),fy,-1∗p_front(i)f,N_back(i),fy,-1∗p_back(i)f,N_front(i),mx,-1∗M_front(i)f,N_back(i),mx,M_back(i)∗enddoIn the above program,front and back denote the two adjacent rail nodes, P and M are the allocated force and moment, L is the length between wheelset location and rail node, and N is the adjacent rail node number. Using this above code, the wheel-rail force can be easily allocated to adjacent rail nodes and the dynamic responses of track structures can be determined.Moreover, for other finite element software such as ADINA and ABAQUS, it is rather difficult to consider track random irregularity in investigating TTI, leading to the calculation processes of train system and track system being usually separate. For instance, precalculated wheel-rail forces, which are calculated with other tools, are applied on the track FE model to investigate the track vibrations.As for ANSYS, the APDL provides a convenient program language to implement secondary development for ANSYS. After the displacements of train and track systems are determined, the wheel-rail forces can be further calculated by considering track irregularity and the relevant program is given as follows. In the program,Irr is the formed irregularity table, del_t is the integration step, i and j denote the number of wheelsets on one vehicle and the integration step, del_z is the elastic compression deformation between wheel and rail at the same position, and X and Zr are the displacements of train system and rail, respectively.! Elastic compression deformation between wheel and raildel_z(1)=X(7)-zr(1)-Irr(j∗del_t∗vel+2∗lc+2∗lt)del_z(2)=X(8)-zr(2)-Irr(j∗del_t∗vel+2∗lc)del_z(3)=X(9)-zr(3)-Irr(j∗del_t∗vel+2∗lt)del_z(4)=X(10)-zr(4)-Irr(j∗del_t∗vel)! Calculation of wheel-rail force∗do,i,1,4∗if,del_z(i),GT,0,thenp(i)=(del_z(i)/GG)∗∗1.5∗elsep(i)=0∗endifF(6+i)=(Mw∗g-p(i))/Mw∗enddo


### 2.4. 3D Train-Track Dynamic System

The model and written code for 2D train-track dynamic system are introduced in Section2.3 and 3D train-track dynamic model is briefly described in this section, which is much more complex.The train submodel is built by considering seven rigid bodies of a car body, two bogies, and four wheelsets. Five DOFs are taken into consideration for each rigid body, describing bounce, sway, roll, yaw, and pitch motions. The detailed equations of motion of all the seven parts can be referred to literature [22, 23].The rails are modeled as Euler beams and discretely supported by fasteners which are simulated as linear spring-damping elements. Three DOFs of each rail are taken into account, describing vertical motion, lateral motion, and torsional motion, and the equations of motion of the rail are given as follows:Vertical motion:(18)EIy∂4zrx,t∂x4+mr∂2zrx,t∂t2=−∑i=1NsFrVitδx−xi+∑j=1NwPjtδx−xwj.Lateral motion:(19)EIz∂4yrx,t∂x4+mr∂2yrx,t∂t2=−∑i=1NsFrLitδx−xi+∑j=1NwQjtδx−xwj.Torsional motion:(20)GJt∂2ϕrx,t∂x2+ρrIr0∂2ϕrx,t∂t2=−∑i=1NsMFitδx−xi+∑j=1NwMwjtδx−xwj.In equations (18)–(20), zr(x, t), yr(x, t), and ϕr(x, t) are vertical, lateral, and torsional displacements of the rail, respectively; mr is the rail mass per unit length; ρr is the rail density; EIy and EIz are the rail bending stiffness to Y-axle and Z-axle, respectively; Ir0 is the torsional inertia of the rail; GJt is the rail torsional stiffness; FrVi(t) and FrLi(t) are the vertical and lateral dynamic forces of the ith fastener; Pj(t) and Qj(t) are the jth wheel-rail vertical and lateral forces; MFi(t) and Mwj(t) are the moments applying on the rails due to forces FrVi(t) and FrLi(t) and due to forces Pj(t) and Qj(t), respectively; xi and xwj are the locations of ith fastener and jth wheelset; Ns and Nw are the numbers of fasteners and wheelsets; and δ(x) is the Dirac delta function.In the 3D train-track dynamic model, the nonlinear Hertzian elastic contact theory is used to calculate the wheel-rail normal contact forces:(21)PNt=1GδNt3/2,where G is the wheel-rail contact constant and δN(t) is the wheel-rail normal elastic compression deformation. The left and right wheel-rail normal elastic compression deformation δNL and δNR can be expressed as(22)δNL=δZLcosδL+ϕw,δNR=δZRcosδR−ϕw,where δL and δR are the left and right wheel-rail contact angles, respectively; ϕw denotes the wheelset roll angle; δZL and δZR are the left and right wheel-rail vertical relative displacements, which depend on the motion of wheelset and rail and the wheel-rail contact state.Obviously, whenδZL or δZR is less than zero indicating that the wheel at the left or right side separates from the rail, the left or right side wheel-rail normal contact force is then equal to zero.Based on Kalker’s linear creep theory, the wheel-rail longitudinal creep forceFx, the lateral creep force Fy, and the spin creep torque Mz can be described by(23)Fx=−f11ξx,Fy=−f22ξy−f23ξsp,Mz=f23ξy−f33ξsp,where fij represents the creep coefficients, and the definitions of the longitudinal, lateral, and spin creepages ξx, ξy, and ξsp are as follows:(24)ξx=Vw1−Vr1V,ξy=Vw2−Vr2V,ξsp=Ωw3−Ωr3V,in which Vw1, Vw2, and Ωw3 are the longitudinal, lateral, and spin speeds of the wheel at contact point, respectively; Vr1, Vr2, and Ωr3 are the speeds of the rail at the contact point; and V is the wheelset speed moving on rails.Since Kalker’s linear creep theory is only suited to the cases with small creepages, the nonlinear modification should be made, for example, by the Shen–Hedrick–Elkins model, for the situations of large creepages. Similar to the 2D train-track model, the 3D wheel-rail relationship is also programmed into ANSYS rather than employing the default contact in the software.The wheel-rail interaction forces mainly include wheel-rail normal contact forces derived by nonlinear Hertzian elastic contact theory according to the elastic compression deformation of wheels and rails at contact points in the normal directions and tangential wheel-rail creep forces, which are calculated first by Kalker’s linear creep theory and then modified with the Shen-Hedrick–Elkins nonlinear model [22, 23].Due to the space limitation, the code for the 3D train-track dynamic model is not listed below, which can be referred in the first author’s doctor degree dissertation [24].


### 2.5. Numerical Integration Method

Through the above explanation, the whole dynamic equations of the train-track system are established, which can be solved by the step-by-step integration method. In this present calculation model, the track is modeled by FEM, while the vehicle system and the wheel-rail relationship are programmed into ANSYS by secondary development with APDL language. The Newmark implicit integration method is selected to calculate the FEM, which ensures the stability of the FEM calculation. In order to improve the computational efficiency, the Zhai method [25], an explicit integration method, is adopted to solve the train model, given as(25)Xn+1=Xn+VnΔt+0.5+ψAnΔt2−ψAn−1Δt2,Vn+1=Vn+1+φAnΔt−φAn−1Δt,where X, V, and A are the generalized displacement, velocity, and acceleration of the system, respectively; ∆t is the time step; the subscripts (n − 1), n, and (n + 1) represent the previous two steps (t = (n − 1) ∆t, t = n∆t) and the present step (t = (n + 1) ∆t); and φ and ψ are independent parameters that control the stability of the algorithm.The frame of Zhai method is implemented in ANSYS as∗do,j,2,outt=t+del_ttime,tnsubst,1!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Determination of wheelset locations!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗do,i,1,10X(i)=X(i)+V(i)∗del_t+(0.5+bosai)∗A(i)∗del_t∗∗2-bosai∗AA(i)∗del_t∗∗2V(i)=V(i)+(1+fai)∗A(i)∗del_t-fai∗AA(i)∗del_t∗enddo∗do,i,1,10AA(i)=A(i)∗enddo!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Determination of rail nodes!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Allocation of wheel-rail forces!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗solve!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Acquirement of rail displacements!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Calculation of wheel-rail forces!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗MOPER,mid1,K,MULT,X∗MOPER,mid2,C,MULT,V∗do,i,1,10A(i)=F(i)-mid1(i)-mid2(i)∗enddo∗enddowhere X, V, and A represent the displacement, velocity, and acceleration vectors of the train system in the current integration step and AA is the acceleration vectors of the train in the previous step. This code only describes the frame of the Zhai method, and the program can be implemented after the annotated parts are written with relevant codes.As seen from the code, in the beginning of each integration step, the displacement and velocity of the train system are calculated, while the acceleration of the train is determined in the end of each integration step.


### 2.6. Procedure of the Proposed Method

The methodology of investigating TTI with ANSYS is shown in Figure3. In this procedure, t, ∆t, and T denote the integration time, integration step, and total calculation time, respectively. According to this procedure, TTI can be investigated in the platform of ANSYS.Figure 3
Procedure of investigating TTI with ANSYS adopting APDL language.


## 2.1. 2D Train Subsystem

The train system contains a series of vehicles, which are evenly placed on the track. Each vehicle is modeled as a mass-spring-damper system consisting of a car body, two bogie frames, four wheelsets, and two-stage suspensions, as shown in Figure1. Each vehicle has 10 DOFs, including the vertical and pitch motions of the carbody and two bogie frames, and the vertical motion of each wheelset. The dynamic equations for all the seven bodies are given as follows.Vertical motion of car body:(1)McZ¨c+2CszZ˙c+2KszZc−CszZ˙t1−KszZt1−CszZ˙t2−KszZt2=Mcg.Pitch motion of car body:(2)Jcβ¨c+2Cszlc2β˙c+2Kszlc2βc+CszlcZ˙t1−CszlcZ˙t2+KszlcZt1−KszlcZt2=0.Vertical motion of front frame:(3)MtZ¨t1+2Cpz+CszZ˙t1+2Kpz+KszZt1−CszZ˙c−KszZc−CpzZ˙w1−CpzZ˙w2−KpzZw1−KpzZw2+Cszlcβ˙c+Kszlcβc=Mtg.Pitch motion of front frame:(4)Jtβ¨t1+2Cpzlt2β˙t1+2Kpzlt2βt1+CpzltZ˙w1−CpzltZ˙w2+KpzlcZw1−KpzlcZw2=0.Vertical motion of rear frame:(5)MtZ¨t2+2Cpz+CszZ˙t2+2Kpz+KszZt2−CszZ˙c−KszZc−CpzZ˙w3−CpzZ˙w4−KpzZw3−KpzZw4−Cszlcβ˙c+Kszlcβc=Mtg.Pitch motion of rear frame:(6)Jtβ¨t2+2Cpzlt2β˙t2+2Kpzlt2βt2+CpzltZ˙w3−CpzltZ˙w4+KpzlcZw3−KpzlcZw4=0.Vertical motion of 1st wheelset:(7)MwZ¨w1+CpzZ˙w1+KpzZw1−CpzZ˙t1−KpzZt1+Cpzltβ˙t1+Kpzltβt1+2p1t−Mwg=F01t.Vertical motion of 2nd wheelset:(8)MwZ¨w2+CpzZ˙w2+KpzZw2−CpzZ˙t1−KpzZt1+Cpzltβ˙t1−Kpzltβt1+2p2t−Mwg=F02t.Vertical motion of 3rd wheelset:(9)MwZ¨w3+CpzZ˙w3+KpzZw3−CpzZ˙t2−KpzZt2+Cpzltβ˙t2−Kpzltβt2+2p3t−Mwg=F03t.Vertical motion of 4th wheelset:(10)MwZ¨w4+CpzZ˙w4+KpzZw4−CpzZ˙t2−KpzZt2+Cpzltβ˙t2−Kpzltβt2+2p4t−Mwg=F04t.Furthermore, the dynamic equations of the train system can be given in matrix form as follows:(11)MZ¨+CZ˙+KZ=F.The detailed expressions of the mass, the stiffness, and the damping matrices of the train system can be referred to literature [8]. Further, this equation can be simplified as(12)Z¨+M−1CZ˙+M−1KZ=M−1F.These matrices of the train system cannot be calculated by ANSYS automatically, which should be written into the software manually. The relevant program code is given as follows.! Definition of the mass, the stiffness, and the damping matrices of the train subsystem∗dim,C,array,10,10∗dim,K,array,10,10∗dim,F,array,10∗dim,p,array,4! Damping matrixC(1,1)=2∗Csz/Mc$C(1,3)=-Csz/Mc$C(1,5)=-Csz/McC(2,2)=2∗Csz∗lc∗∗2/Jc$C(2,3)=Csz∗lc/Jc$C(2,5)=-Csz∗lc/JcC(3,1)=-Csz/Mt$C(3,2)=Csz∗lc/Mt$C(3,3)=(Csz+2∗Cpz)/Mt$C(3,7)=-Cpz/Mt$C(3,8)=-Cpz/MtC(4,4)=2∗Cpz∗lt∗∗2/Jt$C(4,7)=Cpz∗lt/Jt$C(4,8)=-Cpz∗lt/JtC(5,1)=-Csz/Mt$C(5,2)=-Csz∗lc/Mt$C(5,5)=(Csz+2∗Cpz)/Mt$C(5,9)=-Cpz/Mt$C(5,10)=-Cpz/MtC(6,6)=2∗Cpz∗lt∗∗2/Jt$C(6,9)=Cpz∗lt/Jt$C(6,10)=-Cpz∗lt/JtC(7,3)=-Cpz/Mw$C(7,4)=Cpz∗lt/Mw$C(7,7)=Cpz/MwC(8,3)=-Cpz/Mw$C(8,4)=-Cpz∗lt/Mw$C(8,8)=Cpz/MwC(9,5)=-Cpz/Mw$C(9,6)=Cpz∗lt/Mw$C(9,9)=Cpz/MwC(10,5)=-Cpz/Mw$C(10,6)=-Cpz∗lt/Mw$C(10,10)=Cpz/Mw! Stiffness matrixK(1,1)=2∗Ksz/Mc$K(1,3)=-Ksz/Mc$K(1,5)=-Ksz/McK(2,2)=2∗Ksz∗lc∗∗2/Jc$K(2,3)=Ksz∗lc/Jc$K(2,5)=-Ksz∗lc/JcK(3,1)=-Ksz/Mt$K(3,2)=Ksz∗lc/Mt$K(3,3)=(Ksz+2∗Kpz)/Mt$K(3,7)=-Kpz/Mt$K(3,8)=-Kpz/MtK(4,4)=2∗Kpz∗lt∗∗2/Jt$K(4,7)=Kpz∗lt/Jt$K(4,8)=-Kpz∗lt/JtK(5,1)=-Ksz/Mt$K(5,2)=-Ksz∗lc/Mt$K(5,5)=(Ksz+2∗Kpz)/Mt$K(5,9)=-Kpz/Mt$K(5,10)=-Kpz/MtK(6,6)=2∗Kpz∗lt∗∗2/Jt$K(6,9)=Kpz∗lt/Jt$K(6,10)=-Kpz∗lt/JtK(7,3)=-Kpz/Mw$K(7,4)=Kpz∗lt/Mw$K(7,7)=Kpz/MwK(8,3)=-Kpz/Mw$K(8,4)=-Kpz∗lt/Mw$K(8,8)=Kpz/MwK(9,5)=-Kpz/Mw$K(9,6)=Kpz∗lt/Mw$K(9,9)=Kpz/MwK(10,5)=-Kpz/Mw$K(10,6)=-Kpz∗lt/Mw$K(10,10)=Kpz/Mw! Loading vectorF(1)=Mc∗g/Mc$F(2)=0$F(3)=Mt∗g/Mt$F(4)=0$F(5)=Mt∗g/Mt$F(6)=0F(7)=(Mw∗g-p(1))/Mw$F(8)=(Mw∗g-p(2))/Mw$F(9)=(Mw∗g-p(3))/Mw$F(10)=(Mw∗g-p(4))/MwIn the above code,C, K, and F denote damping matrix, stiffness matrix, and load vector, respectively; Csz and Cpz are the damping in the secondary suspension and primary suspension, respectively; Ksz and Kpz are the stiffness in the secondary suspension and primary suspension, respectively; lc and lt denote the half of the distance between two bogies on one vehicle and the half of the distance between two wheelsets on one bogie; Mc, Mt, and Mw represent the mass of the carbody, the bogie frame, and the wheelset, respectively.


## 2.2. 2D Track Subsystem

The track system consists of many parts with different mechanical characteristics. ANSYS provides different kinds of elements and materials, which are enough to accurately simulate different components in track structures. According to the modelling principle in research [19], the rails can be modeled with Timoshenko-beam element BEAM188, the fastener systems are usually simulated by spring-damper element COMBIN14, and the slabs are established adopting 3D solid element SOLID185. The code for establishing rail and fasteners is given by! Definition of real constantsr,1,40e6,40e3! for fastenerr,2,0.007745,0.524e-5,0.3217e-4,0.176,0.1115! for rail! Modeling of railk,1,0,0,0k,2,0,0,Lengthk,3,0,1,0l,1,2latt,1,2,1,3,esize, L_raillmesh,all! Modeling of fasteners∗do,i,1,Length/L_space+1k,10+i,0,-0.1,(i-1)∗L_space∗enddotype,2real,1l, kp(0,0,0), kp(0,-0.1,0)lsel,s,loc,z,0lesize,all,1lmesh,alllgen, Length/L_space+1,all,,L_spaceIn the above code,Length, L_rail, and L_space are the total length of the rail, the length of each element on rail, and the fastener space.


## 2.3. 2D Train-Track Interaction

The train system and the track system are coupled by wheel-rail interaction relationship, which can be described by the nonlinear Hertz contact theory [20, 21]. Therefore, the wheel-rail force is determined by(13)pt=1GδZt3/2.Furthermore,δZt is written as(14)δZt=Zwt−Zrt−Z0t.It should be underlined that the wheel-rail contact relationship is programmed into ANSYS rather than employing the default contact in this software.The rail is modeled by the finite element method; thus, all the wheels will not be at the node locations exactly in each integration step. In most cases, the wheels are located between the two adjacent nodes. Therefore, it is necessary to allocate the wheel-rail force to the two adjacent nodes, as shown in Figure2. The wheel-rail force Pw is applied on any point between the ith node and the (i + 1)th node on the beam element.Figure 2
Allocation of wheel-rail force on the beam element of rail.As shown in Figure2, the wheel-rail force applying on the beam element can be allocated to two forces and two bending moments on the adjacent two nodes, as expressed in the following:(15)pwi=pwb2l+2al3,Mwi=pwab2l2,pwi+1=pwa2l+2bl3,Mwi+1=pwa2bl2,where pwi and pwi+1 are the forces allocated on the ith node and the (i + 1)th node, respectively; Mwi and Mwi+1 are the bending moments allocated on the ith node and the (i + 1)th node, respectively; and l is the length of the beam.Also, in the finite element model, the vibration displacements of all the nodes are calculated. In order to obtain the rail displacement at the wheelset location between the adjacent two nodes, the following equation is employed based on the interpolation theory:(16)Zr=S1Zi+S2Ri+S3Zi+1+S4Ri+1,where Zr is the rail displacement at the wheel location; Zi and Zi+1 are the displacements of the ith node and the (i + 1)th node, respectively; Ri and Ri+1 are the rotation angles of the ith node and the (i + 1)th node, respectively; and S1, S2, S3, and S4 are the coefficients, which can be written as Hermite shape functions as follows:(17)S1=1−3a2l2+2a3l3,S2=a−2a2l+a3l2,S3=3a2l2−2a3l3,S4=−a2l+a3l2.As a result, the rail displacement at any point can be obtained combining equations (16) and (17).The allocation of wheel-rail force in APDL language is written as∗do,i,1,4! Determination of rail nodes applied by wheel-rail force∗if,nint(xw(i)/L_rail),GE,xw(i)/L_rail, thenN_front(i)=node_left(nint(xw(i)/L_rail))N_back(i)=node_left(nint(xw(i)/L_rail)-1)L_front(i)=nint(xw(i)/L_rail)∗L_rail-xw(i)L_back(i)=L_rail-L_front(i)∗elseN_front(i)=node_left(nint(xw(i)/L_rail)+1)N_back(i)=node_left(nint(xw(i)/L_rail))L_back(i)=xw(i)-nint(xw(i)/L_rail)∗L_railL_front(i)=L_rail-L_back(i)∗endifallsel,all,all! Allocation of wheel-rail forcep_front(i)=(L_back(i)+3∗L_front(i))∗p(i)∗L_back(i)∗∗2/L_rail∗∗3p_back(i)=(3∗L_back(i)+L_front(i))∗p(i)∗L_front(i)∗∗2/L_rail∗∗3M_front(i)=p(i)∗L_back(i)∗∗2∗L_front(i)/L_rail∗∗2M_back(i)=p(i)∗L_back(i)∗L_front(i)∗∗2/L_rail∗∗2f,N_front(i),fy,-1∗p_front(i)f,N_back(i),fy,-1∗p_back(i)f,N_front(i),mx,-1∗M_front(i)f,N_back(i),mx,M_back(i)∗enddoIn the above program,front and back denote the two adjacent rail nodes, P and M are the allocated force and moment, L is the length between wheelset location and rail node, and N is the adjacent rail node number. Using this above code, the wheel-rail force can be easily allocated to adjacent rail nodes and the dynamic responses of track structures can be determined.Moreover, for other finite element software such as ADINA and ABAQUS, it is rather difficult to consider track random irregularity in investigating TTI, leading to the calculation processes of train system and track system being usually separate. For instance, precalculated wheel-rail forces, which are calculated with other tools, are applied on the track FE model to investigate the track vibrations.As for ANSYS, the APDL provides a convenient program language to implement secondary development for ANSYS. After the displacements of train and track systems are determined, the wheel-rail forces can be further calculated by considering track irregularity and the relevant program is given as follows. In the program,Irr is the formed irregularity table, del_t is the integration step, i and j denote the number of wheelsets on one vehicle and the integration step, del_z is the elastic compression deformation between wheel and rail at the same position, and X and Zr are the displacements of train system and rail, respectively.! Elastic compression deformation between wheel and raildel_z(1)=X(7)-zr(1)-Irr(j∗del_t∗vel+2∗lc+2∗lt)del_z(2)=X(8)-zr(2)-Irr(j∗del_t∗vel+2∗lc)del_z(3)=X(9)-zr(3)-Irr(j∗del_t∗vel+2∗lt)del_z(4)=X(10)-zr(4)-Irr(j∗del_t∗vel)! Calculation of wheel-rail force∗do,i,1,4∗if,del_z(i),GT,0,thenp(i)=(del_z(i)/GG)∗∗1.5∗elsep(i)=0∗endifF(6+i)=(Mw∗g-p(i))/Mw∗enddo


## 2.4. 3D Train-Track Dynamic System

The model and written code for 2D train-track dynamic system are introduced in Section2.3 and 3D train-track dynamic model is briefly described in this section, which is much more complex.The train submodel is built by considering seven rigid bodies of a car body, two bogies, and four wheelsets. Five DOFs are taken into consideration for each rigid body, describing bounce, sway, roll, yaw, and pitch motions. The detailed equations of motion of all the seven parts can be referred to literature [22, 23].The rails are modeled as Euler beams and discretely supported by fasteners which are simulated as linear spring-damping elements. Three DOFs of each rail are taken into account, describing vertical motion, lateral motion, and torsional motion, and the equations of motion of the rail are given as follows:Vertical motion:(18)EIy∂4zrx,t∂x4+mr∂2zrx,t∂t2=−∑i=1NsFrVitδx−xi+∑j=1NwPjtδx−xwj.Lateral motion:(19)EIz∂4yrx,t∂x4+mr∂2yrx,t∂t2=−∑i=1NsFrLitδx−xi+∑j=1NwQjtδx−xwj.Torsional motion:(20)GJt∂2ϕrx,t∂x2+ρrIr0∂2ϕrx,t∂t2=−∑i=1NsMFitδx−xi+∑j=1NwMwjtδx−xwj.In equations (18)–(20), zr(x, t), yr(x, t), and ϕr(x, t) are vertical, lateral, and torsional displacements of the rail, respectively; mr is the rail mass per unit length; ρr is the rail density; EIy and EIz are the rail bending stiffness to Y-axle and Z-axle, respectively; Ir0 is the torsional inertia of the rail; GJt is the rail torsional stiffness; FrVi(t) and FrLi(t) are the vertical and lateral dynamic forces of the ith fastener; Pj(t) and Qj(t) are the jth wheel-rail vertical and lateral forces; MFi(t) and Mwj(t) are the moments applying on the rails due to forces FrVi(t) and FrLi(t) and due to forces Pj(t) and Qj(t), respectively; xi and xwj are the locations of ith fastener and jth wheelset; Ns and Nw are the numbers of fasteners and wheelsets; and δ(x) is the Dirac delta function.In the 3D train-track dynamic model, the nonlinear Hertzian elastic contact theory is used to calculate the wheel-rail normal contact forces:(21)PNt=1GδNt3/2,where G is the wheel-rail contact constant and δN(t) is the wheel-rail normal elastic compression deformation. The left and right wheel-rail normal elastic compression deformation δNL and δNR can be expressed as(22)δNL=δZLcosδL+ϕw,δNR=δZRcosδR−ϕw,where δL and δR are the left and right wheel-rail contact angles, respectively; ϕw denotes the wheelset roll angle; δZL and δZR are the left and right wheel-rail vertical relative displacements, which depend on the motion of wheelset and rail and the wheel-rail contact state.Obviously, whenδZL or δZR is less than zero indicating that the wheel at the left or right side separates from the rail, the left or right side wheel-rail normal contact force is then equal to zero.Based on Kalker’s linear creep theory, the wheel-rail longitudinal creep forceFx, the lateral creep force Fy, and the spin creep torque Mz can be described by(23)Fx=−f11ξx,Fy=−f22ξy−f23ξsp,Mz=f23ξy−f33ξsp,where fij represents the creep coefficients, and the definitions of the longitudinal, lateral, and spin creepages ξx, ξy, and ξsp are as follows:(24)ξx=Vw1−Vr1V,ξy=Vw2−Vr2V,ξsp=Ωw3−Ωr3V,in which Vw1, Vw2, and Ωw3 are the longitudinal, lateral, and spin speeds of the wheel at contact point, respectively; Vr1, Vr2, and Ωr3 are the speeds of the rail at the contact point; and V is the wheelset speed moving on rails.Since Kalker’s linear creep theory is only suited to the cases with small creepages, the nonlinear modification should be made, for example, by the Shen–Hedrick–Elkins model, for the situations of large creepages. Similar to the 2D train-track model, the 3D wheel-rail relationship is also programmed into ANSYS rather than employing the default contact in the software.The wheel-rail interaction forces mainly include wheel-rail normal contact forces derived by nonlinear Hertzian elastic contact theory according to the elastic compression deformation of wheels and rails at contact points in the normal directions and tangential wheel-rail creep forces, which are calculated first by Kalker’s linear creep theory and then modified with the Shen-Hedrick–Elkins nonlinear model [22, 23].Due to the space limitation, the code for the 3D train-track dynamic model is not listed below, which can be referred in the first author’s doctor degree dissertation [24].


## 2.5. Numerical Integration Method

Through the above explanation, the whole dynamic equations of the train-track system are established, which can be solved by the step-by-step integration method. In this present calculation model, the track is modeled by FEM, while the vehicle system and the wheel-rail relationship are programmed into ANSYS by secondary development with APDL language. The Newmark implicit integration method is selected to calculate the FEM, which ensures the stability of the FEM calculation. In order to improve the computational efficiency, the Zhai method [25], an explicit integration method, is adopted to solve the train model, given as(25)Xn+1=Xn+VnΔt+0.5+ψAnΔt2−ψAn−1Δt2,Vn+1=Vn+1+φAnΔt−φAn−1Δt,where X, V, and A are the generalized displacement, velocity, and acceleration of the system, respectively; ∆t is the time step; the subscripts (n − 1), n, and (n + 1) represent the previous two steps (t = (n − 1) ∆t, t = n∆t) and the present step (t = (n + 1) ∆t); and φ and ψ are independent parameters that control the stability of the algorithm.The frame of Zhai method is implemented in ANSYS as∗do,j,2,outt=t+del_ttime,tnsubst,1!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Determination of wheelset locations!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗do,i,1,10X(i)=X(i)+V(i)∗del_t+(0.5+bosai)∗A(i)∗del_t∗∗2-bosai∗AA(i)∗del_t∗∗2V(i)=V(i)+(1+fai)∗A(i)∗del_t-fai∗AA(i)∗del_t∗enddo∗do,i,1,10AA(i)=A(i)∗enddo!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Determination of rail nodes!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Allocation of wheel-rail forces!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗solve!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Acquirement of rail displacements!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗! Calculation of wheel-rail forces!∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗MOPER,mid1,K,MULT,X∗MOPER,mid2,C,MULT,V∗do,i,1,10A(i)=F(i)-mid1(i)-mid2(i)∗enddo∗enddowhere X, V, and A represent the displacement, velocity, and acceleration vectors of the train system in the current integration step and AA is the acceleration vectors of the train in the previous step. This code only describes the frame of the Zhai method, and the program can be implemented after the annotated parts are written with relevant codes.As seen from the code, in the beginning of each integration step, the displacement and velocity of the train system are calculated, while the acceleration of the train is determined in the end of each integration step.


## 2.6. Procedure of the Proposed Method

The methodology of investigating TTI with ANSYS is shown in Figure3. In this procedure, t, ∆t, and T denote the integration time, integration step, and total calculation time, respectively. According to this procedure, TTI can be investigated in the platform of ANSYS.Figure 3
Procedure of investigating TTI with ANSYS adopting APDL language.


## 3. Validation of the Proposed Calculation Method

The calculation method is introduced in detail, and on this basis, two kinds of validations are conducted to verify its validity.


### 3.1. Validation with Program Written by FORTRAN

Based on vehicle-track dynamics [8], the authors also programmed a procedure of TTI in FORTRAN language on the computing platform of INTER PARALLEL STUDIO to validate the proposed calculation method. The parameters of the vehicle are given in Table 2, and the parameters of the track can be referred in [19].Table 2
Parameters of the adopted vehicle.
Item
Vehicle
Unit
Distance between bogie centers
17.375
m
Bogie wheelbase
2.5
m
Wheel rolling circle diameter
0.92
m
Carbody mass
38.884
t
Bogie frame mass
2953
kg
Wheelset mass
1517
kg
Inertia moment of carbody
1905.3
t·m2
Inertia moment of bogie frame
1293
kg·m2
Inertia moment of wheelset
118
kg·m2
Primary suspension stiffness
0.8865
MN/m
Secondary suspension stiffness
0.195
MN/mSubject to impact loads, such as wheel flat and welding joint, the dynamic responses of the train-track system are illustrated in Figure4.Figure 4
Dynamic responses of the train-track system subject to impact loads: (a) wheel-rail force and (b) rail displacement at the impact location.
(a)
(b)As seen from this result, the wheel-rail forces and rail displacements determined through the two tools are almost the same, indicating the proposed calculation method is effective and able to investigate TTI. It should be noted that the periodic fluctuation of the wheel-rail force from ANSYS is caused by the length of each rail element.


### 3.2. Validation with Commercial Software UM

UM (Universal Mechanism) is developed by the Laboratory of Computational Mechanics in Bryansk State Technical University, which is designed to automate the analysis of mechanical objects representing as a multibody system. By far, the software has been widely adopted in the field of railway engineering worldwide.Excited by lateral and vertical harmonic irregularities, the dynamic responses of the train, whose dynamic parameters can be seen in Table1, are, respectively, calculated by UM and the method proposed in this work, which are listed in Table 3. As seen from the table, the calculated results by UM and the proposed method show good agreement with each other, indicating the proposed calculation method in this work is effective to solve the dynamic problems in the train-track system.Table 3
Calculated responses of CRH3 by UM and the proposed method.
Indicator
Calculated result
By UM
By the proposed method
Vertical acceleration of carbody (m/s2)
±0.133
±0.134
Vertical wheel-rail force (kN)
60.89∼64.11
60.76∼63.97
Lateral acceleration of carbody (m/s2)
±0.142
±0.161
Lateral wheel-rail force (kN)
−1.35∼3.83
−1.39∼4.43


### 3.3. Validation with the Field Test

A 40-day field experimental test of the train-track interaction dynamic was carried out in a certain high-speed railway in 2013 (shown in Figure5), in which the authors were in charge of the experiment on the dynamic performance of the slab ballastless track. In this work, the field test data are adopted to validate this proposed calculation method. In this field experiment, the measured track dynamic indicators consist of the displacement and acceleration of the rail, accelerations of the slab, and base.Figure 5
Field test on a high-speed railway with the slab ballastless track: (a) measurement of the displacement and acceleration of the rail, (b) measurement of the acceleration of the slab, (c) measurement of the acceleration of the base, and (d) the data acquisition system.
(a)
(b)
(c)
(d)The details of the field test, e.g. the track irregularity, the parameters of the train-track system, and the operation speed, can be referred to the authors’ published paper [19].The validation is illustrated in Figure6, from which it can be seen the results obtained from the field test and simulation show good agreement with each other, verifying the validity of the proposed calculation method.Figure 6
Comparison between test results and simulated results: (a) calculated rail acceleration, (b) measured rail acceleration, (c) comparison of rail displacement, and (d) comparison of slab acceleration.
(a)
(b)
(c)
(d)


## 3.1. Validation with Program Written by FORTRAN

Based on vehicle-track dynamics [8], the authors also programmed a procedure of TTI in FORTRAN language on the computing platform of INTER PARALLEL STUDIO to validate the proposed calculation method. The parameters of the vehicle are given in Table 2, and the parameters of the track can be referred in [19].Table 2
Parameters of the adopted vehicle.
Item
Vehicle
Unit
Distance between bogie centers
17.375
m
Bogie wheelbase
2.5
m
Wheel rolling circle diameter
0.92
m
Carbody mass
38.884
t
Bogie frame mass
2953
kg
Wheelset mass
1517
kg
Inertia moment of carbody
1905.3
t·m2
Inertia moment of bogie frame
1293
kg·m2
Inertia moment of wheelset
118
kg·m2
Primary suspension stiffness
0.8865
MN/m
Secondary suspension stiffness
0.195
MN/mSubject to impact loads, such as wheel flat and welding joint, the dynamic responses of the train-track system are illustrated in Figure4.Figure 4
Dynamic responses of the train-track system subject to impact loads: (a) wheel-rail force and (b) rail displacement at the impact location.
(a)
(b)As seen from this result, the wheel-rail forces and rail displacements determined through the two tools are almost the same, indicating the proposed calculation method is effective and able to investigate TTI. It should be noted that the periodic fluctuation of the wheel-rail force from ANSYS is caused by the length of each rail element.


## 3.2. Validation with Commercial Software UM

UM (Universal Mechanism) is developed by the Laboratory of Computational Mechanics in Bryansk State Technical University, which is designed to automate the analysis of mechanical objects representing as a multibody system. By far, the software has been widely adopted in the field of railway engineering worldwide.Excited by lateral and vertical harmonic irregularities, the dynamic responses of the train, whose dynamic parameters can be seen in Table1, are, respectively, calculated by UM and the method proposed in this work, which are listed in Table 3. As seen from the table, the calculated results by UM and the proposed method show good agreement with each other, indicating the proposed calculation method in this work is effective to solve the dynamic problems in the train-track system.Table 3
Calculated responses of CRH3 by UM and the proposed method.
Indicator
Calculated result
By UM
By the proposed method
Vertical acceleration of carbody (m/s2)
±0.133
±0.134
Vertical wheel-rail force (kN)
60.89∼64.11
60.76∼63.97
Lateral acceleration of carbody (m/s2)
±0.142
±0.161
Lateral wheel-rail force (kN)
−1.35∼3.83
−1.39∼4.43


## 3.3. Validation with the Field Test

A 40-day field experimental test of the train-track interaction dynamic was carried out in a certain high-speed railway in 2013 (shown in Figure5), in which the authors were in charge of the experiment on the dynamic performance of the slab ballastless track. In this work, the field test data are adopted to validate this proposed calculation method. In this field experiment, the measured track dynamic indicators consist of the displacement and acceleration of the rail, accelerations of the slab, and base.Figure 5
Field test on a high-speed railway with the slab ballastless track: (a) measurement of the displacement and acceleration of the rail, (b) measurement of the acceleration of the slab, (c) measurement of the acceleration of the base, and (d) the data acquisition system.
(a)
(b)
(c)
(d)The details of the field test, e.g. the track irregularity, the parameters of the train-track system, and the operation speed, can be referred to the authors’ published paper [19].The validation is illustrated in Figure6, from which it can be seen the results obtained from the field test and simulation show good agreement with each other, verifying the validity of the proposed calculation method.Figure 6
Comparison between test results and simulated results: (a) calculated rail acceleration, (b) measured rail acceleration, (c) comparison of rail displacement, and (d) comparison of slab acceleration.
(a)
(b)
(c)
(d)


## 4. Advantages and Disadvantages of the Proposed Calculation Method

Every calculation method certainly has advantages and disadvantages, which will be discussed in this section.The advantages are given as follows:(i)
The track structures can be easily established, even if the structure is abnormal and peculiar(ii)
With the help of finite element software, the nonlinear materials and mechanical behaviors can be easily simulated(iii)
The train and track systems are coupled and calculated simultaneously(iv)
Track irregularity is easily considered(v)
Wind loads and earthquake loads are convenient to be considered in investigating TTI with ANSYS(vi)
The calculation method can also be used in investigating train-bridge interaction, train-track-bridge interaction, train-tunnel interaction, and so onOn the contrary, the disadvantages of the method are listed below:(i)
The computational efficiency is lower than that programmed with FORTRAN, C and so on(ii)
The train model is simpler than the established models in SIMPACK and UM(iii)
Boundary condition is very important in solving the FE model, which should be seriously and accurately considered and appliedCompared to the proposed calculation method with the classified 4 methods (methods A-D) in Section1, the computational accuracy and calculation time of different methods are concluded in Table 4.Table 4
Evaluations of computational accuracy and calculation time of different methods.
Method
Computational accuracy
Calculation time
Description
Level
Description
Level
A
To write vibration equations, complicated structures are simplified in large extent; moreover, inapparent mistakes in the written codes are difficult to find out
Middle
Compiled programs are highly efficient in solving TTI
Fast
B
Rail random irregularity is hard to be accurately considered
Low
Computational efficiency is low due to the implicit integration algorithm in FE software
Slow
C
Nonlinear characteristics of tracks and foundations cannot be considered
Middle
Adopting the algorithm of multibody dynamics and mode superposition method or modal synthesis method, computational efficiency can be ensured
Middle
D
Due to the changes of mass stiffness matrices, nonlinear properties of substructures cannot be considered
Middle
Compiled programs are highly efficient in solving TTI
Fast
Proposed method
Disadvantages of the above methods are settled
High
Hybrid explicit-implicit integration method is adopted to efficiently solve TTI
MiddleAs seen from the above table, the computational accuracy of the proposed method can be improved due to (a) nonlinear characteristics of tracks and foundations can be accurately considered, (b) complicate substructures can be accurately modeled, and (c) track irregularity is easily and accurately considered.


## 5. Application Case

To demonstrate the wide application range of the proposed calculation method, the dynamic responses of a high-speed train running through a bridge-subgrade transition section are investigated using this present method. In previous studies, to explore such a complex problem, the most common method is to simplify the transition section as an additional track irregularity [26], which is different from the actual situation, and the vibrations and dynamic stresses of the track structures cannot be calculated synchronously.A widely used bridge-subgrade transition section with a double-block ballastless track (Figure7) is adopted in this numerical case, whose parameters are listed in Table 5. The dynamic parameters of the train can be found in literature [19]. In the following simulation, the running speed is set to 350 km/h, the length of the transition section is 25 m, and the settlement difference between bridge section and subgrade section is 80 mm.Figure 7
Bridge-subgrade transition section with the double-block ballastless track.Table 5
Key parameters of the bridge-subgrade transition section.
Item
Value
Unit
Stiffness of fastener
30
MN/m
Damping of fastener
3 × 104
N·s/m
Elasticity modulus of slab
3.4 × 104
MPa
Elasticity modulus of support base
0.5 × 104
MPa
Elasticity modulus of subgrade top layer
180
MPa
Elasticity modulus of subgrade bottom layer
110
MPa
Elasticity modulus of embankment
30
MPaAdopting the above parameters, a 2D finite element model of this transition section is displayed in Figure8. In the model, the rail and bridge are modeled by BEAM 188, and the track and infrastructure are simulated with PLANE 42. On this basis, the TTI is investigated according to the calculation method in Section 2, whose simulated results are illustrated in Figures 9–11.Figure 8
FE model of bridge-subgrade transition section in ANSYS.Figure 9
Deformation of rail caused by settlement.Figure 10
Structural dynamic stress of the infrastructure at a certain time.Figure 11
Dynamic responses of the train-transition section system: (a) wheel-rail force and (b) dynamic stress in subgrade bottom layer.
(a)
(b)The rail deformation caused by settlement is displayed in Figure9. As can be clearly seen, the rail moves upward nearby the interface of subgrade and bridge. The deformation curve is not a straight line but an irregular curve, indicating that simplifying this rail deformation as straight lines or sine curve in some published papers [26] is not accurate enough. When investigating TTI in bridge-subgrade transition section, the train system and track model should be coupled and simulated simultaneously.As seen from Figure10, the maximum stress of the infrastructure appears in the transition section under settlement situation, indicating that the transition section is more sensitive to settlement than other sections and this section needs more maintenance. Moreover, the deformations of track and infrastructure due to settlement can also be seen from this figure.Figure11(a) shows the wheel-rail force of the 3rd wheelset when a high-speed vehicle runs through the bridge-subgrade transition section. As seen from the result, when the 1st wheelset runs into the transition section (zone I), the whole vehicle starts to lean backward, resulting in the increase of the wheel-rail force of the 3rd wheelset. Then, after the center of mass of the carbody moves through the location on the deformed rail with the maximum slope (zone II), the wheel-rail force of the 3rd wheelset begins to decrease. Furthermore, when the 3rd wheelset runs into the transition section (zone III), the wheel-rail force increases again. And finally, in zone IV, the force reduces to the normal value after the whole vehicle runs through the transition section. Additionally, Figure 11(b) displays the dynamic stress in the subgrade bottom layer, in which the nonzero initial stress is caused by the settlement. When the vehicle is passing by this area, its influence on the subgrade can be clearly obtained employing this present calculation method.


## 6. Conclusion

Based on secondary development of ANSYS adopting APDL, this paper has proposed an alternative calculation method for investigating train-track interaction (TTI). In this method, the train system and wheel-rail relationship have been programmed into ANSYS based on multibody dynamics, while the track system is established in ANSYS with the help of different kinds of elements. An explicit-implicit hybrid integration method has been employed to solve the large complex dynamic system. Finally, the proposed calculation method has been validated. From this present research, the following conclusions can be reached:(1)
The proposed calculation method is effective to investigate TTI. The whole train-track system is coupled and calculated synchronously.(2)
This paper also provides a guidance to calculate a coupled mechanical system with different integration methods in ANSYS.(3)
This calculation method has a wide application range, which can also be used to study train-track-bridge interaction, train-tunnel interaction, and so on.(4)
The computational efficiency of this present method is lower than that of programming with FORTRAN, C, and so on, indicating that new calculation methods to speed the calculation in ANSYS, such as parallel computing, should be further explored.



---

*Source: 1859261-2019-09-11.xml*