## Abstract

Train safety and operational efficiency can be improved by investigating the dynamics of the train under varying conditions. Longitudinal train dynamics (LTD) simulations performed for such purposes, usually by utilising a nonlinear time-domain model. This paper covers two modes of LTD results corresponding to the time domain and frequency domain analysis. Time-domain solutions are essential to evaluate the full response used for parameter optimisation and controller design studies while frequency domain solutions can provide significant but straightforward clues regarding system dynamics. An advanced draft gear model, which works under a four-stage process is constructed considering all structural components, geometric relationships, friction modelling and dynamic characteristics such as hysteresis, stiffening, state transition, locked unloading, softening. Then, this model is parametrically reduced and implemented into an LTD simulation. The simulation in the time domain is conducted assuming a locomotive connected with a nine wagon via “ode3” fixed-step solver. The transfer function among the first wagon acceleration (output) and the locomotive force (input) estimated through system identification methodology. Then, the identification results interpreted by investigating step-response characteristic and best response giving parameter set is selected. Next, the modal and spectral analysis performed to reveal the behaviour of the in-train forces and the effects of vibration. This paper discusses a reliable methodology for the longitudinal dynamic analysis of the multi-bodied train in time and frequency domain and clarifies in-train vibration behaviour under the existence of sophisticated draft gear.

## 1 Introduction

The study of longitudinal train dynamics (LTD) is one of the principal analysis performed in the design of all kinds of railway vehicles (freight car, high-speed train, locomotive, EMU, DEMU, rolling stock, heavy-haul or passenger train etc.) [1]. The investigation of the wagon’s motion in the direction of the railroad is the main subject in LTD simulation covering the overall movement of the train and the relative movements resulting from the looseness of the connections between the rolling stocks. Reducing longitudinal vibrations is indispensable to enhance comfort, sustain stability and avoid derailment. Thus, LTD calculations have often been made by engineers for such purposes. Also, LTD analysis needed with the aim of designing useful simulation models where locomotive performance can assess. Such kind of models primarily includes locomotive and wagon parameters, track characteristics, and train handlings (throttle, brake, dynamic, etc.) as input parameters.

The response of nonlinear draft gear and end-of-car cushioning units in impact are also critical issues, which should be considered in the LTD simulations. Draft gears are highly nonlinear spring-damper sub-systems in which include Coulomb damping. Moreover; if two wagons are coupled, their nonlinearities increase due to a dead band existing among couplers. As such nonlinearities make mathematical models comparatively complicated, therefore solutions in the time domain are preferred. However; frequency domain analysis, which can provide quick and understandable results, can also be feasible under some simplifications and assumptions.

The main characteristics of the LTD classified as problem genre (single degree of freedom per vehicle, nonlinear rigid multi-body system dynamics, etc.), force inputs (coupler force, traction force, dynamic brake force, air brake force, curving resistance, propulsion resistance, and gravitational component etc.) and main applications (evaluation of coupler forces, train speed, and train position, traction and adhesion control, coupler system and brake system design etc.). The effects related to LTD comprise the following elements: steady forces (traction or braking, coupler forces and curving, propulsion, gravitational resistances), impact forces and low-frequency vibration forces [2]. The locomotive produces traction and brake forces and represents the primary inputs, which cause the dynamic behaviour of the train. It is necessary to perform LTD analysis to prove that the train operating parameters, load distributions, connection elements are selected optimally and that the train guarantees excellent dynamic performance.

While the trains are getting more prolonged and more massive, the computational complexity increases, therefore, advanced simulations with more efficient calculation methods and massive hardware should be conducted [3]. Numerical solver methods and their step size also discussed in the LTD simulations. Runge-Kutta, Euler, Adams-Bashforth methods have been widely implemented to solve LTD problem. In LTD simulations, wagon-wagon or wagon-locomotive connections are usually provided integrating pair coupled draft gear or buffer models. The draft gears have a nonlinear and complex nature, which is generalised as look-up tables with various patterns. They can be classified as friction, rubber, hydraulic and types of their combinations.

Forces in draft gears are related to draft gear deflections. Angles at the connections determine how the wagon body and running gear react to these forces. Therefore, a detailed mathematical model becomes essential for examining the force-deflection characteristics of the draft gears. Such a model should contain a force-displacement hysteresis loop to realise Coloumb energy absorption [4]. The transition from kinetic to static friction included in a dynamic model, which is developed to simulate the draft gear behaviour in rail-car impacts during the impact. The complexity of the draft gear model can be chosen according to the accuracy required for the simulation. However, there exists a trade-off between efficiency and computational effort or consumed time.

The advanced models allow modelling a broader spectrum of dynamic behaviour and include significant features in designing a draft gear [5]. LTD simulations perform high accuracy with developing computing availabilities. Advanced component modelling and vehicle connection mechanisms serve a better solution with the aid of parallel computing.

Advanced friction draft gear dynamics should be modelled due to the occurrence of larger in-train forces and more complex force patterns to assess the LTD for longer, more substantial and faster heavy haul trains. An auto-coupler having friction draft gear is a mechanism which connects wagons and locomotive. Drawbars can be replaced with auto-couplers to reduce the coupling slack and to lower the dynamic in-train forces. Friction draft gears consist of springs, wedges and plates enclosed by a housing. Their functions are transmitting the longitudinal forces and damping longitudinal in-train forces like friction dampers, respectively. A review for dynamic modelling of friction draft gear has been written by the authors [6, 7].

If a railway vehicle is accelerated or decelerated along the curved track, a coupling effect occurs between the LTD and the lateral train dynamic. Due to the excessive lateral force acting on the wheels, the derailment is expected with the increase of the derailment index. The breaking wave propagation time increases proportionally to the train’s length. When a free rotation angle of coupler occurs, the coupler system deflects from the centre line, and the vehicle body is subjected to a horizontal force at the coupler position. The pressing force of the coupler may grow large enough leading to locomotive derailment [8].

A sequential, parallel and hybrid schemes were proposed to simulate the LTD and compared concerning their features of the computing scheme in a distributed power train (2 locomotives+ 105 wagons+ 2 locomotives + 105 wagons). A parallel computing scheme, which executes multiple processes simultaneously, outperforms the others regarding improving computational efficiency [9]. In [10], the feasibility of the parallel computing scheme has represented by a simulation of a long heavy haul train having 214 vehicles (2 locomotives + 105 wagons + 2 locomotives + 105 wagons). The draft gear types are friction draft gear between wagon and locomotives, polymer draft gear between locomotive and locomotive, respectively. A 3-hour train trip simulated with 216 cores, whose accumulated computing time of all cores was about 253 days, while the wall-clock time was about 29 hour.

A computer program to simulate the LTD has designed including distributed/remote power, electronically controlled pneumatic brakes, new brake valves, M-901G draft gears and cushion units with an active draft stroke, etc. [11]. By utilising a measured functional relationship for draft gear characteristics, a nonlinear LTD model was developed to estimate the coupler forces [12]. The reliability of a coupled train/track model verified through a commercial multi-body analysis software “SIMPACK” and dynamical performance, the vibration frequency of the components, ride comfort and curving behaviour of the model were discussed [13]. A map for all types of in-train forces in long trains (steady, impact, sustained longitudinal oscillation, etc.) was created by using a wireless measurement system in real-time [14]. A simulation model, which covers power generation (internal combustion engine, electric motors and generator) and transmission in the diesel-electric locomotive components in detail, was built to calculate the fuel consumption. The locomotive performance and operational conditions in each throttle position were tested experimentally in dynamometer from the viewpoint of efficiency maps [15]. In [16], a non-inertial longitudinal train force model with a nonlinear coupler developed.

This model allows three-dimensional motions of the vehicle body and captures kinematic degrees of freedom that are not caught using existing simpler models.

In this part of the introduction, studies conducted in the frequency domain examined. The longitudinal vibration of the train consists of the lower frequency, which expands and contract of an overall railroad and the higher frequency part, which causes shock impact. Low-frequency oscillations were classified into two different modes, namely cyclic vibration and sustained longitudinal wave [17]. Continued longitudinal wave, which was underdamped, is arisen if the whole railway vehicle remains in a single stress state (tensile or compression). Cyclic fluctuations correspond to the square wave oscillations originated from run-in/run-out behaviour. Draft gear’s in-train characteristics affect the vibrational behaviour of the train. The LTD simulations and frequency response simulations conducted in [18] revealed that three types of draft gears gave similar coupler forces but different draft gear deflections and impact characteristics at different frequencies. The mode shapes and rates for an individual draft gear were found utilising the finite element approach in [19]. Which draft pad should be checked for fatigue response under the cyclic loading of 10 Hz, was clarified within a range of frequency values. This work supplies a tool for detecting new frequencies of draft gear, whose modal frequencies are influenced by the crack due to parameters like aspect ratio, crack width, and break orientation [20].

Although the frequency domain solutions are not practical due to nonlinear system parameters, this does not necessarily mean that frequency domain investigations are useless. Indeed, they can quickly provide the system response parameters such as natural frequencies in the lower and upper limits so that engineers have more insight while choosing appropriate time step values and analysing time-domain results.

As a summary, LTD has been studied by many researchers, most of whom have dealt with the solution methodology for various design needs and objectives. However, there is still a need for a robust and compact method, including system identification approaches which can provide sufficiently accurate results in a short time to enhance the controller design performance.

In this study, the dynamic simulation of a locomotive connected to nine wagons through a nonlinear draft gear was performed in “Matlab&Simulink”. A speed input trajectory pattern fed into the system and the acceleration of the locomotive was controlled via a discrete-PID controller to minimise the error between the reference speed and actual speed of the train. The simulation scenario includes traction and braking phases while turning in a sharply curved track and climbing a slope with a full load condition. After the simulation phase, which uses full-length track information (curvatures, gradients and speed limits), train information and control information, the results were analysed both in a time and frequency domain. The time-domain analysis comprises the system identification process, which derives the discrete transfer functions among the first wagon acceleration as output and traction force as input under the conditions of the various coupler parameter space and the step response characteristics of each parameter set. The simulation results of the coupling parameter set, which demonstrates the best time response behaviour in terms of the system dynamics, were examined and the wagon couplers, which are subjected to the most critical transitions were determined by considering the forces exerted on the wagon couplers in four different track/train phases (acceleration, curving, slope climbing, deceleration).

The frequency analysis starts with collecting natural frequencies in the longitudinal direction for two different coupler behaviour. Then, the previously obtained time-domain results transformed into the frequency domain, and the meaningful inferences about spectral effects for each track phases are interpreted.

Thanks to the reduced model built in the paper, an equally complex system can be simulated with a small computational load, high computational efficiency and accuracy. The other contribution related to the article is that the results of the LTD simulation were interpreted from the frequency domain perspective.

The rest of the paper organised as follows. Section 2 provides a mathematical modelling of the LTD. The model reduction process and the information about the traction curve are described in Section 3 and Section 4, respectively. Results obtained from the LTD simulation presented in Section 5. Section 6 represents the limitation and discussion section. Lastly, the outcomes of work were discussed and concluded in Section 7.

## 2 Mathematical modelling of the LTD

### 2.1 General description of the LTD model

This section reveals the equations of motion related to the LTD. The LTD takes account of the movement of all railway vehicles along the way including locomotives and carriages. The longitudinal behaviour of the trains is a function of the characteristics of the train control inputs coming from the locomotive, train brake inputs, track topography, track curvature and wagon connections. The longitudinal dynamic behaviour of the train is defined as a system of differential equations. For this purpose, mathematical equations are established, and simulations are conducted under the assumption that wagons do not exhibit any movement in the lateral and vertical directions [21].

In Figure 1, the schematic representation of the generalised train model is given.

The stiffness and damping coefficients of the wagon couplings (draft gear) are expressed regarding complex dynamical characteristic functions

### 2.2 Friction draft gear model

Friction draft gears have friction and velocity-dependent stiffness nature. Nonlinear hysteresis, which leads to discontinuities among loading and unloading force-displacement curves, occurs due to the friction damping. Velocity-dependent friction, slack, limiting stiffness, pre-load and translational characteristics should be considered in a draft gear dynamics model. A wedge-spring approach can be used to model the force-displacement characteristics of friction draft gears. The friction draft gears, which are utilised in LTD simulation and depicted in Figure 2, share the same structure as in [5].

The working procedure of this kind of draft gear can be separated into four different stages, which express different force-displacement behaviour. These stages are named as loading stage 1, loading stage 2, unloading stage 1, unloading stage 2, respectively. The term “loading stage” corresponds to the compression of the draft gear, while “unloading stage” means that the draft gear is being released. When the draft gear is entirely compressed, a noticeable clearance among the spring seat and the bottoms of moveable plates remains. Such kind of clearance also exists in a wholly released condition among the follower and the tops of movable plates. The values of these clearances can be predicted by the structural design and pre-load characteristic of the draft gear.

In loading stage 1, the loading process is initialised, but the follower is not in contact with the movable plates. During this stage, the elements, which are related to the friction draft gear forces, are the central wedge, wedge shoes, the spring seat and springs, respectively. The wedge group is responsible for the energy-dissipation mechanism. In loading stage 2, the central wedge and movable plates are pushed by the follower. The forces exerted on the movable plates, and a central wedge is the components of the draft gear force. Both friction elements are responsible for the draft gear force and energy dissipation. In the unloading stage 1, the unloading process begins, but the spring seat remains untouched by the movable plates. The plate group does not affect the draft gear force. In the unloading stage 2, both the spring seat and the movable plates forced by the mainspring. All friction elements are dissipating energy and contributing the draft gear force having an advanced model which contains the geometrical and dynamical characteristics of all draft gear sub-systems (friction mechanisms, springs, and housing).

This model was constructed by conducting quasi-static force analyses and kinematic relationships for all internal components. The force-displacement characteristics of the draft gear model separated into four stages, which are represented in [5] and given as in Eq. (1).

where

where

corresponds to the four different stages (loading_1, loading_2, unloading_1, unloading_2) as explained above; *F*_{fdg,i} is the external force applied on the central wedge; *F _{sm}* and

*F*are the mainspring force and release spring force, respectively;

_{sr}*α*,

*β*and

*y*are wedge shoe angles as shown in Figure 2;

*μ*

_{1},

*μ*

_{2}and

*μ*

_{3}are the coefficients of friction, as shown in Figure 2;

*μ*

_{4}refers to the friction coefficient of the moveable plates;

*k*and

_{m}*k*are mainspring stiffness and release spring stiffness, respectively;

_{r}*x*

_{m0}and

*x*

_{r0}are the pre-deflection of the mainspring and the release spring, respectively;

*x*is the friction draft gear deflection;

_{fdg}*x*is the release spring deflection.

_{rdg}The exponential friction model integrated into the force-displacement characteristic presented in this study. In the exponential friction model, the coefficient of friction is dependent on the velocity. At the same time, parameters such as displacement, acceleration, roughness, hardness and normal stress can also be included in the friction models to obtain more accurate results. A lubrication effect must also be added to the friction model to prevent detrimental stick of the wedge shoes. The friction model used for the lubricated contact surfaces of each wedge shoe described as in Eq. (2).

where *v*_{1} is the relative velocity between the wedge shoe and stationary inner plate, *v*_{f0} is the initial impact velocity of the two adjacent vehicles; *S*(.) is the modification factor given by utilising a two-dimensional look-up table, which is the function of impact velocity and deflection; *h*_{1}, *h*_{2} and *h*_{3} are specific model parameters, which reflect kinetic friction, part of static friction; the transition from static to kinetic friction, respectively.

The friction model corresponding to the all non-lubricated coefficients of friction can be calculated as in Eq. (3).

where *v _{j}* are the relative velocity components between the central wedge and wedge shoe, spring seat and wedge shoe, movable plate and a stationary inner plate, respectively.

The relative velocities among various elements must be known to calculate the friction model. The interactions between the component velocities obtained by evaluating the kinematic equations originated from the geometrical relationship as in Eq. (4).

where *v _{fdg}* is the friction draft gear working velocity.

The advanced model presented in this paper utilises a transitional characteristic named as “locked stiffness”, which first recorded in [22]. The forces leading longitudinal vibrations of the train discovered as a third type of in-train force. The locked stiffness is defined as the longitudinal stiffness of the whole train when train oscillates longitudinally. The locked stiffness cannot be captured from the impact tests or drop-hammer tests. The locked stiffness can be stated mathematically as in Eq. (5)-(7).

where *x _{t}*,

*x*

_{t−Δt},

*v*,

_{t}*v*

_{t−Δt}are draft gear deflections and relative velocities of adjacent vehicles of the current and previous time-step, respectively.

*f*and

_{l}*f*are force-displacement characteristics for loading and unloading processes, respectively.

_{u}*k*is the locked stiffness. In every time-step, the values in Eq. (5) and (6) must be evaluated, then, Eq. (7) determines which result will be produced as the final draft gear force.

_{locked}### 2.3 Train modelling

The dynamic equations of the generalised nonlinear longitudinal train model with nonlinear friction draft gear wagon coupling are stated as in Eq. 8 [23].

where *ẍ _{i}*,

*ẋ*,

_{i}*x*correspond to the acceleration

_{i}*m*) of each element, respectively;

*m*

_{wagon,i}is the mass (

*kg*);

*F*

_{r,i}is the resultant resistance force (

*kN*) of the

*i*. element;

*b*is the

_{i}*i*. damping coefficient

*m*

_{dyn,loco}is the dynamic mass of the locomotive at zero slip state (

*kg*);

*kN*);

*F*

_{r,loco}is the resultant resistance force of the locomotive (

*kN*).

## 3 Model reduction

Traction motor model used in the simulations is 3-phase squirrel-cage induction motor (2180 V DC bus voltage, 850 kW, 270/310 A, 1283/2484 rpm. The motors are controlled using field-oriented control with the indirect vector control scheme. The details in the traction motor modelling are given in our previous paper [24].

Air brake has a brake cylinder force which is converted employing rigging factors and shoe friction coefficients into a retardation force. Pneumatic control is utilised in the electronically controlled brake system. The brake designs vary pipe exhausts systems. The pneumatic control regulations done to the brakes via the brake pipe take time to propagate along the train. As the control is via a pressure wave, the system is constrained to sonic speed. The main problem in LTD simulation is to determine the braking force and the operation timing. The solution methods include a look-up table of brake forces against time, empirical models and fluid dynamics-based models [25]. The train air brake system model used in the simulation is adapted from [26].

The traction/dynamic brake force term

The highly nonlinear coupler dynamics can be analytically transformed into a simplified but also nonlinear model given in Eq. (9) [27].

where *P _{i}* is the preload,

*b*is the damping ratio

_{i}*k*

_{i3}is the linear stiffness

*k*

_{i2}is the nonlinear stiffness

*ν*is the critical velocity constant,

*k*

_{i1}is the recoil stiffness

*Δx*,

_{i}*Δx*

^{.i}are relative displacement and velocity of the wagons,

*ϵ*is the coupling clearance.

The terms *k*_{i1} and *k*_{i3} are the main cause of hysteresis behaviour. *k*_{i2} moreover, *ν* are responsible for the force peaking that can be derived both from experimental and theoretical analysis of draft gears [28].

The dynamic equations of the generalised nonlinear longitudinal train model with simplified wagon coupling stated as in Eq. (10).

The resultant resistance forces depend on the braking regulations, speed, track curvature and wagon/locomotive design parameters. The resultant resistance force of the locomotive and wagons defined as in Eq. (11).

where *F*_{p,loco} and *F*_{p,i} are the propulsion resistances which are represented as the sum of the aerodynamic and rolling resistances (*kN*); *F*_{crv,loco} and *F*_{crv,i} are the resistance force elements originating from the track curvature (*kN*); *F*_{grv,loco} and *F*_{grv,i} are the gravitational force elements due to the track gradient (*kN*); *F*_{brake,loco} and *F*_{brake,i} are the braking resistance forces caused by pneumatic braking (*kN*), respectively.

A quadratic formula has been utilised to approximate the propulsion resistance in Eq. (12) [29].

where *A*(*N*), *B**C*

*A* and *B* are related to the mass, and mechanical resistance and *C*account for air resistance.

Gradient force can be given as in Eq. (13).

where *m _{dyn}* is the dynamic mass of the vehicle system (

*kg*);

*g*is the acceleration due to gravity

*κ*is the gradient in the form scaled as ”1:

*κ*” (for example a grade of 5% is calculated as

*κ*=

where *h* is the dimensionless parameter, which changes from 500 to 1200 depending upon the railway vehicle; *R _{y}*is the radius of the curved track (

*m*).

The complex model has seven-dimensional parameter space such as [*h*_{1}, *h*_{2}, *h*_{3}, *k _{locked}*,

*k*,

_{m}*k*,

_{r}*S*] and four different stages. The simplified model has five-dimensional parameter space such as [

*b*,

*k*

_{1},

*k*

_{2},

*k*

_{3},

*ν*] moreover, two different stages.

To validate the reducibility of the complex dynamical model into a simplified analytical model, the parameters of the complex model and simplified model matched by conducting a simulation. “Simulink Design Optimization Toolbox” was utilised in this experiment. Single friction draft gear corresponding to the sophisticated and simplified model was excited during 14.6 sec., with a sinusoidal input having various amplitude and frequency patterns. The frequency and amplitude components of the excitation signal given as; 10^{6} N, ^{6} N, *π* rad/sec along 4 seconds, 10^{6} N, 2**π* rad/sec along 2 seconds, 10^{7} N, 10**π* rad/sec along 0.4 seconds, 10^{7}N, 20**π* rad/sec along 0.2 seconds, respectively.

First of all, the friction draft gear forces of the complex model were evaluated as the target outputs of the simplified model w.r.t to the parameters given above. Then the parameters of the simplified model were sought during the optimisation process. The optimisation options selected as, a genetic algorithm using a robust cost with 0.0001 parameter tolerance and 0.0001 function tolerance, sum squared error cost function, 1000 max.iteration. The simulation runs on parallel pools for practical computing.

The lowest and the highest values in the parameter space of the complex model and simplified model after the optimisation process are listed in Table 1.

Complex model |
---|

h_{1}(low − L,0.05; high − H, 0.0525) |

h_{2}(low − L, 0.01; high − H, 0.0105) |

h_{3}(low − L, 0.4; high − H, 0.42) |

S(low − L, 0.03; high − H, 0.0315) |

Simplified model |

ν(low − L, 0.1018; high − H, 0.7974) |

The output of the simplified model and sophisticated model for low and high parameter set are represented in Figure 3a and 3b, respectively.

## 4 Traction curve

The traction curve defines the evolution of the traction forces, which should be supplied for each speed value of each vehicle (tram, subway, EMU, DEMU, etc.) that produces the traction. Two fundamental values determine the traction curve;

Maximum adhesion force

Maximum traction force

These two constraints limit the traction force for different speed values of the vehicle. The limiting factor for the traction force is the maximum adhesion force. The maximum adhesion force defined as in Eq. (15).

where *ẋ*_{1} is the locomotive speed *m _{stc}* =

*n*.

*m*

_{wagon,i}+

*m*

_{stc,loco}is the static mass of the locomotive&wagon system (

*kg*).

The coefficient of adhesion is determined by several parameters; namely the material properties, rail surface, wheel profile, axle pressures and vertical movements. According to UIC standards, the maximum value of the coefficient of adhesion is 0.33. It is also shown by experimental methods that the ratio of adhesion is dependent on the speed. For this reason, the “Curtius-Kniffler” formulation, which varies with speed, is chosen to be utilised for determining the adhesion coefficient. The maximum traction power (*P*_{τ, max}*in**kW*)is defined by using a maximum traction force (*F*_{τ, max} = *m _{dyn}*.

*a*+

_{max}*F*), without considering the curving and gravitational resistances except for acceleration and braking, and the speed at which the vehicle reaches its maximum power

_{r}Maximum traction power calculated as in Eq. (16).

where *a _{max}* is the maximum available acceleration

## 5 LTD Simulation

LTD simulations were conducted for a locomotive with nine wagons. “ode3” (Bogacki-Shampine) with a fixed-step size (fundamental sample time) 0.0001 s (0.1 ms) was selected to evaluate the equations of motion-related with LTD via “Matlab&Simulink”. Transient responses of the LTD were analysed between 0-270 sec.

The computing time for all simulation cases lasts approximately averaged 854.7 seconds. The data precision was selected as double to maintain accuracy, and the truncation errors are prevented.

The locomotive-wagon system operates in a sharply curved track (*R _{y}* = 300

*m*) from

*t*= 75

*s*. to = 185

*s*. At

*t*= 155 − 225

*s*, the system is climbing %ffl 5 slope with a full load. The operational speed at maximum motor power (

*v*) was considered as 16.72

_{Pmax}*a*) moreover, deceleration are 0.976

_{max}The simulations performed for 16 cases in the parameter space including the permutations of the damping ratio, recoil stiffness, nonlinear stiffness and linear stiffness, respectively to understand the effects of the stiffness and damping parameters concerning the time and frequency response of the locomotive-wagon system. The settings and the numerical values used in the simulation represented in APPENDIX-I.

The desired acceleration of the locomotive controlled via a parallel form of discrete-time PID controller (*P* + *N*) is selected as 100 and the sample time is *T _{s}* = 1

*e*− 4. The parameters of the PID controller and the root mean square (RMS) values of the controller error for all 16 cases are given in Table 2.

Stiffness and damping parameter space | PID Controller parameters | Control results | ||
---|---|---|---|---|

Permutation | Proportional | Derivative | Integral | Controller Error (RMS) |

LLLL | 12,292 | 0,0146 | 24,5850 | 0,0634 |

LLLH | 34,300 | 0,0147 | 16,8600 | 0,0870 |

LLHL | 33,080 | 0,0125 | 16,6151 | 0,0819 |

LLHH | 29,230 | 0,0146 | 15,8453 | 0,1002 |

LHLL | 31,220 | 0,0162 | 16,2440 | 0,0624 |

LHLH | 47,720 | 0,0175 | 19,5436 | 0,0633 |

LHHL | 31,220 | 0,0175 | 16,2440 | 0,0709 |

LHHH | 22,730 | 0,0160 | 14,5466 | 0,0841 |

HLLL | 25,790 | 0,0128 | 15,1579 | 0,0716 |

HLLH | 19,170 | 0,0200 | 13,8335 | 0,0851 |

HLHL | 25,790 | 0,0220 | 16,1579 | 0,0798 |

HLHH | 35,420 | 0,0350 | 17,0842 | 0,0851 |

HHLL | 41,240 | 0,0250 | 18,2471 | 0, 0574 |

HHLH | 53,520 | 0,0450 | 20,7045 | 0,0589 |

HHHL | 55,930 | 0,0250 | 21,1859 | 0,0591 |

HHHH | 57,080 | 0,0125 | 21,4153 | 0,0628 |

The trajectories of the reference and actual locomotive speed in the simulation for the case (HHLL), which the lowest controller error was obtained, represented in Figure 4a. The entire trajectory of the error signal between the reference and actual locomotive speed conducted in the simulation depicted in Figure 4b.

As can be seen from Figure 4, the controller performance is satisfactory.

### 5.1 Time-domain analysis

After completing the simulation of all cases, data which are related to the force, speed, acceleration information of each element recorded. The system identification process [31] was made to evaluate the transfer function between the acceleration of the first wagon (output) and the locomotive force (input) for examining the step response of each case via “System Identification Toolbox”.

Identification conditions is set as discrete-time feed through, estimated number of poles and zeros are selected 4 and 3, respectively, the estimation method is the nonlinear least squares with a line search method, fit frequency range is [0-0.01571 rad/s], initial conditions are selected automatically, the initialization method is IV, termination tolerance is 1e-7, maximum number of iterations is 1000. Half partitioned data set, which belong to the acceleration of the second wagon and locomotive traction force as train data and validation data, respectively. The system identification results presented in Table 3.

Stiffness & damping parameter space | System identification results (Ts=1e-4) | ||||||
---|---|---|---|---|---|---|---|

Permutation | FPE | MSE | Fit to estimation (%) | Function evaluation # | Cost | Norm of step | First-order optimality |

LLLL | 9.95e-7 | 8.15e-16 | 99.5% | 172 | 8.15e-16 | 3.63e-12 | 0,4580 |

LLLH | 6.43e-17 | 2.94e-21 | 99.2% | 532 | 8.05e-22 | 2.61e-20 | 0,0691 |

LLHL | 3.87e-7 | 3.09e-16 | 99.3% | 260 | 3.09e-16 | 1.24e-17 | 0,8690 |

LLHH | 9.20e-16 | 7.23e-21 | 99.6% | 463 | 2.01e-21 | 6.54e-18 | 0,4670 |

LHLL | 1.23e-15 | 5.62e-21 | 99.7% | 677 | 1.60e-21 | 1.03e-21 | 0,0238 |

LHLH | 1.37e-18 | 3.85e-21 | 99.5% | 961 | 8.66e-22 | 1.51e-21 | 0.0525 |

LHHL | 1.24e-17 | 2.52e-21 | 99.9% | 659 | 6.19e-22 | 1.67e-21 | 0,0750 |

LHHH | 4.91e-16 | 6.98e-21 | 99.2% | 451 | 1.96e-21 | 2.09e-21 | 0,1690 |

HLLL | 3.57e-18 | 5.41e-21 | 99.3% | 567 | 1.46e-21 | 1.05e-21 | 0,0516 |

HLLH | 8.12e-7 | 2.87e-16 | 99.5% | 204 | 2.87e-16 | 6.28e-16 | 0,5700 |

HLHL | 6.63e-7 | 3.28e-16 | 99.3% | 182 | 3.28e-16 | 7.22e-12 | 0,7100 |

HLHH | 7.26e-17 | 4.38e-21 | 99.6% | 394 | 1.04e-21 | 6.22e-14 | 0,0437 |

HHLL | 3.87e-20 | 3.88e-21 | 99.2% | 1028 | 8.53e-22 | 2.50e-21 | 0,0199 |

HHLH | 3.93e-7 | 3.97e-16 | 99.4% | 299 | 3.96e-16 | 5.13e-17 | 0,7170 |

HHHL | 5.12e-18 | 5.95e-21 | 99.7% | 755 | 1.60e-21 | 4.11e-21 | 0,0962 |

HHHH | 2.17e-15 | 4.93e-20 | 99.3% | 556 | 9.28e-21 | 1.32e-21 | 0,1970 |

The list of the estimated transfer function is proposed in APPENDIX-II. After obtaining the estimated transfer function, the unit step response time characteristics (rise time, settling time, overshoot, peak time, peak, settling minimum, settling maximum) and steady-state error in unit ramp response, which correlated with the integral of error in step response, were evaluated and the results demonstrated in Table 4.

Stiffness & damping parameter space | Time-domain characteristics (input: locomotive force, output: acceleration of first wagon) (Ts=1e-4) (unit step response) | Ramp response | ||||||
---|---|---|---|---|---|---|---|---|

Permutation | Rise time | Settling time | Overshoot (%) | Peak time | Peak | Settling min. | Settling max. | Integral of error (ramp steady-state error) |

LLLL | 0,0984 | 16,0424 | 10,0458 | 0,8856 | 1,10 | 0,8924 | 1,1005 | 0,0241 |

LLLH | 0,0630 | 9,6578 | 8.65e-5 | 0,8985 | 1,00 | 1,0000 | 1,0000 | 0,0000 |

LLHL | 0,0908 | 14,1356 | 6,9475 | 0,8993 | 1,07 | 0,9337 | 1,0695 | 0,0233 |

LLHH | 0,0717 | 16,0063 | 2.89e-4 | 0,8996 | 1,00 | 1,0000 | 1,0000 | 0,0895 |

LHLL | 0,0622 | 12,9081 | 3.45e-4 | 0,8866 | 1,00 | 1,0000 | 1,0000 | 0,0000 |

LHLH | 0,0629 | 11,9214 | 1.24e-5 | 0,8962 | 1,00 | 1,0000 | 1,0000 | 0,0889 |

LHHL | 0, 0562 | 9, 3024 | 3.85e-5 | 0, 8959 | 1, 00 | 1, 0000 | 1, 0000 | 0, 0000 |

LHHH | 0,0724 | 16,0856 | 2.08e-4 | 0,9029 | 1,00 | 1,0000 | 1,0000 | 0,0895 |

HLLL | 0,0646 | 13,3827 | 1.94e-5 | 0,8873 | 1,00 | 1,0000 | 1,0000 | 0,0890 |

HLLH | 0,0990 | 18,1009 | 9,4985 | 0,8846 | 1,10 | 0,0990 | 1,0950 | 0,0240 |

HLHL | 0,0880 | 13,4079 | 8,7312 | 0,9025 | 1,09 | 0,9094 | 1,0873 | 0,0261 |

HLHH | 0,0648 | 12,0161 | 8.73e-5 | 0,9011 | 1,00 | 1,0000 | 1,0000 | 0,0000 |

HHLL | 0,0571 | 9,3533 | 2.10e-6 | 0,8840 | 1,00 | 1,0000 | 1,0000 | 0,1595 |

HHLH | 0,0838 | 11,6350 | 6,8591 | 0,8860 | 1,07 | 0,9303 | 1,0686 | 0,0241 |

HHHL | 0,0714 | 16,2612 | 2.27e-5 | 0,8977 | 1,00 | 1,0000 | 1,0000 | 0,0893 |

HHHH | 0,0856 | 54,2275 | 3.87e-4 | 1,9809 | 1,00 | 1,0000 | 1,0000 | 0,0899 |

If the spring stiffness increases, results in higher overshoot with more rigid values are obtained. The case (LHHL) performs better due to the fast response (lowest rise time, settling time) and minimum overshoot for safety and reliable operation.

After running the (LHHL) case in the simulation, the position difference (*Δx _{i}* =

*x*

_{1}−

*x*

_{i+1},

*for i*= 1, …, 9) moreover, the velocity difference (

*Δx*

^{.i}=

*x*

^{.1}−

*x*

^{.i+1},

*for*

*i*= 1, …, 9) of wagons concerning the locomotive were visualised for acceleration (0-10 sec.), curving (75-85 sec.), slope climbing (185-195 sec.) and deceleration (120-130 sec.) phase in Figure 5a, 5b, respectively. These phases were chosen by considering the critical changes in the force and acceleration responses at the simulation.

The locomotive input (traction) force is shown in Figure 6a. The forces applied on wagon couplers are calculated during simulation, and the results are demonstrated for acceleration, curving, slope climbing, deceleration phase in Figure 6b, respectively.

The Euclidean norm (2-norm) and maximum absolute sum (infinite-norm) of the coupler forces, wagon accelerations, locomotive force and controller error were evaluated for all phases and tabulated in Table 5, respectively.

Coupler number | Euclidean norm of coupler force (sum(abs(X). ^2)^ (1/2)) | |||
---|---|---|---|---|

Acceleration phase | Curving phase | Slope climbing phase | Deceleration phase | |

1 | 1,0881 e7 | 6,5859 e6 | 3,3520 e6 | 1,7959 e7 |

2 | 8,1790 e6 | 6,5733 e6 | 3,3167 e6 | 1,8004 e7 |

3 | 7,6021 e6 | 6,5628 e6 | 3,2948 e6 | 1,8046 e7 |

4 | 7,5245 e6 | 6,5538 e6 | 3,2699 e6 | 1,8085 e7 |

5 | 7,2580 e6 | 6,5479 e6 | 3,2549 e6 | 1,8123 e7 |

6 | 7,1546 e6 | 6,5472 e6 | 3,2426 e6 | 1,8158 e7 |

7 | 8,4859 e6 | 6,5412 e6 | 3,2293 e6 | 1,8191 e7 |

8 | 1,0050 e7 | 6,5153 e6 | 3,2854 e6 | 1,8222 e7 |

9 | 8,5480 e6 | 6,7601 e6 | 3,6738 e6 | 1,8252 e7 |

Coupler Number | Infinite Norm of Coupler Force (max(abs(X))) | |||

Acceleration Phase | Curving Phase | Slope Climbing Phase | Deceleration Phase | |

1 | 2,0407 e5 | 8,1541 e4 | 6,4030 e4 | 6,6105 e4 |

2 | 1,1723 e5 | 7,4836 e4 | 5,5003 e4 | 6,6070 e4 |

3 | 1,0271 e5 | 7,0911 e4 | 4,9925 e4 | 6,6724 e4 |

4 | 1,0200 e5 | 6,8055 e4 | 4,6354 e4 | 6,6999 e4 |

5 | 9,5294 e4 | 6,5795 e4 | 4,3614 e4 | 6,7428 e4 |

6 | 9,5384 e4 | 6,3923 e4 | 4,1412 e4 | 6,8059 e4 |

7 | 2,0198 e5 | 6,2410 e4 | 4,5772 e4 | 6,8313 e4 |

8 | 2,2348 e5 | 6,3988 e4 | 4,4586 e4 | 7,0401 e4 |

9 | 1,4802 e5 | 8,7190 e4 | 6,6110 e4 | 7,1669 e4 |

Wagon number | Euclidean norm of wagon acceleration (sum(abs(X).^2)^(1/2)) | |||

Acceleration phase | Curving phase | Slope climbing phase | Deceleration phase | |

1 | 103, 6771 | 73,1766 | 37,2448 | 199,5401 |

2 | 87,7673 | 73,0369 | 36,8519 | 200,0473 |

3 | 83,2560 | 72,9204 | 36,6084 | 200,5109 |

4 | 82,5249 | 72,8202 | 36,3326 | 200,9500 |

5 | 80,1823 | 72,7540 | 36,1655 | 201,3623 |

6 | 79,0968 | 72,7463 | 36,0284 | 201,7515 |

7 | 80,4709 | 72,6796 | 35,8807 | 202,1201 |

8 | 78,9832 | 72,3927 | 36,5045 | 202,4688 |

9 | 80,0325 | 75, 1122 | 40, 8194 | 202, 8003 |

Wagon number | Infinite norm of wagon acceleration ((max(abs(X)))) | |||

Acceleration phase | Curving phase | Slope climbing phase | Deceleration phase | |

1 | 0,9759 | 0,9060 | 0,7114 | 0,7345 |

2 | 0,9511 | 0,8315 | 0,6111 | 0,7341 |

3 | 0,9448 | 0,7879 | 0,5547 | 0,7414 |

4 | 0,9445 | 0,7562 | 0,5150 | 0,7444 |

5 | 0,9416 | 0,7311 | 0,4846 | 0,7492 |

6 | 0,9414 | 0,7103 | 0,4601 | 0,7562 |

7 | 0,9748 | 0,6934 | 0,5086 | 0,7590 |

8 | 0, 9760 | 0,7110 | 0,4954 | 0,7822 |

9 | 0,9646 | 0, 9688 | 0, 7346 | 0, 7963 |

Euclidean norm of locomotive force (sum(abs(X).^2)^(1/2)) | ||||

Acceleration phase | Curving phase | Slope climbing phase | Deceleration phase | |

1,2907 e8 | 6,2586 e7 | 4,3145 e7 | 1,7524 e8 | |

Infinite norm of locomotive force ((max(abs(X)))) | ||||

Acceleration phase | Curving phase | Slope climbing phase | Deceleration phase | |

7,3159 e5 | 3,5971 e5 | 1,6093 e5 | 6,1303 e5 | |

Euclidean norm of controller error (sum(abs(X).^2)^(1/2)) | ||||

Acceleration phase | Curving phase | Slope climbing phase | Deceleration phase | |

88, 1869 | 30,1110 | 29,4740 | 43,0090 | |

Infinite norm of controller error ((max(abs(X)))) | ||||

Acceleration phase | Curving phase | Slope climbing phase | Deceleration phase | |

1, 1752 | 0,7748 | 0,9014 | 0,6665 |

The first coupler element is subjected to more force from the viewpoint of the Euclidean norm approach in the acceleration phase, while the last coupler element is exposed to more strength in the curving, slope climbing and deceleration phases. The maximum coupling force applied to the eight-coupler component in the acceleration phase, and the maximum force applied to the final coupler element in all other steps. The same correlation pattern observed for the wagon acceleration trajectories. Euclidean norm of the locomotive force, which applied in deceleration phase is much more than involved in all other stages, and infinite norm of the locomotive force, which applied in acceleration phase is much more than used in all different stages.

From the control point of view, it can be inferred that the tracking speed reference in the acceleration phase is more laborious than all other phases. Generally, it can be concluded that the nearer couplers to the locomotive are subjected to the higher coupler forces at acceleration phase and farthest coupler to the train is subjected to the higher coupler forces at deceleration phase, curving phase and slope climbing phase. If the speed of the system increases, the most pressing of coupler forces also increases.

In conclusion, the situation of the 1., 8., and 9. coupler is quite critical than the other couplers; therefore it should be focused on these couplers through the help of the frequency domain analysis. The traction curve for locomotive and wagon system drawn for all phases in Figure 7.

The traction force remains under the curve of the adhesion force in all phases, which means that the excessive wheel slip or skid is prevented during train operation.

### 5.2 Frequency-domain analysis

Frequency domain analysis is sometimes implemented as an alternative method to solve and to investigate engineering problems. It is generally quicker than time-domain analysis as the time dependency of the terms are converted into harmonic counterparts by assuming the dynamics response of the system is the superposition of the eigenvectors (mode shapes) of the system with distinct oscillation periods (natural frequencies). Also, while analysing the dynamic systems, investigating the problems through frequency domain makes it more practical to relate the system parameters, namely stiffness, inertia (mass) and damping with the system response. However, frequency analysis is a linear method meaning the system properties are assumed to be constant through time, and the contact conditions, as well as the boundary conditions, needs to be settled and unchanging. Consequently, nonlinear systems are not categorically suitable for implementing frequency analysis directly. One can either include iterative schemes to calculate the response (in rare problems) as accurately as time domain analysis or simply use frequency analysis to investigate the problem for fixed (e.g. lowest and highest values) system properties to understand the problem boundaries.

Dynamic response of the system was further investigated in the frequency domain to understand the effect of resonances clearly. As the crucial stage of spectrum analysis is first to check out the natural frequencies of the significant mode shapes, modal analysis of the system performed in “ANSYS Mechanical” which is an engineering software based on finite element methodology. The primary question herein is how to handle the nonlinear stiffness behaviour of the couplers since the modal analysis in the frequency domain is a linear theory.

In this study, the stiffness behaviour was considered for two situations as given in Eq. (1). For the first modal case (*Δx*^{.i} < 0), the stiffness parameter of the couplers is equal to *k*_{i1} to check the lower frequency range, whereas it is taken as *k*_{i2} + *k*_{i3} to consider the upper limits for the second case (*Δx*^{.i} > 0).

The natural frequency results are shown below in Figure 8.

For low stiffness behaviour (first case), the natural frequency range (related with couplers) starts at 0.23 Hz and ends at 1.45 Hz, while the corresponding values are at 0.65 Hz and 4.06 Hz, respectively for high stiffness behaviour. In other words, the resonance modes in a longitudinal direction can be excited in the broader frequency range between 0.23 Hz and 4.06 Hz.

In the next step, the coupler forces in the time domain corresponding to the 1., 8. and 9. couplers were transformed into spectral components in a frequency domain. The sampling frequency is 10 kHz, and the number of data is equal to 100, 000, which means the spectral resolution is 0.1 Hz. The corresponding results are shown in Figure 9. As the dynamic response of this system exists in the lower frequency range, the upper limit of abscissa is fixed to 10 Hz.

According to the results, the spectral content of the coupler forces is in the frequency range between 0 and 4 Hz. This outcome is expected since the previously made modal analysis has shown that the significant resonant modes are up to 4 Hz.

The results may further be investigated either per phase or per coupler number. In the acceleration phase, there are two significant peaks (valid for 1. and 8. couplers) which appears at around 0.2 Hz and 3.3 Hz. This means that the dynamic response of the coupler force is mostly dominated by the modes at (or near) these frequency values. Moreover, according to the frequency values in Figure 8, the peaks are supposed to be related to low stiffness *k*_{i1} moreover, high stiffness *k*_{i2} + *k*_{i3} cases, respectively. In other words, the dynamic response of the acceleration phase includes both low and high stiffness behaviour. Similar deductions for the other phases can also be made. Table 6 summarises these significant inferences.

Phase type | Coupler 1 | Coupler 8 | Coupler 9 |
---|---|---|---|

Acceleration | Two distinct peaks, botd low and high stiffness modes | Two distinct peaks, botd low and high stiffness modes | One single peak, low stiffness mode |

Curving | Several distinct peaks, botd low and high stiffness modes | Two distinct peaks, botd low and high stiffness modes | One single peak, low stiffness mode |

Deceleration | One single peak, low stiffness mode | One single peak, low stiffness mode | One single peak, low stiffness mode |

Slope climbing | Several distinct peaks, botd low and high stiffness modes | Two distinct peaks, botd low and high stiffness modes | One single peak, low stiffness mode |

Phase type | Coupler 1 | Coupler 8 | Coupler 9 |
---|---|---|---|

Acceleration | Max: 9691 N | Max: 9481 N | Max: 8337 N |

RMS: 4066 N | RMS: 3519 N | RMS: 2792 N | |

Curving | Max: 2190 N | Max: 2544 N | Max: 2700 N |

RMS: 1010 N | RMS: 965 N | RMS: 1261 N | |

Deceleration | Max: 6490 N | Max: 8617 N | Max: 8617 N |

RMS: 1408 N | RMS: 1957 N | RMS: 2000 N | |

Slope climbing | Max: 2090 N | Max: 2785 N | Max: 2860 N |

RMS: 913 N | RMS: 898 N | RMS: 1162 N |

## 6 Limitations and discussion

This paper suggests a computer-oriented mathematical model obtain the transient response and coupling force of a long train subject to the various train operating conditions. The experimental approaches based on the field data is not included. However, the design steps were handled carefully by comparing the models build in the experimentally validated numerical results [32, 33, 34]

Moreover, the results obtained in [35] exhibit high consistency with the simplified model built in this paper from the viewpoint of the coupler force patterns. This train model developed using Simpack consists of a group of 50 freight cars with two locomotives, including all suspension stages. The model is capable of investigating the influence of the lateral vehicle dynamics on the coupler forces. The vehicle has three rigid bodies, two representing bogies and one the car body. The two bogies were constrained to the track centre line using a rigid joint permitting the motion along the track longitudinal direction. The car body was constrained using spherical joints granting all the rotations, which were modelled using a high stiffness bushing element. The force element has three translational and three rotational. The roll motion of the car body was also included applying two anti-roll stiffnesses modelled with two bushing elements among each bogie frame and car body. The anti-roll stiffnesses allow a roll angle less than 6^{∘} with a lateral acceleration of 1

The authors were also not investigated the longitudinal-lateral interaction dynamics of the railway vehicle originated from the coupler rotation in this paper; therefore, they limited the study within the constraints of the decoupled longitudinal dynamic mode. The structural characteristics and compressed stability of coupler affect the dynamic performance and running safety. For a long distributed power train, the longitudinal motions of the front and rear sides are not synchronised. However, couplers are necessitated to transmit the occurred longitudinal forces along the train centre line. Otherwise, significant lateral components of the longitudinal compressive forces, which cause rail overturning, gauge widening, track lateral motion, or derailment, are produced and are exerted to the wheelsets [36]. The coupler rotation angle decreases in a stepped pattern with an increased friction coefficient of the coupler-tail surfaces. If the train is modelled as in [36], the rotation angle is around 7^{∘} when *μ* ≤ 0.2, around 3^{∘} when 0.25 ≤ *μ* ≤ 0.45, and is less than 1^{∘} when *μ* = 0.5. The maximum values of the derailment coefficients and wheelset lateral forces vary similarly with the changes that occurred for the rotation angle. When *μ* ≤ 0.2, large coupler deflection and significant lateral components of the longitudinal coupler compressive forces are generated and transferred to the wheelsets via the locomotive secondary and primary suspensions leading to large derailment coefficients and wheelset lateral forces. When *μ* ≥ 0.25, the lateral components of the coupler reduce dramatically as the coupler rotation angle keep below 3^{∘}, and the derailment coefficients and wheelset lateral forces stay at a low level. Coupler destabilisation takes place as the friction coefficient increases. Long trains may induce instability, derailment, structural failures to the coupling system. The couplers among two wagons once connected allow a certain rotation along the vertical axis providing the train runs on a curve. Couplers exhibit high stiffness and small inertia as compared to the train inertia, therefore, including the effect of the coupler inertia in the dynamic model cause high-frequency oscillations which lowers the computational efficiency. Longitudinal train forces during braking, traction, and curving rely on the coupler design (degrees of freedom (DOF), impact force absorption, energy dissipation, i.e.). A tangent track was utilised in [37] to investigate the effect of the geometric nonlinearities resulting from the coupler rotational kinematic DOF. The investigation was performed to reveal the effect of angle among the coupler that can affect the lateral forces on the track and the vehicle derailment safety.

The performances of couplers with small and large rotation angle on a model consisting of four locomotives and single mass were compared in [38]. The results implied that couplers with low rotation angle demonstrate better performance in terms of derailment safety when running on curves. The behaviour of a train, which is composed of five detailed locomotives and wagons and additional simplified vehicles, on S-shaped curve with a narrow curve radius was also represented. It was deduced that small coupler angle less than 2^{∘} induces higher lateral forces reducing the safety against derailment.

The wheel-rail contact condition (e.g., dry, wet, and slippery rail conditions) affecting the traction dynamics was addressed in detail with numerical results characterising the general and structural characteristics of electromechanical systems in our previous papers [39, 24].

## 7 Conclusive summary and comments

As a result, the equivalence of the proposed simplified model with a more realistic and sophisticated model was established. A model simplification method explanation from an advanced friction draft gear model and the process by which simplification is achieved was introduced. The non-linear train system is simulated with a full non-linear hysteresis, and the trajectory tracking control was performed.

The effects of the system parameters, namely the stiffness and damping of the couplers revealed regarding the time-frequency responses after conducting full simulation for all possible cases. By comparing the value of the controller error, which is the difference between the reference and actual locomotive speed, the best-combined value of the system parameters proposed. On the other hand, the transfer function between the acceleration of the first wagon and the locomotive force estimated via system identification procedure and the unit step response time characteristics for each parameterised case was determined. It can be concluded that one specific case for system parameters (LHHL) demonstrates better performance concerning fast response and minimum overshoot for safety and reliable operation. By considering two different norm criteria; namely the most pressing and cumulative force exposure, the coupler forces, wagon accelerations, locomotive force and controller error have been evaluated under the condition of four different phases (acceleration, curving, slope climbing and deceleration). The results were utilised to select the first three critical couplers.

In the second section, the natural frequencies of the system were calculated by assuming the dynamic behaviour can be separated into low and high stiffness regions. By this means the significant longitudinal modes which are responsible for overshoots in the dynamic response of the system are determined together with the corresponding frequency range. Subsequently, spectral contents of the critical coupler forces in time domain were analysed for specified track sections by checking peak amplitudes and corresponding frequency values which are closely related with the low and high stiffness characteristics of the nonlinear couplers. Consequently, it can be deduced that modal and spectral analysis provides vital information about the dynamic response of such systems and its relationship with the coupler parameters.

This paper has new consideration for extended train dynamics analysis by using not only the time domain simulation but also from the viewpoint of the frequency domain, which is helpful for better understanding the in-train vibration behaviour.

## References

[1] Cole C., Longitudinal Train Dynamics, Handbook of Railway Vehicle Dynamics, 2006, 239-277, Taylor & Francis Group, LLC.10.1201/9781420004892.ch9Search in Google Scholar

[2] Belforte P., et al., Numerical and Experimental Approach for the Evaluation of Severe Longitudinal Dynamics of Heavy Freight Trains, Vehic. Syst. Dynam., 2008, 46(1), 937-955.10.1080/00423110802037180Search in Google Scholar

[3] Wu Q. et al., Longitudinal Train Dynamics: an Overview, Vehic. Syst. Dynam., 2016, 54(12), 1688-1714.10.1080/00423114.2016.1228988Search in Google Scholar

[4] Hsu, T.-K., and D. A. Peters., A Simple Dynamic Model For Simulating Draft-Gear Behavior In Rail-Car Impacts, J. Eng. Industry, 1978, 100(4), 492.10.1115/1.3439467Search in Google Scholar

[5] Wu Q., et al., Advanced Dynamic Modelling for Friction Draft Gears, Vehic. Syst. Dynam., 2015, 53(4), 475-492.10.1080/00423114.2014.1002504Search in Google Scholar

[6] Wu Q., et al., A Review of Dynamics Modelling of Friction Draft Gear, Vehic. Syst. Dynam., 2014, 52(6), 733-758.10.1080/00423114.2014.894199Search in Google Scholar

[7] Wu Q., Spiryagin M., Cole C., A Dynamic Model of Friction Draft Gear, Volume 6: 10th Int. Conf. Multibody Systems, Nonlin. Dyn. Contr., 2015, DETC2014-34159, V006T10A020.10.1115/DETC2014-34159Search in Google Scholar

[8] Ma W., et al., Coupler Dynamic Performance Analysis of Heavy Haul Locomotives, Vehic. Syst. Dynam., 2012, 50(9), 1435-1452.10.1080/00423114.2012.667134Search in Google Scholar

[9] Wu Q., Cole C., Computing Schemes for Longitudinal Train Dynamics: Sequential, Parallel and Hybrid, J. Comput. Nonlin. Dyn., 2015, 10(6), 064502.10.1115/1.4029716Search in Google Scholar

[10] Wu Q., et al., Parallel Computing Scheme for Three-Dimensional Long Train System Dynamics, J. Comput. Nonlin. Dyn., 2017, 12(4), 044502.10.1115/1.4035484Search in Google Scholar

[11] Andersen D.R., et al., Train Energy and Dynamics Simulator (TEDS): A State-of-the-Art Longitudinal Train Dynamics Simulator, ASME 2012 Rail Transport. Div. Fall Techn. Conf., 2012.10.1115/RTDF2012-9418Search in Google Scholar

[12] Genin J., et al., Longitudinal Track-Train Dynamics: A New Approach, J. Dyn. Syst., Meas. Contr., 1974, 96(4), 466-469.10.1115/1.3426847Search in Google Scholar

[13] Ling L., et al., A 3D Model for Coupling Dynamics Analysis of High-Speed Train/Track System, J. Zhejiang Univ. Sci. A, 2014, 15(12), 964-983.10.1007/978-981-10-5610-9_18Search in Google Scholar

[14] Davydov Y., Keyno M., Longitudinal Dynamics in Connected Trains, Proc. Eng., 2016, 165, 1490-1495.10.1016/j.proeng.2016.11.884Search in Google Scholar

[15] Saadat M., et al., Longitudinal Dynamics and Energy Flow Modelling for Diesel-Electric Locomotives, Int. J. Heavy Vehic. Syst., 2016, 23(2), 155.10.1504/IJHVS.2016.075504Search in Google Scholar

[16] Shabana A.A., et al., Use of the Non-Inertial Coordinates in the Analysis of Train Longitudinal Forces, J. Comput. Nonlin. Dyn., 2011, 7(1), 011001.10.1115/1.4004122Search in Google Scholar

[17] Duncan I.B., Webb P.A., The Longitudinal Behaviour of Heavy Haul Trains Using Remote Locomotives, 4th Int. Heavy Haul Conf., Brisbane, 1989, 587-590.Search in Google Scholar

[18] Wu Q., et al., Advanced Dynamic Modelling for Friction Draft Gears, Vehic. Syst. Dynam., 2015, 53(4), 475-492.10.1080/00423114.2014.1002504Search in Google Scholar

[19] Harak S.S., et al., Dynamic Analysis of Draft Gear and Draft Pad of Freight Wagon Due to Localized Defects Using FEM, Int. J. Acoust. Vibrat., 2016, 21, 3, 3422.10.20855/ijav.2016.21.3422Search in Google Scholar

[20] Harak S.S., et al., Effect of Multiple Location Defects on the Dynamics of Draft Gear Used in Freight Railway Wagon, Int. J. Vehic. Struct. Syst., 2015, 7(3), 107-113.10.4273/ijvss.7.3.04Search in Google Scholar

[21] Ansari M., et al., Longitudinal Dynamics of Freight Trains., Int. J. Heavy Vehic. Syst., 2009, 16(1/2), 102.10.1504/IJHVS.2009.023857Search in Google Scholar

[22] Wu Q., et al., Longitudinal Dynamics and Energy Analysis for Heavy Haul Trains., J. Modern Transport., 2014, 22(3), 127-136.10.1007/s40534-014-0055-xSearch in Google Scholar

[23] Garg V.K., Dukkipati R.V., Dynamics of Railway Vehicle Systems, 1984, Toronto, Acad. Press.10.1016/B978-0-12-275950-5.50016-XSearch in Google Scholar

[24] Uyulan Ç., Metin G., Modeling, Simulation and Re-Adhesion Control of an Induction Motor-Based Railway Electric Traction System., Proc. Instit. Mech. Eng. Part I: J. Syst. Control Eng., 2017, 232(1), 3-11.10.1177/0959651817732487Search in Google Scholar

[25] Cole C., et al., Modelling, Simulation and Applications of Longitudinal Train Dynamics, Vehic. Syst. Dynam., 2017, 55(10), 1498-1571.10.1080/00423114.2017.1330484Search in Google Scholar

[26] Wei W., et al., An Air Brake Model for Longitudinal Train Dynamics Studies, Vehic. Syst. Dynam., 2016, 55(4), 517-533.10.1080/00423114.2016.1254261Search in Google Scholar

[27] Spiryagin M., et al., Design and Simulation of Rail Vehicles, 2017, CRC Press, Taylor & Francis Group.Search in Google Scholar

[28] Spiryagin M., Design and Simulation of Heavy Haul Locomotives and Trains, 2017, CRC Press, Taylor & Francis Group.10.1201/9781315369792Search in Google Scholar

[29] Rochard B.P., Schmid F., A Review of Methods to Measure and Calculate Train Resistances, Proc. Instit. Mech. Eng. Part F: J. Rail Rapid Transit, 2000, 214(4), 185-199.10.1243/0954409001531306Search in Google Scholar

[30] Iwnicki S., et al., Handbook of Railway Vehicle Dynamics, 2020, CRC Press, Taylor & Francis Group.10.1201/9780429469398Search in Google Scholar

[31] Ljung L., System Identification: Theory for the User, 2012, Prentice Hall PTR.Search in Google Scholar

[32] Giorgio D., et al., Numerical and Experimental Investigation of Heavy Freight Train Dynamics, ASME 2007 Int. Mech. Eng. Congr. Expo. Vol. 16: Transport. Syst., IMECE2007-42693, 439-448.Search in Google Scholar

[33] Liu P., et al., Establishment and Verification of Three-Dimensional Dynamic Model for Heavy-Haul Train-Track Coupled System, Vehic. Syst. Dynam., 2016, 54(11), 1511-1537.10.1080/00423114.2016.1213862Search in Google Scholar

[34] Qi Z., et al., Simulation of Longitudinal Dynamics of Long Freight Trains in Positioning Operations, Vehic. Syst. Dynam., 2012, 50(9), 1409-1433.10.1080/00423114.2012.661063Search in Google Scholar

[35] Bosso N., Zampieri N., Long Train Simulation Using a Multibody Code, Vehic. Syst. Dynam., 2016, 55(4), 552-570.10.1080/00423114.2016.1267373Search in Google Scholar

[36] Zhang Z,, et al., Compressed Stability Analysis of the Coupler and Buffer System of Heavy-Haul Locomotives, Vehic. Syst. Dynam., 2015, 53(6), 833-855.10.1080/00423114.2015.1023318Search in Google Scholar

[37] Massa A., et al., Numerical Study of the Noninertial Systems: Application to Train Coupler Systems, Nonlin. Dyn., 2011, 68(1-2), 215-233.10.1007/s11071-011-0220-2Search in Google Scholar

[38] Ma W., et al., Coupler Dynamic Performance Analysis of Heavy Haul Locomotives, Vehic. Syst. Dyn., 2012, 50(9), 1435-1452.10.1080/00423114.2012.667134Search in Google Scholar

[39] Uyulan C., et al., Modeling, Simulation and Slip Control of a Railway Vehicle Integrated with Traction Power Supply, Cogent Eng., 2017, 4(1), 1312680.10.1080/23311916.2017.1312680Search in Google Scholar

[40] Ansari M., et al., Effects of Coupler Specifications and Operational Conditions on the Longitudinal Freight Train Dynamics., Volume 3: 19th Int. Conf. Design Theory and Methodology; 1st Int. Conf. Micro- and Nanosystems; and 9th Int. Conf. Advanced Vehicle Tire Technologies, Parts A and B, 2007, DETC2007-35809, 1279-1286.10.1115/DETC2007-35809Search in Google Scholar

## Appendix-I

Parameters and their numerical values used in the LTD simulation [40, 35, 15, 21].

Parameter | Numerical value |

b (_{i}low) | 999.8706 |

b (_{i}high) | 2512 |

k_{i1} (low) | 3.6086*10^{6} |

k_{i1} (high) | 8.265776*10^{6} |

k_{i2} (low) | 1.0079*10^{5} |

k_{i2} (high) | 1.856364*10^{6} |

k_{i3} (low) | 3.4748*10^{6} |

k_{i3} (high) | 7.998176*10^{6} |

m_{dyn,loco} | 74151 kg |

m_{wagon,i} | 90000 kg |

ϵ | 0.025 m |

A | 8.202*10^{3}N |

B | 0.10656*10^{3} |

C | 0.0119322*10^{3} |

κ | 200 |

R_{y} | 300 m |

h | 540 |

n | 9 |

g | 9.81 |

m_{stc} | 8.7663*10^{5}kg |

m_{dyn} | 8.8405*10^{5}kg |

ν | 0.6 |

## Appendix-II: List of the estimated transfer functions

**Received:**2019-07-11

**Accepted:**2019-11-12

**Published Online:**2020-02-07

© 2020 Caglar Uyulan and Ersen Arslan, published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.