An Alternative Solution of Train-Track Dynamic Interaction

,


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][3][4][5].As a matter of fact, running trains cause track vibrations, which in turn aggravate the train vibrations.erefore, 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][7][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.ese 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, MAT-LAB, and C, which are then solved with numerical integration methods. is 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][9][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].
is 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 SIM-PACK and UM, while the track structures are simulated using finite element software.en, the modal properties of tracks are exported from finite element software into multibody software to make the two different systems coupled [15,16]. is 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.
en, 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.
rough the above methods, several generalized conclusions are reached by the authors: (i) e 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.

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][7][8][9][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.e 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.

2D Train Subsystem.
e 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 e dynamic equations for all the seven bodies are given as follows.
Vertical motion of car body: Pitch motion of car body: (2) 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: e 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 ese matrices of the train system cannot be calculated by ANSYS automatically, which should be written into the software manually.e 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 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.

2D Track Subsystem.
e 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 Timoshenkobeam element BEAM188, the fastener systems are usually simulated by spring-damper element COMBIN14, and the slabs are established adopting 3D solid element 4 Shock and Vibration SOLID185.e 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

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.

2D Train-Track Interaction.
e 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].erefore, the wheel-rail force is determined by Furthermore, δZ(t) 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.
e 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.erefore, it is necessary to allocate the wheel-rail force to the two adjacent nodes, as shown in Figure 2. e wheel-rail force P w is applied on any point between the ith 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 p i w and p i+1 w are the forces allocated on the ith node and the (i + 1)th node, respectively; M i w and M i+1 w 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: where Z r is the rail displacement at the wheel location; Z i and Z i+1 are the displacements of the ith node and the (i + 1)th node, respectively; R i and R i+1 are the rotation angles of the ith 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:

Shock and Vibration 5
As a result, the rail displacement at any point can be obtained combining equations ( 16) and ( 17). e allocation of wheel-rail force in APDL language is written as * do,i,1,4 !Determination of rail nodes applied by wheel-rail force 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.

3D Train-Track Dynamic System.
e 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.
e 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.e detailed equations of motion of all the seven parts can be referred to literature [22,23].
e rails are modeled as Euler beams and discretely supported by fasteners which are simulated as linear springdamping elements.ree 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: 6

Shock and Vibration
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 ith fastener; P j (t) and Q j (t) are the jth 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 ith fastener and jth 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 wheelrail normal contact forces: where G is the wheel-rail contact constant and δN(t) is the wheel-rail normal elastic compression deformation.e left and right wheel-rail normal elastic compression deformation δN L and δN R can be expressed as where δ 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 by where 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.
e wheel-rail interaction forces mainly include wheelrail 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 traintrack dynamic model is not listed below, which can be referred in the first author's doctor degree dissertation [24].

Numerical Integration Method.
rough the above explanation, the whole dynamic equations of the train-track system are established, which can be solved by the step-bystep 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.e 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 Shock and Vibration 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.
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.is 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.

Procedure of the Proposed Method.
e 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.

Validation of the Proposed
Calculation Method e calculation method is introduced in detail, and on this basis, two kinds of validations are conducted to verify its validity.

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

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.

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 8 Shock and Vibration 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.e 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].
e 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.

Advantages and Disadvantages of the Proposed Calculation Method
Every calculation method certainly has advantages and disadvantages, which will be discussed in this section.e advantages are given as follows: (i) e 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) e 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) e calculation method can also be used in investigating train-bridge interaction, train-trackbridge interaction, train-tunnel interaction, and so on      10

Shock and Vibration
On the contrary, the disadvantages of the method are listed below: (i) e computational efficiency is lower than that programmed with FORTRAN, C and so on (ii) e 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.

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. e dynamic parameters of the train can be found in literature Shock and Vibration 11 [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.
e 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.e deformation curve is not a straight line but an irregular curve, indicating that simplifying this rail deformation as straight     Shock and Vibration 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 bridgesubgrade 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.en, 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.

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) e proposed calculation method is effective to investigate TTI. e whole train-track system is coupled and calculated synchronously.( 2) is paper also provides a guidance to calculate a coupled mechanical system with different integration methods in ANSYS.(3) is calculation method has a wide application range, which can also be used to study train-trackbridge interaction, train-tunnel interaction, and so on.(4) e 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.

Figure 2 :
Figure 2: Allocation of wheel-rail force on the beam element of rail.

Figure 4 :
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.

Figure 5 :
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.

Figure 6 :
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.

Figure 9 :
Figure 9: Deformation of rail caused by settlement.

Figure 10 :
Figure 10: Structural dynamic stress of the infrastructure at a certain time.

Figure 11 :
Figure 11: Dynamic responses of the train-transition section system: (a) wheel-rail force and (b) dynamic stress in subgrade bottom layer.

Table 1 :
Main notations adopted in this work.

Table 2 :
Parameters of the adopted vehicle.

Table 3 :
Calculated responses of CRH3 by UM and the proposed method.

Table 5 :
Key parameters of the bridge-subgrade transition section.

Table 4 :
Evaluations of computational accuracy and calculation time of different methods.