#### 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.

#### 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 Section 2, 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, Figure 1 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.

##### 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 Figure 1. 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:

Pitch motion of car body:

Vertical motion of front frame:

Pitch motion of front frame:

Vertical motion of rear frame:

Pitch motion of rear frame:

Vertical motion of 1^{st} wheelset:

Vertical motion of 2^{nd} wheelset:

Vertical motion of 3^{rd} wheelset:

Vertical motion of 4^{th} wheelset:

Furthermore, the dynamic equations of the train system can be given in matrix form as follows:

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

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 matrix* *C(1,1)=2∗Csz/Mc$C(1,3)=-Csz/Mc$C(1,5)=-Csz/Mc* *C(2,2)=2∗Csz∗lc∗∗2/Jc$C(2,3)=Csz∗lc/Jc$C(2,5)=-Csz∗lc/Jc* *C(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/Mt* *C(4,4)=2∗Cpz∗lt∗∗2/Jt$C(4,7)=Cpz∗lt/Jt$C(4,8)=-Cpz∗lt/Jt* *C(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/Mt* *C(6,6)=2∗Cpz∗lt∗∗2/Jt$C(6,9)=Cpz∗lt/Jt$C(6,10)=-Cpz∗lt/Jt* *C(7,3)=-Cpz/Mw$C(7,4)=Cpz∗lt/Mw$C(7,7)=Cpz/Mw* *C(8,3)=-Cpz/Mw$C(8,4)=-Cpz∗lt/Mw$C(8,8)=Cpz/Mw* *C(9,5)=-Cpz/Mw$C(9,6)=Cpz∗lt/Mw$C(9,9)=Cpz/Mw* *C(10,5)=-Cpz/Mw$C(10,6)=-Cpz∗lt/Mw$C(10,10)=Cpz/Mw* *! Stiffness matrix* *K(1,1)=2∗Ksz/Mc$K(1,3)=-Ksz/Mc$K(1,5)=-Ksz/Mc* *K(2,2)=2∗Ksz∗lc∗∗2/Jc$K(2,3)=Ksz∗lc/Jc$K(2,5)=-Ksz∗lc/Jc* *K(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/Mt* *K(4,4)=2∗Kpz∗lt∗∗2/Jt$K(4,7)=Kpz∗lt/Jt$K(4,8)=-Kpz∗lt/Jt* *K(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/Mt* *K(6,6)=2∗Kpz∗lt∗∗2/Jt$K(6,9)=Kpz∗lt/Jt$K(6,10)=-Kpz∗lt/Jt* *K(7,3)=-Kpz/Mw$K(7,4)=Kpz∗lt/Mw$K(7,7)=Kpz/Mw* *K(8,3)=-Kpz/Mw$K(8,4)=-Kpz∗lt/Mw$K(8,8)=Kpz/Mw* *K(9,5)=-Kpz/Mw$K(9,6)=Kpz∗lt/Mw$K(9,9)=Kpz/Mw* *K(10,5)=-Kpz/Mw$K(10,6)=-Kpz∗lt/Mw$K(10,10)=Kpz/Mw* *! Loading vector* *F(1)=Mc∗g/Mc$F(2)=0$F(3)=Mt∗g/Mt$F(4)=0$F(5)=Mt∗g/Mt$F(6)=0* *F(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))/Mw*

In 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 constants* *r,1,40e6,40e3* *! for fastener* *r,2,0.007745,0.524e-5,0.3217e-4,0.176,0.1115* *! for rail* *! Modeling of rail* *k,1,0,0,0* *k,2,0,0,Length* *k,3,0,1,0* *l,1,2* *latt,1,2,1,3,* *esize, L_rail* *lmesh,all* *! Modeling of fasteners* *∗do,i,1,Length/L_space+1* *k,10+i,0,-0.1,(i-1)∗L_space* *∗enddo* *type,2* *real,1* *l, kp(0,0,0), kp(0,-0.1,0)* *lsel,s,loc,z,0* *lesize,all,1* *lmesh,all* *lgen, Length/L_space+1,all,,L_space*

In 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

Furthermore, is written as

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 Figure 2. The wheel-rail force *P*_{w} is applied on any point between the *i*th node and the (*i* + 1)th node on the beam element.

As shown in Figure 2, 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:where and are the forces allocated on the *i*th node and the (*i* + 1)th node, respectively; and are the bending moments allocated on the *i*th 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:where *Z*_{r} is the rail displacement at the wheel location; *Z*_{i} and *Z*_{i+1} are the displacements of the *i*th node and the (*i* + 1)th node, respectively; *R*_{i} and *R*_{i+1} are the rotation angles of the *i*th node and the (*i* + 1)th node, respectively; and *S*_{1}, *S*_{2}, *S*_{3}, and *S*_{4} are the coefficients, which can be written as Hermite shape functions as follows:

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, then* *N_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)* *∗else* *N_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_rail* *L_front(i)=L_rail-L_back(i)* *∗endif* *allsel,all,all* *! Allocation of wheel-rail force* *p_front(i)=(L_back(i)+3∗L_front(i))∗p(i)∗L_back(i)∗∗2/L_rail∗∗3* *p_back(i)=(3∗L_back(i)+L_front(i))∗p(i)∗L_front(i)∗∗2/L_rail∗∗3* *M_front(i)=p(i)∗L_back(i)∗∗2∗L_front(i)/L_rail∗∗2* *M_back(i)=p(i)∗L_back(i)∗L_front(i)∗∗2/L_rail∗∗2* *f,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)* *∗enddo*

In 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 *Z*_{r} are the displacements of train system and rail, respectively. *! Elastic compression deformation between wheel and rail* *del_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,then* *p(i)=(del_z(i)/GG)∗∗1.5* *∗else* *p(i)=0* *∗endif* *F(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 Section 2.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:

Lateral motion:

Torsional motion:

In equations (18)–(20), *z*_{r}(*x*, *t*), *y*_{r}(*x*, *t*), and *ϕ*_{r}(*x*, *t*) are vertical, lateral, and torsional displacements of the rail, respectively; *m*_{r} is the rail mass per unit length; *ρ*_{r} is the rail density; *EI*_{y} and *EI*_{z} are the rail bending stiffness to *Y*-axle and *Z*-axle, respectively; *I*_{r0} is the torsional inertia of the rail; *GJ*_{t} is the rail torsional stiffness; *F*_{rVi}(*t*) and *F*_{rLi}(*t*) are the vertical and lateral dynamic forces of the *i*th fastener; *P*_{j}(*t*) and *Q*_{j}(*t*) are the *j*th wheel-rail vertical and lateral forces; *M*_{Fi}(*t*) and *M*_{wj}(*t*) are the moments applying on the rails due to forces *F*_{rVi}(*t*) and *F*_{rLi}(*t*) and due to forces *P*_{j}(*t*) and *Q*_{j}(*t*), respectively; *x*_{i} and *x*_{wj} are the locations of *i*th fastener and *j*th wheelset; *N*_{s} and *N*_{w} 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: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 *δN*_{L} and *δN*_{R} can be expressed aswhere *δ*_{L} and *δ*_{R} are the left and right wheel-rail contact angles, respectively; *ϕ*_{w} denotes the wheelset roll angle; *δZ*_{L} and *δZ*_{R} 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 *δZ*_{L} or *δZ*_{R} 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 force *F*_{x}, the lateral creep force *F*_{y}, and the spin creep torque *M*_{z} can be described bywhere *f*_{ij} represents the creep coefficients, and the definitions of the longitudinal, lateral, and spin creepages *ξ*_{x}, *ξ*_{y}, and *ξ*_{sp} are as follows:in which *V*_{w1}, *V*_{w2}, and Ω_{w3} are the longitudinal, lateral, and spin speeds of the wheel at contact point, respectively; *V*_{r1}, *V*_{r2}, 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 aswhere **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,out* *t=t+del_t* *time,t* *nsubst,1* *!*∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ *! Determination of wheelset locations* *!*∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ *∗do,i,1,10* *X(i)=X(i)+V(i)∗del_t+(0.5+bosai)∗A(i)∗del_t∗∗2-bosai∗AA(i)∗del_t∗∗2* *V(i)=V(i)+(1+fai)∗A(i)∗del_t-fai∗AA(i)∗del_t* *∗enddo* *∗do,i,1,10* *AA(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,10* *A(i)=F(i)-mid1(i)-mid2(i)* *∗enddo* *∗enddo*where *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 Figure 3. 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.

#### 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].

Subject to impact loads, such as wheel flat and welding joint, the dynamic responses of the train-track system are illustrated in Figure 4.

**(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 Table 1, 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.

##### 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 Figure 5), 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.

**(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 Figure 6, 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.

**(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 on

On 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 applied

Compared to the proposed calculation method with the classified 4 methods (methods A-D) in Section 1, the computational accuracy and calculation time of different methods are concluded in Table 4.

As 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 (Figure 7) 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.

Adopting the above parameters, a 2D finite element model of this transition section is displayed in Figure 8. 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.

**(a)**

**(b)**

The rail deformation caused by settlement is displayed in Figure 9. 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 Figure 10, 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.

Figure 11(a) shows the wheel-rail force of the 3^{rd} wheelset when a high-speed vehicle runs through the bridge-subgrade transition section. As seen from the result, when the 1^{st} 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 3^{rd} 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 3^{rd} wheelset begins to decrease. Furthermore, when the 3^{rd} 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.

#### Data Availability

No data were used to support this study.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest in preparing this article.

#### Acknowledgments

This work was supported by the Basic Natural Science and Frontier Technology Research Program of the Chongqing Municipal Science and Technology Commission (Grant numbers: cstc2018jcyjAX0271 and cstc2019jcyj-msxmX0777), the National Science Fund for Distinguished Young Scholars of China (Grant number: 51425801), the Open Research Fund of State Key Laboratory of Traction Power (Grant number: TPL1901), the China Postdoctoral Science Foundation (Grant number: 2019M650236), the Open Research Fund of MOE Key Laboratory of High-speed Railway Engineering (Grant number: 2017-HRE-01), the Open Research Fund of Chongqing Key Laboratory of Railway Vehicle System Integration and Control (Grant number: CKLURTSIC-KFKT-201804), and the Project of Fund Cultivation of Chongqing Jiaotong University (Grant number: 2018PY14).