The assessment of the integration methods for the rail vehicle ride dynamics solution

The article is devoted to the analysis of the dynamic characteristics of the mechanical system composed of rigid bodies mutually connected with flexible bindings that represent springs and dampers. The analysed system represents the rail vehicle, which is distinguished by the contact link between the wheel and rail, or wheelsets and the track. Currently, analyses of similar type are carried out basically by using commercially distributed software packages, where you enter inputs (which the input forms allow) and the user may achieve results respecting the handling procedures. In accordance with the different approaches to solutions of partial solutions of the issue, the results of the simulations can differ. In the article the theoretical analysis is conducted the parameters of calculation are set out and the results are compared with the results obtained from the calculation by SIMPACK.


Introduction
The article is devoted to the analysis of dynamic properties of multibody mechanical model representing a rail vehicle.The motivation for the creation of the article is the fact that there exists a steadily rising demand for detecting motion characteristics of transport machinery with the simulation calculation on computers [1].Compared with the standard procedures, the transport machines used to be designed, a prototype was produced, the prototype was tested and adjusted on the basis of the results of this method, has distinct advantages in lower implementation costs [2,3].The price of the computer and the software comes out for sure to be cheaper than to make a new bogie, a wagon or locomotive.And it is even faster.Such considerations do not often take into account the quality of the people that such a simulation carried out.
The quality of working staff doing the computational analyses lies in education, skills with creation of relevant computation models in a particular program environment.In addition to skills when working with the running program, it is necessary to know the issue of equipment, which are the subject of simulation analyses, their standard or specific technical parameters, terms and conditions of the external action of the environment on

Mathematic model
The procedure of mathematic model creation is processed in the following text: rail vehicle equations of motion, that describe its movement in the straight track and in track arcs.The procedure was algorithmized and written into the programme language code.This programme code represents a complete software solution.By means of this programme it is possible to gain the results of the simulated ride of the rail vehicle with specified parameters of the mechanical system of the vehicle and track.Contact ellipses and the normal stress course are calculated by means of the Hertz method [6,7,8,9].
For middle wheelset radius, middle contact angle value, track curvature and angle velocity of the wheelset revolution is valid: ; ; where: r 0 ...the average radius of the wheel [m], r l ...an immediate radius of the left wheel [m], r r ...an immediate radius of the right wheel [m], δ 0 ...middle angle of tangent planes [rad], δ l ...angle of the left wheel tangent plane [rad], δ r ...angle of the right wheel tangent plane [rad], a...half distance of the contact circles [m], C...bending curvature of the track [m -1 ] (In the case of right arc is valid the sign "+", in the case of left arc is valid the sign: " -"), R...radius of a track [m], γ 0 ... angle of attack, when the bogie in the track is in a chordal position [rad] (For the first (front) wheelset is valid the sign: " -", for pre second (rear) wheelset is valid the sign: " + "), Ω... angular speed of the wheelsets [rad.s - ], v 0 ... half distance of the contact circles [m].

Fig. 1. Representation of the rolling stock chassis passing arc track
In the Fig. 1 there is a rail vehicle schematically depicted.It moves from the left to the right side direction in the track arc with given velocity.A vehicle approached to the track stretch of rails in arc under the angle of attack.
For the relative slip vector is valid: where: s i,1 ... slip in the contact points of the left or the right wheel in the direction 1 (=x), 2 (= in the transverse direction of the tangential plane) and 3 (in the direction perpendicular to the tangential plane, [-], [-], [m -1 ],  ... the angular speed of the approaching wheel of the first wheelset during the rotation about the z-axis [rad.s - ], γ... wheelset turning [rad],  ... the angular speed of the wheel sets in rotation around the x-axis [rad.s - ].
For combined modulus of elasticity is valid: where: G... combined modulus of elasticity in shear [Pa], G W ... modulus of elasticity in wheel shear [Pa], G R ... modulus of elasticity in rail shear [Pa].
For tangential forces and moments in the contact points is valid: where: F ti,1 ... tangential force at the contact point left or right wheels, in the (lengthwise) direction of the 1 axis of the contact plane [N], F ti,2 ... tangential force at the contact point left or right wheels, in the (transverse) direction of the 2 axis of the contact plane [N], l i ... spin moment of the left or right wheel around an axis 3 perpendicular to the contact plane of the rail and wheel [N.m].
For normal forces is valid: where: f nl ... the normal force at the contact point of the left wheel [N], I 2 ... the moment of inertia of the wheel when turning around an axis of rotation [kg.m 2 ], Q l ... the wheel force on the left wheel [N], Q r ... the wheel force on the right wheel [N].
where: f nr ... the normal force in the contact point of the right wheel [N].Depiction of forces acting on a wheelset is in Fig. 2.

Fig. 2. Depiction of forces acting on a wheelset
For the left wheel force is valid: Where: Q l ... the wheel force on the left wheel [N], y r ... the coordinate of the contact point of the right wheel and rail [m], y ... transverse displacement of the wheel [m], F z pl ... the force in place of axle box (primary suspension) of the left wheel [N], F z pr ... the force in place of axle box (primary suspension) of the right wheel [N], a p ... half distance of the bearings boxes (primary suspension on one wheelset) [m], For right wheel force is valid: where: Q r ... the wheel force on the right wheel [N], y l ... the coordinate of the contact point of the left wheel and rail [m].
Contact forces and moments were computed by means of FASTSIM method.

Wheelset equations of motion
We will express equations of motion in the matrix shape (10) In addition to the forces acting on the wheelset, the force acts on each vehicle when operating during the transition as a result of the effects of gravity and centripetal acceleration.(11):

Vehicle equations of motion
The vehicle model consisting of two two-axis bogies (four wheelsets à four degrees of freedom, two bogie frames à six degrees of freedom) and one body of wagon (six degrees of freedom) has 34 degrees of freedom.In the following we presuppose, that the model consists of rigid bodies connected with flexible bindings.The rigid bodies are defined by mass matrix and the flexible bindings by stiffness matrix and damping matrix.
We can express the vehicle equation of motion in a matrix shape:

vector of external loads [-].
A vector of the external load forms the right side of a matrix equation as well as operating force.For the calculation, it is possible to use different types of methods.In the following there are depicted all stabil methods, that were tested for equation of motion solution.

Central-difference method
Initialize   0 X ,   0 X and   0 X .Select time step Δt and calculate integration constants: , damping [C] matrixes and {f} load vector: Form effective stiffness matrix: Calculate effective force vector at a time t + Δt: Solve for displacements at time t + Δt: Calculate   X and  

Fourth Order Runge-Kutta Method
One of the methods that can be used for such simulations with the advantage of the method is Runge-Kutta.To get results, we have used the Runge-Kutta methodin this case 4-th order.
Transformed into first order matrix differential equation as follows: then the equation becomes: Assuming the duration of each time step is h , ti=ih , i=1,2,…then the computational formulas over the time interval [ti,ti+1] are as follow: , where we assume: Computational formula is:
Form stiffness [K], mass [M], damping [C] matrixes and {f} load vector: Form effective stiffness matrix: effective force vector at a time t + Δt: Solve for displacements at time t + Δt: Calculate   X and  

Park Stiffly stable method
Select time step Δt and calculate integration constants: Calculate effective force vector at a time t + Δt: Solve for displacements at time t + Δt:

Computational model of a vehicle
The processing computational model is characterized by parameters of the vehicle, its geometric and mass parameters, parameters of flexible linkages, dampers.Others are parameters of the contact links and wheel track and the track parameters.All necessary parameters are set out deliberately to make them clearly referred to the input data, which were used to obtain the results using our own programme and a commercial package.

Fig. 4. Schematic drawing of the vehicle
In accordance with the specification, our model is a simple vehicle compounded from two two-wheelsets bogies and one wagon body (Fig. 4).All flexible bindings are linear, bodies are taken as rigid bodies.The data of the vehicle are introduced in the following tables from Tab. 1 to Tab. 4.

Geometric characteristics of a wheelset and a track
Geometric characteristics of wheelset and track can help to predict the movement of a wheelset in a track.At the same time, they depict the characteristics of reciprocal attitude of wheel tread profile and rail head profile.These characteristics represent an important input parameter for a mechanical system excitation due to the wheels of wheelset rolling in rails of track.Both functions input into the computation directly.From all geometric characteristics we are introducing Delta R and Tangent Gamma functions only.From Fig. 5 and Fig. 6 it is visible, that the most important curve changes are when the transversal movement of a wheel tread profile along a rail head profile reaches cca 5 mm.
Free track gauge channel is exhausted in this point and a wheel is via its flange climbing the rail head side.In the Tangent Gamma function curve this region is characterized by a steep change of course.

Track parameters
The model track consists of connected geometrically different track parts: straight line track sections, transition curves, superelevation ramps and from the track in arcs.It arises from Fig. 7, that length of the track is 3 km.Geometry of the track in horizontal plane is depicted by a bold black line.It has the vertical axis description on the right side.

Computation results
Here are compared two different physical variables for the first or third wheelset.In the following graphs we can see the course of wheelset lateral shift (displacement) in (y), slips (S) in the separate directions, guiding forces (Y), sum of guiding forces (SY), wheel forces (Q).At the same time there are evaluated and graphically depicted the forces in the contact plane of wheel and rail: normal forces (N), tangential forces (T) in special directions.In Fig. 8 there is a view at a vehicle model in SIMPACK programme.The look of the model is demonstrative only.We utilized the basic depicting tools, that are a part of the programme environment.

Fig. 8. View at a vehicle model in SIMPACK programme
The results obtained from the computations executed by means of our own software, where are algorhitmized are the equations of motions and auxiliary relations, presented in this paper.They are compared with results obtained by means of computations performed with SIMPACK programme package.
The results gained by means our own programme are depicted with red colour.The results from SIMPACK are of black colour.

Shift and rotating of the wheelset
We may take as the same the wheelset shift in the y-direction obtained by both of procedures (Fig. 9).

Fig. 9. Shift of the wheelset in the y-axis direction
We may take very similar a revolution of the wheelset around z-axis by both of procedures.Small differences appear at transitions from one arc to another arc (Fig. 10).Slip parameters are important for tangential forces calculation.This is the reason, that the "overshootings" during the transition of curves courses that are visible at slips in the xdirections of right (Fig. 11) and left wheels (Fig. 12), appear in the tangential forces calculation (Fig. 17 and Fig. 20) too.

Normal forces
Normal forces in their results represent the force acting in the normal axis direction to the contact plane in the rail/wheel contact point.The course of their values is very similar again (Fig. 15 and Fig. 16).

Tangential forces
When we evaluated the tangential forces, the differences appear in the places, where vehicles come from one track geometry into the other track geometry (straight, a ramp, an arc).The course of the tangential forces in the longitudinal x-direction gained by means of SIMPACK programme is of more steep inclination of the transitions (Fig. 17 -Fig.20).This course comes back to the stable state when a vehicle runs in a straight line track or in an arc via "overshooting".In this cases, there are short-timed detected higher values of forces.The "overshoootings" forces courses follow from the "overshooting" courses of slips.

Guiding forces
The guiding forces course is an important parameter for results comparison (Fig. 21 -Fig.23).This parameter gives not only evidence of lateral force influence of the wheel to the rail, but in this case it gives the view at the result of processing of the contact parameters of the wheelset and track, contact patches and stresses and assumptions of the force acting between a vehicle and a track.The close course of the individual forces verified results obtained with the mathematical expression investigated the movement state of a vehicle on the basis of the results obtained from the programme, which is now widely used by reputable professionals.Or when viewed from the opposite side: the results produced by an automated computing system with a closed core and prescribed (and also some optional) calculation procedures, we would be able to get the solution to the system of equations derived in this article.

Wheel forces
The computations results of the wheel forces at the left and the right wheels are very similar to each other, see Fig. 24 and Fig. 25.We tested the methods in the way, that all model bodies were moved (deflected from centered position) 0.003 m in transversal y-direction (Tab.5).The graphical depiction of wagon body position deflection in y-axis direction is shown in Figs. 26, 27, 28, 29, 30.The difference is in computation with various damping parameters (γ).The Integration methods stability comparison is in Table 5.The parameters for computations in accordance to equations (47) and (48) when γ = 0 are in the Table 6.The parameters for computations in accordance to equations ( 47) and (48) when γ = 0.15 are in the Table 7.
Table 7. Wagon body position deflection in y-axis direction and calculated damping parameters when = 0.15 The parameters for computations in accordance to equations ( 47) and (48) when γ = 0.2 are in the Table 8.The parameters for computations in accordance to equations (47) and (48) when γ = 0.3 are in the Table 9.

Conclusion
At present, the analyses of the dynamic characteristics of railway vehicles are carried out basically by using commercially distributed software packages, where you enter inputs (which the input forms allow) and the user may achieve results respecting the handling procedures.In view of the different approaches to solutions of partial solutions, the results of the simulations can differ.
In the paper the theoretical analysis of issue solution is performed, the computation parameters are specified and the results of computations performed with own programme are compared with the results gained from programme SIMPACK.The evaluated results are for the first or the third wheelset.In the figures with graphical representation of the courses: wheelset shift in the lateral direction (y), wheel-set rotation, slips (S), forces in the contact plane of rail and wheel: normal (N), tangential (T), guiding forces (Y), sum of guiding forces (SY), wheel forces (Q) there can be seen a very good coincidence at the computation with programme SIMPACK and our programme, here presented theory is processed in this.We analysed various time integration methods for judgement of the suitability, stability and numerical damping of these methods.
Finally, we can state, that the calculations met the goal, that we had set to ourselves.The results obtained by means of commercial package in general confirm the rightness of the theoretical analysis and realized algorithm in our own programme code.At the same time the theories in the investigated issue confirm the trustworthiness of the results obtained by commercial package.In the following period it can enable to us to model the flexible bindings and dampers with characteristics, that program SIMPACK does not directly enables and to use other contact theories of the rail and wheel contact with the nonelliptic contact patch.
: U ... cant of the outer rail in the track arc [m], 2.a ... distance of contact circles [m].The accelerations distribution is depicted in Fig.3

Fig. 5 .
Fig. 5. View at the dependence of the Delta R function on the lateral shift of the wheel profile after rail profile.

Fig. 6 .
Fig. 6.View at the dependence of the Tangent Gamma function on the lateral shift of the wheel profile after rail profile.

Fig. 7 .
Fig. 7. Track definitionThe track consists of four arcs with the radius of R 293 meters, three of them are left and three of them are right.The axis is on the left side.In the bottom part of graph there are depicted the superelevations of the outside rail 0.144 m.The vertical axis for this graph is on the right side.This track has neither climbing nor descending.For the computation transition curves and superelevation ramps of the Bloss were used.

Fig. 11 .
Fig. 11.Slip in the x-direction on the right wheel of the first wheelset

Fig. 12 .Fig. 13 .
Fig. 12. Slip in the x-direction on the left wheel of the first wheelset

Fig. 14 .
Fig. 14.Slip in the y-direction on the left wheel of the first wheelset

Fig. 15 .Fig. 16 .
Fig. 15.Normal force on the right wheel of the first wheelset

Fig. 24 .
Fig. 24.Wheel force on right wheel of the first wheelset

Fig. 30 .
Fig. 30.Wagon body position deflection in y-axis direction calculated by HHT method when γ = 0.3 The work was supported by the Cultural and Educational Grant Agency of the Ministry of Education of the Slovak Republic in project No. KEGA 077ŽU-4/2017: Modernization of the Vehicles and engines study program.The work was also supported by the project No. APVV-0842-11: Equivalent railway operation load simulator on the roller rig and VEGA No. 1/0927/15: Research of the use of alternative fuels and hybrid drives on traction vehicles with aim to reduce fuel consumption and air pollutants production.Research-Educational Centre of Rail Vehicles (VVCKV) Conferences 157, 03013 (2018) https://doi.org/10.1051/matecconf/201815703013MMS 2017

Table 1 .
Mass vehicle parameters

Table 2 .
Flexible bindings parameters

Table 4 .
Geometrical parameters of the vehicle and track

Table 5 .
Integration methods stability comparison

Table 8 .
Wagon body position deflection in y-axis direction and calculated damping parameters when = 0.2

Table 9 .
Wagon body position deflection in y-axis direction and calculated damping parameters when γ = 0.3