Simulation and time-frequency analysis of the longitudinal train dynamics coupled with a nonlinear friction draft gear

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.


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 simpli cations and assumptions.
The main characteristics of the LTD classi ed 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 e ects 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 e cient 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 bu er models. The draft gears have a nonlinear and complex nature, which is generalised as look-up tables with various patterns. They can be classi ed as friction, rubber, hydraulic and types of their combinations.
Forces in draft gears are related to draft gear de ections. 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-de ection 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-o between e ciency and computational effort or consumed time.
The advanced models allow modelling a broader spectrum of dynamic behaviour and include signi cant 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 autocoupler 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 e ect 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 de ects 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 e ciency [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 veri ed through a com-mercial 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 dieselelectric 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 e ciency 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. Lowfrequency oscillations were classi ed into two di erent 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 uctuations correspond to the square wave oscillations originated from run-in/run-out behaviour. Draft gear's intrain characteristics a ect 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 di erent draft gear de ections and impact characteristics at di erent frequencies. The mode shapes and rates for an individual draft gear were found utilising the nite element approach in [19]. Which draft pad should be checked for fatigue response under the cyclic loading of 10 Hz, was clari ed within a range of frequency values. This work supplies a tool for detecting new frequencies of draft gear, whose modal frequencies are in uenced 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 identi cation approaches which can provide su ciently 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 fulllength 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 identi cation process, which derives the discrete transfer functions among the rst 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 di erent track/train phases (acceleration, curving, slope climbing, deceleration).
The frequency analysis starts with collecting natural frequencies in the longitudinal direction for two di erent coupler behaviour. Then, the previously obtained timedomain results transformed into the frequency domain, and the meaningful inferences about spectral e ects 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 e ciency 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, respec-tively. 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.

Mathematical modelling of the LTD . 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 de ned as a system of di erential 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 sti ness and damping coe cients of the wagon couplings (draft gear) are expressed regarding complex dynamical characteristic functions

. Friction draft gear model
Friction draft gears have friction and velocity-dependent sti ness nature. Nonlinear hysteresis, which leads to discontinuities among loading and unloading forcedisplacement curves, occurs due to the friction damping. Velocity-dependent friction, slack, limiting sti ness, 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 di erent stages, which express di erent 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 a ect 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 quasistatic 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).
corresponds to the four di erent stages (loading_1, load-ing_2, unloading_1, unloading_2) as explained above; F fdg,i is the external force applied on the central wedge; Fsm and Fsr are the mainspring force and release spring force, respectively; α, β and γ are wedge shoe angles as shown in Figure 2; µ , µ and µ are the coe cients of friction, as shown in Figure 2; µ refers to the friction coe cient of the moveable plates; km and kr are mainspring sti ness and release spring sti ness, respectively; x m and x r are the pre-de ection of the mainspring and the release spring, respectively; x fdg is the friction draft gear de ection; x rdg is the release spring de ection.
The exponential friction model integrated into the force-displacement characteristic presented in this study.
In the exponential friction model, the coe cient 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 e ect 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 is the relative velocity between the wedge shoe and stationary inner plate, v f is the initial impact velocity of the two adjacent vehicles; S(.) is the modi cation factor given by utilising a two-dimensional look-up table, which is the function of impact velocity and de ection; h , h and h are speci c model parameters, which reect kinetic friction, part of static friction; the transition from static to kinetic friction, respectively.
The friction model corresponding to the all nonlubricated coe cients 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 sti ness", which rst recorded in [22]. The forces leading longitudinal vibrations of the train discovered as a third type of intrain force. The locked sti ness is de ned as the longitudinal sti ness of the whole train when train oscillates longitudinally. The locked sti ness cannot be captured from the impact tests or drop-hammer tests. The locked sti ness can be stated mathematically as in Eq. (5)- (7).
where x t , x t−∆t , v t , v t−∆t are draft gear de ections and relative velocities of adjacent vehicles of the current and previous time-step, respectively. f l and fu are forcedisplacement characteristics for loading and unloading processes, respectively. k locked is the locked sti ness. 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 nal draft gear force.

. 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 i correspond to the acceleration m s , the velocity m s , is the displacement (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 i is the i. damping coe cient Ns m ; m dyn,loco is the dynamic mass of the locomotive at zero slip state (kg); F tr db ,loco is the traction or dynamic brake force of the locomotive (kN) ; F r,loco is the resultant resistance force of the locomotive (kN).

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 eld-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 coe cients 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 uid dynamics-based models [25]. The train air brake system model used in the simulation is adapted from [26].
The traction/dynamic brake force term F tr db ,loco must be repetitively updated according to the changes in locomotive speed and driver control demands. Due to the complexities occurred in controlling the complex locomotivewagon system with nonlinear friction draft gear, the model must be simpli ed, and parameters must be reduced to make it control by a proper PID controller.
The highly nonlinear coupler dynamics can be analytically transformed into a simpli ed but also nonlinear model given in Eq. (9) [27].
where P i is the preload, b i is the damping ratio ( Ns m ), k i is the linear sti ness N m , k i is the nonlinear sti ness N m , ν is the critical velocity constant, k i is the recoil sti ness N m , and ∆x i , ∆x .
i are relative displacement and velocity of the wagons, ϵ is the coupling clearance.
The terms k i and k i are the main cause of hysteresis behaviour. k i 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 simpli ed 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 de ned 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 Ns m and C Ns m are regression coe cients obtained by tting test data to the Davis equation. A and B are related to the mass, and mechanical resistance and Caccount 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 m s and κ is the gradient in the form scaled as " : κ" (for example a grade of 5% is calculated as κ = . = in the formula) Curving resistance obtained as in Eq. (14).
where h is the dimensionless parameter, which changes from 500 to 1200 depending upon the railway vehicle; Ryis the radius of the curved track (m).
The complex model has seven-dimensional parameter space such as [h , h , h , k locked , k m , kr , S] and four di erent stages. The simpli ed model has ve-dimensional parameter space such as [b, k , k , k , ν] moreover, two di erent stages.
To validate the reducibility of the complex dynamical model into a simpli ed analytical model, the parameters of the complex model and simpli ed model matched by conducting a simulation. "Simulink Design Optimization Toolbox" was utilised in this experiment. Single friction draft gear corresponding to the sophisticated and simpli ed 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; N, π rad/sec along 8 seconds, N, π rad/sec along 4 seconds, N, * π rad/sec along 2 seconds, N, * π rad/sec along 0.4 seconds, N, * π 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 simplied model w.r.t to the parameters given above. Then the parameters of the simpli ed 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 simpli ed model after the optimisation process are listed in Table 1.
The output of the simpli ed model and sophisticated model for low and high parameter set are represented in Figure 3a and 3b, respectively.
whereẋ is the locomotive speed km h ; m stc = n.m wagon,i + m stc,loco is the static mass of the loco-motive&wagon system (kg).
The coe cient of adhesion is determined by several parameters; namely the material properties, rail surface, wheel pro le, axle pressures and vertical movements. Ac-cording to UIC standards, the maximum value of the coefcient 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-Kni er" formulation, which varies with speed, is chosen to be utilised for determining the adhesion coe cient. The maximum traction power (Pτ, max in kW)is de ned by using a maximum traction force (Fτ,max = m dyn .amax + Fr), without considering the curving and gravitational resistances except for acceleration and braking, and the speed at which the vehicle reaches its maximum power v Pmax in m s [30]. Maximum traction power calculated as in Eq. (16).
where amax is the maximum available acceleration m s .

LTD Simulation
LTD simulations were conducted for a locomotive with nine wagons. "ode3" (Bogacki-Shampine) with a xedstep 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 (Ry = m) from t = s. to = s. At t = − s, the system is climbing % 5 slope with a full load. The operational speed at maximum motor power v Pmax was considered as .   Table 2.
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.

. 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 identi cation process [31] was made to evaluate the transfer function between the acceleration of the rst wagon (output) and the locomotive force (input) for examining the step response of each case via "System Identi cation Toolbox".
Identi cation 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, t 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 identication results presented in Table 3. 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.
If the spring sti ness 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.
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 (in nite-norm) of the coupler forces, wagon accelerations, locomotive force and controller error were evaluated for all phases and tabulated in Table 5, respectively.
The rst 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 nal 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 in nite norm of the locomotive force, which applied in acceleration phase is much more than used in all di erent 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 fre-quency 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.

. 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 sti ness, 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 xed (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 e ect of resonances clearly. As the crucial stage of spectrum analysis is rst to check out the natural frequencies of the signi cant mode shapes, modal analysis of the system performed in "ANSYS Mechanical" which is an engineering software based on nite element methodology. The primary question herein is how to handle the nonlinear stiness behaviour of the couplers since the modal analysis in the frequency domain is a linear theory.
In this study, the sti ness behaviour was considered for two situations as given in Eq. (1). For the rst modal case (∆x . i < ), the sti ness parameter of the couplers is equal to k i to check the lower frequency range, whereas it is taken as k i + k i to consider the upper limits for the second case (∆x . i > ). The natural frequency results are shown below in Figure 8.
For low sti ness behaviour ( rst 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 sti ness 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 xed 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 signi cant 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 signi cant 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 sti ness k i moreover, high sti ness k i + k i cases, respectively. In other words, the dynamic response of the acceleration phase includes both low and high sti ness behaviour. Sim-ilar deductions for the other phases can also be made. Table 6 summarises these signi cant inferences.
Finally, the couplers and phases can be compared in terms of maximum and average (RMS between 0-10 Hz) coupler force values. Corresponding results are given in Table 7.

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 eld 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 simpli ed 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 in uence 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 sti ness 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 sti nesses modelled with two bushing elements among each bogie frame and car body. The anti-roll sti nesses allow a roll angle less than • with a lateral acceleration of m s . The cars were coupled in groups of two employing a steel bar, which connects the draft gears without slack. Couplers include pointto-point nonlinear spring-damper element with hysteresis. Due to the assumptions stated in the simpli ed model, the adopted elements do not react with torque in case of relative rotations among couplers. The locomotive model is consists of 27 rigid bodies:6 wheelsets each one connected to two axle-boxes using revolute joints, six motors and gearboxes, two bogie frames and a car body, respectively. Each axle-box was connected to the bogie frame utilising two springs and a traction rod with a point-to-point spring damper. The secondary suspension is composed of four rubber springs connecting the bogie frame to the coach and by two trailing rods embedded anti-symmetrically to the bogie centre [35].
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 a ect 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, signi cant 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 coe cient of the coupler-tail surfaces. If the train is modelled as in [36], the rotation angle is around • when µ ≤ . , around • when . ≤ µ ≤ . , and is less than • when µ = . . The maximum values of the derailment coe cients and wheelset lateral forces vary similarly with the changes that occurred for the rotation angle. When µ ≤ . , large coupler de ection and signi cant 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 coe cients and wheelset lateral forces. When µ ≥ . , the lateral components of the coupler reduce dramatically as the coupler rotation angle keep below • , and the derailment coe cients and wheelset lateral forces stay at a low level. Coupler destabilisation takes place as the friction coe cient 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 sti ness and small inertia as compared to the train inertia, therefore, including the e ect of the coupler inertia in the dynamic model cause high-frequency oscillations which lowers the computational e ciency. 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 e ect of the geometric nonlinearities resulting from the coupler rotational kinematic DOF. The investigation was performed to reveal the e ect of angle among the coupler that can a ect 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 ve detailed locomotives and wagons and additional simpli ed vehicles, on S-shaped curve with a narrow curve radius was also represented. It was deduced that small coupler angle less than • induces higher lateral forces reducing the safety against derailment.
The wheel-rail contact condition (e.g., dry, wet, and slippery rail conditions) a ecting 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].

Conclusive summary and comments
As a result, the equivalence of the proposed simpli ed model with a more realistic and sophisticated model was established. A model simpli cation method explanation from an advanced friction draft gear model and the process by which simpli cation is achieved was introduced.
The non-linear train system is simulated with a full nonlinear hysteresis, and the trajectory tracking control was performed.
The e ects of the system parameters, namely the stiness 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 di erence 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 rst wagon and the locomotive force estimated via system identi cation procedure and the unit step response time characteristics for each parameterised case was determined. It can be concluded that one speci c case for system parameters (LHHL) demonstrates better performance concerning fast response and minimum overshoot for safety and reliable operation. By considering two di erent 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 di erent phases (acceleration, curving, slope climbing and deceleration). The results were utilised to select the rst 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 sti ness re-gions. By this means the signi cant 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 speci ed track sections by checking peak amplitudes and corresponding frequency values which are closely related with the low and high sti ness 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 intrain vibration behaviour.