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

• Caglar Uyulan and Ersen Arslan
From the journal Nonlinear Engineering

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

Figure 1

The schematic representation of the generalised train model

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

fcw,ixi,xi1,x.i,x.i1,fcw,nxn,xn1,x.n,x.n1.

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

Figure 2

Frction draft gear model [5]

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

Ffdg,i=ψi.Fsmψi1.Fsr(i=1,?,4)(1)

where

ψ1=1+tanβ+arctanμ3.tany+arctanμ11tanα+arctanμ2.tany+arctanμ1ψ2=ψ1+21μ1tany.μ4.ψ11μ1+tanyψ3=1+tanβarctanμ3.tanyarctanμ11tanαarctanμ2.tanyarctanμ1ψ4=tanyμ1.ψ3tany.12μ1μ4+2μ1μ4ψ3+2μ4ψ32μ4μ1Fsm=kmxm0+cosα.cosyβcosα+y.cosβ.xfdgFsr=krxr0+xrdg

where

(i=1,?,4)

corresponds to the four different stages (loading_1, loading_2, unloading_1, unloading_2) as explained above; Ffdg,i is the external force applied on the central wedge; Fsm and Fsr are the mainspring force and release spring force, respectively; α, β 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; km and kr are mainspring stiffness and release spring stiffness, respectively; xm0 and xr0 are the pre-deflection of the mainspring and the release spring, respectively; xfdg is the friction draft gear deflection; xrdg is the release spring deflection.

xrdg=xfdg.cosα.cosyβcosα+y.cosβ1

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

μ1v1,vf0,xfdg=h1+h2.exph3.v1S(vf0,xfdg)(2)

where v1 is the relative velocity between the wedge shoe and stationary inner plate, vf0 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; h1, h2 and h3 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).

μjvj=h1+h2.exph3vj(j=2,4)(3)

where vj 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).

v1=cos(α)cos(α+y).vfdgv2=sin(y)cos(α+y).vfdgv3=cosα.sin(y)cosα+y.cosβ.vfdgv4=cosα.cos(yβ)cosα+y.cos(β).vfdg,i=4vfdg,i=2(4)

where vfdg 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).

Ffdg,i(xt,vt)=flxt,vt,xtxtΔtfuxt,vt,xt<xtΔt(5)
Ffdg,ixt,vt=Ffdg,ixtΔt,vtΔt+klocked(xtxtΔt(6)
fuxt,vtFfdg,i(xt,vt)flxt,vt(7)

where xt, xtΔt, vt, vtΔt are draft gear deflections and relative velocities of adjacent vehicles of the current and previous time-step, respectively. fl and fu are force-displacement characteristics for loading and unloading processes, respectively. klocked 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.

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

mdyn,loco.x¨1+Ffdg,12,1=Ftrdb,locoFr,locoΔx.1<0,Δx1<ϵmdyn,loco.x¨1+Ffdg,34,1=Ftrdb,locoFr,locoΔx.1>0,Δx1>ϵmwagon,i.x¨i+Ffdg,12,i=Fr,iΔx.i<0,Δxi<ϵmwagon,i.x¨i+Ffdg,34,i=Fr,iΔx.i>0,Δxi>ϵmwagon,n.x¨nFfdg,12,n=Fr,nΔx.n<0,Δxn<ϵmwagon,n.x¨nFfdg,34,n=Fr,nΔx.n>0,Δxn>ϵ
(i=2,...,n1)(8)

where i, i, xi correspond to the acceleration ms2, the velocity ms, is the displacement (m) of each element, respectively; mwagon,i is the mass (kg); Fr,i is the resultant resistance force (kN) of the i. element; bi is the i. damping coefficient Nsm; mdyn,loco is the dynamic mass of the locomotive at zero slip state (kg); Ftrdb,loco is the traction or dynamic brake force of the locomotive (kN); Fr,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 Ftrdb,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 locomotive-wagon system with nonlinear friction draft gear, the model must be simplified, 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 simplified but also nonlinear model given in Eq. (9) [27].

fcw,i=Pi+ki1.Δxi+bi.(Δx.i)Δx.i<0,Δxi<ϵfcw,i=Pi+ki3.Δxi+ki2.Δxi.eΔx.iν+bi.(Δx.i)Δx.i>0,Δxi>ϵ(9)

where Pi is the preload, bi is the damping ratio (Nsm), ki3 is the linear stiffness Nm, ki2 is the nonlinear stiffness Nm, ν is the critical velocity constant, ki1 is the recoil stiffness Nm, and Δxi, Δx.i are relative displacement and velocity of the wagons, ϵ is the coupling clearance.

The terms ki1 and ki3 are the main cause of hysteresis behaviour. ki2 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).

mdyn,loco.x¨1+b1.(x.1x.2)+k11.(x1x2)=Ftrdb,locoFr,locoΔx.1<0,Δx1<ϵ
mdyn,loco.x¨1+b1.(x.1x.2)+k13.(x1x2)+k12.(x1x2).e((x˙1x˙2)ν)=Ftrdb,locoFr,locoΔx.1>0,Δx1>ϵ
mwagon,i.x¨i+bi.(x.ix.i+1)+bi1.(x.ix.i1)+ki1.(xixi+1)+k(i1)1.(xixi1)=Fr,iΔx.i<0,Δxi<ϵ
mwagon,i.x¨i+bi.(x.ix.i+1)+bi1.(x.ix.i1)+ki3.(xixi+1)+ki2.(xixi+1).e((x˙ix˙i+1)ν)
+k(i1)3.(xixi1)+k(i1)2.(xixi1).e((x˙ix˙i1)ν)=Fr,iΔx.i>0,Δxi>ϵ
mwagon,n.x¨n+bn1.(x.nx.n1)+k(n1)1.(xnxn1)=Fr,n4Δx.n<0,Δxn<ϵ
mwagon,n.x¨n+bn1.(x.nx.n1)+k(n1)3.(xnxn1)+k(n1)2.(xnxn1).e((x˙nx˙n1)ν)=Fr,n(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).

Fr,loco=Fp,loco+Fcrv,loco+Fgrv,loco+Fbrake,loco
Fr,i=Fp,i+Fcrv,i+Fgrv,i+Fbrake,i(11)

where Fp,loco and Fp,i are the propulsion resistances which are represented as the sum of the aerodynamic and rolling resistances (kN); Fcrv,loco and Fcrv,i are the resistance force elements originating from the track curvature (kN); Fgrv,loco and Fgrv,i are the gravitational force elements due to the track gradient (kN); Fbrake,loco and Fbrake,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].

Fp,i=A+Bx.i+Cx.i2(12)

where A(N), BNsm and CNs2m2 are regression coefficients obtained by fitting 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).

Fgrv,i=mdyngκ(13)

where mdyn is the dynamic mass of the vehicle system (kg); g is the acceleration due to gravity ms2 and κ is the gradient in the form scaled as ”1:κ” (for example a grade of 5% is calculated as κ = 10.05 = 20 in the formula) Curving resistance obtained as in Eq. (14).

Fcrv,i=0.1hRy(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 [h1, h2, h3, klocked, km, kr, S] and four different stages. The simplified model has five-dimensional parameter space such as [b, k1, k2, k3, ν] 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; 106 N, π2 rad/sec along 8 seconds, 106 N, π rad/sec along 4 seconds, 106 N, 2*π rad/sec along 2 seconds, 107 N, 10*π rad/sec along 0.4 seconds, 107N, 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.

Table 1

The parameter space corresponding to the complex and simplified models

Complex model
h1(lowL,0.05; highH, 0.0525)
h2(lowL, 0.01; highH, 0.0105)
h3(lowL, 0.4; highH, 0.42)
klockedlowL,4106Nm;highH,8106Nm
kmlowL,1106Nm;highH,9106Nm
krlowL,1.5106Nm;highH,7.5106Nm
S(lowL, 0.03; highH, 0.0315)
Simplified model
bi(lowL,999.8706Nsm;highH,2512Nsm)
ki1(lowL,3.6086106Nm;highH,8.265776106Nm)
ki2(lowL,1.0079105Nm;highH,1.856364106Nm)
ki3lowL,3.4748106Nm;highH,7.998176106Nm
ν(lowL, 0.1018; highH, 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.

Figure 3

Friction draft gear forces a) for the low parameter set; b) for the high parameter set

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

2. 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).

Fa,maxx˙1=mstc.g.(0.161+7.5x˙1+44)(15)

where 1 is the locomotive speed kmh; mstc = n.mwagon,i + mstc,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τ, maxinkW)is defined by using a maximum traction force (Fτ, max = mdyn.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 vPmaxinms [30].

Maximum traction power calculated as in Eq. (16).

Pτ,max=Fτ,maxvPmax(16)

where amax is the maximum available acceleration ms2.

## 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 (Ry = 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 (vPmax) was considered as 16.72 ms. Maximum allowable operational acceleration (amax) moreover, deceleration are 0.976 ms2 moreover, − 0.976 ms2, respectively. The coupler force is the most critical variable in the LTD.

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 + I.Tsz1+D.N1+N.Ts.zz1) to track the reference locomotive speed. The parameters of the PID controller was tuned by considering the reference tracking which its response time and transient behaviour become 1.02 second, 0.75, respectively. The integrator and filter methods are selected as forward and backward Euler, respectively. The filter coefficient (N) is selected as 100 and the sample time is Ts = 1e − 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.

Table 2

The parameters of the PID controller and RMS values of the controller error for all cases in the stiffness and damping parameter space

Stiffness and damping parameter spacePID Controller parametersControl results
PermutationProportionalDerivativeIntegralController Error (RMS)
LLLL12,2920,014624,58500,0634
LLLH34,3000,014716,86000,0870
LLHL33,0800,012516,61510,0819
LLHH29,2300,014615,84530,1002
LHLL31,2200,016216,24400,0624
LHLH47,7200,017519,54360,0633
LHHL31,2200,017516,24400,0709
LHHH22,7300,016014,54660,0841
HLLL25,7900,012815,15790,0716
HLLH19,1700,020013,83350,0851
HLHL25,7900,022016,15790,0798
HLHH35,4200,035017,08420,0851
HHLL41,2400,025018,24710, 0574
HHLH53,5200,045020,70450,0589
HHHL55,9300,025021,18590,0591
HHHH57,0800,012521,41530,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.

Figure 4

a) Trajectories of the reference and actual locomotive speed in the simulation; b) Error signal

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.

Table 3

System identification results

Stiffness & damping parameter spaceSystem identification results (Ts=1e-4)
PermutationFPEMSEFit to estimation (%)Function evaluation #CostNorm of stepFirst-order optimality
LLLL9.95e-78.15e-1699.5%1728.15e-163.63e-120,4580
LLLH6.43e-172.94e-2199.2%5328.05e-222.61e-200,0691
LLHL3.87e-73.09e-1699.3%2603.09e-161.24e-170,8690
LLHH9.20e-167.23e-2199.6%4632.01e-216.54e-180,4670
LHLL1.23e-155.62e-2199.7%6771.60e-211.03e-210,0238
LHLH1.37e-183.85e-2199.5%9618.66e-221.51e-210.0525
LHHL1.24e-172.52e-2199.9%6596.19e-221.67e-210,0750
LHHH4.91e-166.98e-2199.2%4511.96e-212.09e-210,1690
HLLL3.57e-185.41e-2199.3%5671.46e-211.05e-210,0516
HLLH8.12e-72.87e-1699.5%2042.87e-166.28e-160,5700
HLHL6.63e-73.28e-1699.3%1823.28e-167.22e-120,7100
HLHH7.26e-174.38e-2199.6%3941.04e-216.22e-140,0437
HHLL3.87e-203.88e-2199.2%10288.53e-222.50e-210,0199
HHLH3.93e-73.97e-1699.4%2993.96e-165.13e-170,7170
HHHL5.12e-185.95e-2199.7%7551.60e-214.11e-210,0962
HHHH2.17e-154.93e-2099.3%5569.28e-211.32e-210,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.

Table 4

Time-domain characteristics of the estimated transfer functions

Stiffness & damping parameter spaceTime-domain characteristics (input: locomotive force, output: acceleration of first wagon) (Ts=1e-4) (unit step response)Ramp response
PermutationRise timeSettling timeOvershoot (%)Peak timePeakSettling min.Settling max.Integral of error (ramp steady-state error)
LLLL0,098416,042410,04580,88561,100,89241,10050,0241
LLLH0,06309,65788.65e-50,89851,001,00001,00000,0000
LLHL0,090814,13566,94750,89931,070,93371,06950,0233
LLHH0,071716,00632.89e-40,89961,001,00001,00000,0895
LHLL0,062212,90813.45e-40,88661,001,00001,00000,0000
LHLH0,062911,92141.24e-50,89621,001,00001,00000,0889
LHHL0, 05629, 30243.85e-50, 89591, 001, 00001, 00000, 0000
LHHH0,072416,08562.08e-40,90291,001,00001,00000,0895
HLLL0,064613,38271.94e-50,88731,001,00001,00000,0890
HLLH0,099018,10099,49850,88461,100,09901,09500,0240
HLHL0,088013,40798,73120,90251,090,90941,08730,0261
HLHH0,064812,01618.73e-50,90111,001,00001,00000,0000
HHLL0,05719,35332.10e-60,88401,001,00001,00000,1595
HHLH0,083811,63506,85910,88601,070,93031,06860,0241
HHHL0,071416,26122.27e-50,89771,001,00001,00000,0893
HHHH0,085654,22753.87e-41,98091,001,00001,00000,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 (Δxi = x1xi+1, for i = 1, …, 9) moreover, the velocity difference (Δx.i = x.1x.i+1, fori = 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.

Figure 5

a) Position difference of wagons, b) Velocity difference of carriages relative to the locomotive

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.

Figure 6

a) The tractive force; b) pressing forces exerted on couplers changing with time

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.

Table 5

The norms of the coupler forces, wagon accelerations, locomotive force and controller error for each phase

Coupler numberEuclidean norm of coupler force (sum(abs(X). ^2)^ (1/2))
Acceleration phaseCurving phaseSlope climbing phaseDeceleration phase
11,0881 e76,5859 e63,3520 e61,7959 e7
28,1790 e66,5733 e63,3167 e61,8004 e7
37,6021 e66,5628 e63,2948 e61,8046 e7
47,5245 e66,5538 e63,2699 e61,8085 e7
57,2580 e66,5479 e63,2549 e61,8123 e7
67,1546 e66,5472 e63,2426 e61,8158 e7
78,4859 e66,5412 e63,2293 e61,8191 e7
81,0050 e76,5153 e63,2854 e61,8222 e7
98,5480 e66,7601 e63,6738 e61,8252 e7
Coupler NumberInfinite Norm of Coupler Force (max(abs(X)))
Acceleration PhaseCurving PhaseSlope Climbing PhaseDeceleration Phase
12,0407 e58,1541 e46,4030 e46,6105 e4
21,1723 e57,4836 e45,5003 e46,6070 e4
31,0271 e57,0911 e44,9925 e46,6724 e4
41,0200 e56,8055 e44,6354 e46,6999 e4
59,5294 e46,5795 e44,3614 e46,7428 e4
69,5384 e46,3923 e44,1412 e46,8059 e4
72,0198 e56,2410 e44,5772 e46,8313 e4
82,2348 e56,3988 e44,4586 e47,0401 e4
91,4802 e58,7190 e46,6110 e47,1669 e4
Wagon numberEuclidean norm of wagon acceleration (sum(abs(X).^2)^(1/2))
Acceleration phaseCurving phaseSlope climbing phaseDeceleration phase
1103, 677173,176637,2448199,5401
287,767373,036936,8519200,0473
383,256072,920436,6084200,5109
482,524972,820236,3326200,9500
580,182372,754036,1655201,3623
679,096872,746336,0284201,7515
780,470972,679635,8807202,1201
878,983272,392736,5045202,4688
980,032575, 112240, 8194202, 8003
Wagon numberInfinite norm of wagon acceleration ((max(abs(X))))
Acceleration phaseCurving phaseSlope climbing phaseDeceleration phase
10,97590,90600,71140,7345
20,95110,83150,61110,7341
30,94480,78790,55470,7414
40,94450,75620,51500,7444
50,94160,73110,48460,7492
60,94140,71030,46010,7562
70,97480,69340,50860,7590
80, 97600,71100,49540,7822
90,96460, 96880, 73460, 7963
Euclidean norm of locomotive force (sum(abs(X).^2)^(1/2))
Acceleration phaseCurving phaseSlope climbing phaseDeceleration phase
1,2907 e86,2586 e74,3145 e71,7524 e8
Infinite norm of locomotive force ((max(abs(X))))
Acceleration phaseCurving phaseSlope climbing phaseDeceleration phase
7,3159 e53,5971 e51,6093 e56,1303 e5
Euclidean norm of controller error (sum(abs(X).^2)^(1/2))
Acceleration phaseCurving phaseSlope climbing phaseDeceleration phase
88, 186930,111029,474043,0090
Infinite norm of controller error ((max(abs(X))))
Acceleration phaseCurving phaseSlope climbing phaseDeceleration phase
1, 17520,77480,90140,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.

Figure 7

Traction curve for locomotive and wagon system; a) Acceleration phase; b) Curving phase; c) Slope climbing phase; d) Deceleration phase

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 ki1 to check the lower frequency range, whereas it is taken as ki2 + ki3 to consider the upper limits for the second case (Δx.i > 0).

The natural frequency results are shown below in Figure 8.

Figure 8

Natural frequencies of the system in the longitudinal direction

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.

Figure 9

Coupler forces in the frequency domain; blue colour: coupler no:1, redcolour: coupler no:8 and orange colour: coupler no:9

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 ki1 moreover, high stiffness ki2 + ki3 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.

Table 6

The significant inferences about spectral results of three couplers for each phase

Phase typeCoupler 1Coupler 8Coupler 9
AccelerationTwo distinct peaks, botd low and high stiffness modesTwo distinct peaks, botd low and high stiffness modesOne single peak, low stiffness mode
CurvingSeveral distinct peaks, botd low and high stiffness modesTwo distinct peaks, botd low and high stiffness modesOne single peak, low stiffness mode
DecelerationOne single peak, low stiffness modeOne single peak, low stiffness modeOne single peak, low stiffness mode
Slope climbingSeveral distinct peaks, botd low and high stiffness modesTwo distinct peaks, botd low and high stiffness modesOne single peak, low stiffness mode

Table 7

The significant inferences about spectral results of three couplers for each phase

Phase typeCoupler 1Coupler 8Coupler 9
AccelerationMax: 9691 NMax: 9481 NMax: 8337 N
RMS: 4066 NRMS: 3519 NRMS: 2792 N
CurvingMax: 2190 NMax: 2544 NMax: 2700 N
RMS: 1010 NRMS: 965 NRMS: 1261 N
DecelerationMax: 6490 NMax: 8617 NMax: 8617 N
RMS: 1408 NRMS: 1957 NRMS: 2000 N
Slope climbingMax: 2090 NMax: 2785 NMax: 2860 N
RMS: 913 NRMS: 898 NRMS: 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 ms2. The cars were coupled in groups of two employing a steel bar, which connects the draft gears without slack. Couplers include point-to-point nonlinear spring-damper element with hysteresis. Due to the assumptions stated in the simplified 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 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 bi (low) 999.8706 Nsm bi (high) 2512 Nsm ki1 (low) 3.6086*106Nm ki1 (high) 8.265776*106Nm ki2 (low) 1.0079*105Nm ki2 (high) 1.856364*106Nm ki3 (low) 3.4748*106Nm ki3 (high) 7.998176*106Nm mdyn,loco 74151 kg mwagon,i 90000 kg ϵ 0.025 m A 8.202*103N B 0.10656*103Nsm C 0.0119322*103Ns2m2 κ 200 Ry 300 m h 540 n 9 g 9.81 ms2 mstc 8.7663*105kg mdyn 8.8405*105kg ν 0.6 ms

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

HLLLL(z)=0.972+0.5093z1+0.3467z2+0.5368z31+0.3974z1+0.5325z2+0.4519z3+0.02797z4HLLLH(z)=1+0.1871z1+0.2661z2+0.3232z31+0.1871z1+0.2661z2+0.3232z3+2.571107z4HLLHL(z)=0.9799+0.3041z1+0.3391z2+0.4271z31+0.2239z1+0.4595z2+0.3469z3+0.02006z4HLLHH(z)=1+0.4581z1+0.471z2+0.48z31+0.4581z1+0.471z2+0.48z3+7.962107z4HLHLL(z)=1+0.4798z1+0.4286z2+0.3648z31+0.4798z1+0.4286z2+0.3648z3+9.275107z4HLHLH(z)=1+0.2387z1+0.294z2+0.3344z31+0.2387z1+0.294z2+0.3344z3+3.647108z4HLHHL(z)=1+0.1734z1+0.2141z2+0.2498z31+0.1734z1+0.2141z2+0.2498z3+1.142107z4HLHHH(z)=1+0.4723z1+0.4739z2+0.4933z31+0.4723z1+0.4739z2+0.4933z3+5.715107z4HHLLL(z)=1+0.3813z1+0.4262z2+0.383z31+0.3813z1+0.4262z2+0.383z3+5.382108z4HHLLH(z)=0.9734+0.4587z1+0.3876z2+0.5393z31+0.3525z1+0.5469z2+0.4331z3+0.02655z4HHLHL(z)=0.9755+0.4287z1+0.3183z2+0.4468z31+0.3307z1+0.4562z2+0.3489z3+0.02449z4HHLHH(z)=1+0.3256z1+0.3438z2+0.3678z31+0.3256z1+0.3438z2+0.3678z3+2.487107z4HHHLL(z)=1+0.328z1+0.341z2+0.2834z31+0.328z1+0.341z2+0.2834z3+5.754109z4HHHLH(z)=0.9805+0.3583z1+0.2993z2+0.3961z31+0.2802z1+0.4164z2+0.3181z3+0.01951z4HHHHL(z)=1+0.3671z1+0.4512z2+0.4662z31+0.3671z1+0.4512z2+0.4662z3+6.404108z4HHHHH(z)=1+0.9142z1+0.8917z2+0.8144z31+0.9142z1+0.8917z2+0.8144z36.806107z4