Simulation and time-frequency analysis of the longitudinal train dynamics coupled with a nonlinear friction 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,xi−1,x.i,x.i−1,fcw,nxn,xn−1,x.n,x.n−1.

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]

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

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

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

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

where

ψ1=1+tanβ+arctanμ3.tany+arctanμ11−tanα+arctanμ2.tany+arctanμ1ψ2=ψ1+21−μ1tany.μ4.ψ1−1μ1+tanyψ3=1+tanβ−arctanμ3.tany−arctanμ11−tanα−arctanμ2.tany−arctanμ1ψ4=tany−μ1.ψ3tany.1−2μ1μ4+2μ1μ4ψ3+2μ4ψ3−2μ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; F_fdg,i is the external force applied on the central wedge; F_sm and F_sr are the mainspring force and release spring force, respectively; α, β and y are wedge shoe angles as shown in Figure 2; μ₁, μ₂ and μ₃ are the coefficients of friction, as shown in Figure 2; μ₄ refers to the friction coefficient of the moveable plates; k_m and k_r are mainspring stiffness and release spring stiffness, respectively; x_m0 and x_r0 are the pre-deflection of the mainspring and the release spring, respectively; x_fdg is the friction draft gear deflection; x_rdg 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.exp−h3.v1∗S(vf0,xfdg)(2)

where v₁ is the relative velocity between the wedge shoe and stationary inner plate, v_f0 is the initial impact velocity of the two adjacent vehicles; S(.) is the modification factor given by utilising a two-dimensional look-up table, which is the function of impact velocity and deflection; h₁, h₂ and h₃ 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.exp−h3vj(j=2,4)(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).

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 v_fdg is the friction draft gear working velocity.

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

Ffdg,i(xt,vt)=flxt,vt,xt≥xt−Δtfuxt,vt,xt<xt−Δt(5)

Ffdg,ixt,vt=Ffdg,ixt−Δt,vt−Δt+klocked(xt−xt−Δt(6)

fuxt,vt≤Ffdg,i(xt,vt)≤flxt,vt(7)

where x_t, x_t−Δt, v_t, v_t−Δt are draft gear deflections and relative velocities of adjacent vehicles of the current and previous time-step, respectively. f_l and f_u are force-displacement characteristics for loading and unloading processes, respectively. k_locked 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.

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

where P_i is the preload, b_i is the damping ratio (Nsm), k_i3 is the linear stiffness Nm, k_i2 is the nonlinear stiffness Nm, ν is the critical velocity constant, k_i1 is the recoil stiffness Nm, and Δx_i, Δx^._i are relative displacement and velocity of the wagons, ϵ is the coupling clearance.

The terms k_i1 and k_i3 are the main cause of hysteresis behaviour. k_i2 moreover, ν are responsible for the force peaking that can be derived both from experimental and theoretical analysis of draft gears [28].

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

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

where F_p,loco and F_p,i are the propulsion resistances which are represented as the sum of the aerodynamic and rolling resistances (kN); F_crv,loco and F_crv,i are the resistance force elements originating from the track curvature (kN); F_grv,loco and F_grv,i are the gravitational force elements due to the track gradient (kN); F_brake,loco and F_brake,i are the braking resistance forces caused by pneumatic braking (kN), respectively.

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

where A(N), 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).

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

where h is the dimensionless parameter, which changes from 500 to 1200 depending upon the railway vehicle; R_yis the radius of the curved track (m).

The complex model has seven-dimensional parameter space such as [h₁, h₂, h₃, k_locked, k_m, k_r, S] and four different stages. The simplified model has five-dimensional parameter space such as [b, k₁, k₂, k₃, ν] moreover, two different stages.

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

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

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

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;

1. Maximum adhesion force

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

where ẋ₁ is the locomotive speed kmh; m_stc = n.m_wagon,i + m_stc,loco is the static mass of the locomotive&wagon system (kg).

The coefficient of adhesion is determined by several parameters; namely the material properties, rail surface, wheel profile, axle pressures and vertical movements. According to UIC standards, the maximum value of the coefficient of adhesion is 0.33. It is also shown by experimental methods that the ratio of adhesion is dependent on the speed. For this reason, the “Curtius-Kniffler” formulation, which varies with speed, is chosen to be utilised for determining the adhesion coefficient. The maximum traction power (P_{τ, max}inkW)is defined by using a maximum traction force (F_{τ, max} = m_dyn.a_max + F_r), 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).

where a_max 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 (R_y = 300 m) from t = 75 s. to = 185 s. At t = 155 − 225 s, the system is climbing %ffl 5 slope with a full load. The operational speed at maximum motor power (v_Pmax) was considered as 16.72 ms. Maximum allowable operational acceleration (a_max) 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.Tsz−1+D.N1+N.Ts.zz−1) 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 T_s = 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.

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 space	System identification results (Ts=1e-4)
Permutation	FPE	MSE	Fit to estimation (%)	Function evaluation #	Cost	Norm of step	First-order optimality
LLLL	9.95e-7	8.15e-16	99.5%	172	8.15e-16	3.63e-12	0,4580
LLLH	6.43e-17	2.94e-21	99.2%	532	8.05e-22	2.61e-20	0,0691
LLHL	3.87e-7	3.09e-16	99.3%	260	3.09e-16	1.24e-17	0,8690
LLHH	9.20e-16	7.23e-21	99.6%	463	2.01e-21	6.54e-18	0,4670
LHLL	1.23e-15	5.62e-21	99.7%	677	1.60e-21	1.03e-21	0,0238
LHLH	1.37e-18	3.85e-21	99.5%	961	8.66e-22	1.51e-21	0.0525
LHHL	1.24e-17	2.52e-21	99.9%	659	6.19e-22	1.67e-21	0,0750
LHHH	4.91e-16	6.98e-21	99.2%	451	1.96e-21	2.09e-21	0,1690
HLLL	3.57e-18	5.41e-21	99.3%	567	1.46e-21	1.05e-21	0,0516
HLLH	8.12e-7	2.87e-16	99.5%	204	2.87e-16	6.28e-16	0,5700
HLHL	6.63e-7	3.28e-16	99.3%	182	3.28e-16	7.22e-12	0,7100
HLHH	7.26e-17	4.38e-21	99.6%	394	1.04e-21	6.22e-14	0,0437
HHLL	3.87e-20	3.88e-21	99.2%	1028	8.53e-22	2.50e-21	0,0199
HHLH	3.93e-7	3.97e-16	99.4%	299	3.96e-16	5.13e-17	0,7170
HHHL	5.12e-18	5.95e-21	99.7%	755	1.60e-21	4.11e-21	0,0962
HHHH	2.17e-15	4.93e-20	99.3%	556	9.28e-21	1.32e-21	0,1970

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

Table 4

Time-domain characteristics of the estimated transfer functions

Stiffness & damping parameter space	Time-domain characteristics (input: locomotive force, output: acceleration of first wagon) (Ts=1e-4) (unit step response)							Ramp response
Permutation	Rise time	Settling time	Overshoot (%)	Peak time	Peak	Settling min.	Settling max.	Integral of error (ramp steady-state error)
LLLL	0,0984	16,0424	10,0458	0,8856	1,10	0,8924	1,1005	0,0241
LLLH	0,0630	9,6578	8.65e-5	0,8985	1,00	1,0000	1,0000	0,0000
LLHL	0,0908	14,1356	6,9475	0,8993	1,07	0,9337	1,0695	0,0233
LLHH	0,0717	16,0063	2.89e-4	0,8996	1,00	1,0000	1,0000	0,0895
LHLL	0,0622	12,9081	3.45e-4	0,8866	1,00	1,0000	1,0000	0,0000
LHLH	0,0629	11,9214	1.24e-5	0,8962	1,00	1,0000	1,0000	0,0889
LHHL	0, 0562	9, 3024	3.85e-5	0, 8959	1, 00	1, 0000	1, 0000	0, 0000
LHHH	0,0724	16,0856	2.08e-4	0,9029	1,00	1,0000	1,0000	0,0895
HLLL	0,0646	13,3827	1.94e-5	0,8873	1,00	1,0000	1,0000	0,0890
HLLH	0,0990	18,1009	9,4985	0,8846	1,10	0,0990	1,0950	0,0240
HLHL	0,0880	13,4079	8,7312	0,9025	1,09	0,9094	1,0873	0,0261
HLHH	0,0648	12,0161	8.73e-5	0,9011	1,00	1,0000	1,0000	0,0000
HHLL	0,0571	9,3533	2.10e-6	0,8840	1,00	1,0000	1,0000	0,1595
HHLH	0,0838	11,6350	6,8591	0,8860	1,07	0,9303	1,0686	0,0241
HHHL	0,0714	16,2612	2.27e-5	0,8977	1,00	1,0000	1,0000	0,0893
HHHH	0,0856	54,2275	3.87e-4	1,9809	1,00	1,0000	1,0000	0,0899

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

After running the (LHHL) case in the simulation, the position difference (Δx_i = x₁ − x_i+1, for i = 1, …, 9) moreover, the velocity difference (Δx^._i = x^.₁ − x^._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 number	Euclidean norm of coupler force (sum(abs(X). ^2)^ (1/2))
	Acceleration phase	Curving phase	Slope climbing phase	Deceleration phase
1	1,0881 e7	6,5859 e6	3,3520 e6	1,7959 e7
2	8,1790 e6	6,5733 e6	3,3167 e6	1,8004 e7
3	7,6021 e6	6,5628 e6	3,2948 e6	1,8046 e7
4	7,5245 e6	6,5538 e6	3,2699 e6	1,8085 e7
5	7,2580 e6	6,5479 e6	3,2549 e6	1,8123 e7
6	7,1546 e6	6,5472 e6	3,2426 e6	1,8158 e7
7	8,4859 e6	6,5412 e6	3,2293 e6	1,8191 e7
8	1,0050 e7	6,5153 e6	3,2854 e6	1,8222 e7
9	8,5480 e6	6,7601 e6	3,6738 e6	1,8252 e7
Coupler Number	Infinite Norm of Coupler Force (max(abs(X)))
	Acceleration Phase	Curving Phase	Slope Climbing Phase	Deceleration Phase
1	2,0407 e5	8,1541 e4	6,4030 e4	6,6105 e4
2	1,1723 e5	7,4836 e4	5,5003 e4	6,6070 e4
3	1,0271 e5	7,0911 e4	4,9925 e4	6,6724 e4
4	1,0200 e5	6,8055 e4	4,6354 e4	6,6999 e4
5	9,5294 e4	6,5795 e4	4,3614 e4	6,7428 e4
6	9,5384 e4	6,3923 e4	4,1412 e4	6,8059 e4
7	2,0198 e5	6,2410 e4	4,5772 e4	6,8313 e4
8	2,2348 e5	6,3988 e4	4,4586 e4	7,0401 e4
9	1,4802 e5	8,7190 e4	6,6110 e4	7,1669 e4
Wagon number	Euclidean norm of wagon acceleration (sum(abs(X).^2)^(1/2))
	Acceleration phase	Curving phase	Slope climbing phase	Deceleration phase
1	103, 6771	73,1766	37,2448	199,5401
2	87,7673	73,0369	36,8519	200,0473
3	83,2560	72,9204	36,6084	200,5109
4	82,5249	72,8202	36,3326	200,9500
5	80,1823	72,7540	36,1655	201,3623
6	79,0968	72,7463	36,0284	201,7515
7	80,4709	72,6796	35,8807	202,1201
8	78,9832	72,3927	36,5045	202,4688
9	80,0325	75, 1122	40, 8194	202, 8003
Wagon number	Infinite norm of wagon acceleration ((max(abs(X))))
	Acceleration phase	Curving phase	Slope climbing phase	Deceleration phase
1	0,9759	0,9060	0,7114	0,7345
2	0,9511	0,8315	0,6111	0,7341
3	0,9448	0,7879	0,5547	0,7414
4	0,9445	0,7562	0,5150	0,7444
5	0,9416	0,7311	0,4846	0,7492
6	0,9414	0,7103	0,4601	0,7562
7	0,9748	0,6934	0,5086	0,7590
8	0, 9760	0,7110	0,4954	0,7822
9	0,9646	0, 9688	0, 7346	0, 7963
	Euclidean norm of locomotive force (sum(abs(X).^2)^(1/2))
Acceleration phase	Curving phase	Slope climbing phase		Deceleration phase
1,2907 e8	6,2586 e7	4,3145 e7		1,7524 e8
	Infinite norm of locomotive force ((max(abs(X))))
Acceleration phase	Curving phase	Slope climbing phase		Deceleration phase
7,3159 e5	3,5971 e5	1,6093 e5		6,1303 e5
	Euclidean norm of controller error (sum(abs(X).^2)^(1/2))
Acceleration phase	Curving phase	Slope climbing phase		Deceleration phase
88, 1869	30,1110	29,4740		43,0090
	Infinite norm of controller error ((max(abs(X))))
Acceleration phase	Curving phase	Slope climbing phase		Deceleration phase
1, 1752	0,7748	0,9014		0,6665

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

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

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

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 k_i1 to check the lower frequency range, whereas it is taken as k_i2 + k_i3 to consider the upper limits for the second case (Δx^._i > 0).

The natural frequency results are shown below in Figure 8.

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 k_i1 moreover, high stiffness k_i2 + k_i3 cases, respectively. In other words, the dynamic response of the acceleration phase includes both low and high stiffness behaviour. Similar deductions for the other phases can also be made. Table 6 summarises these significant inferences.

Table 6

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

Phase type	Coupler 1	Coupler 8	Coupler 9
Acceleration	Two distinct peaks, botd low and high stiffness modes	Two distinct peaks, botd low and high stiffness modes	One single peak, low stiffness mode
Curving	Several distinct peaks, botd low and high stiffness modes	Two distinct peaks, botd low and high stiffness modes	One single peak, low stiffness mode
Deceleration	One single peak, low stiffness mode	One single peak, low stiffness mode	One single peak, low stiffness mode
Slope climbing	Several distinct peaks, botd low and high stiffness modes	Two distinct peaks, botd low and high stiffness modes	One single peak, low stiffness mode

Table 7

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

Phase type	Coupler 1	Coupler 8	Coupler 9
Acceleration	Max: 9691 N	Max: 9481 N	Max: 8337 N
	RMS: 4066 N	RMS: 3519 N	RMS: 2792 N
Curving	Max: 2190 N	Max: 2544 N	Max: 2700 N
	RMS: 1010 N	RMS: 965 N	RMS: 1261 N
Deceleration	Max: 6490 N	Max: 8617 N	Max: 8617 N
	RMS: 1408 N	RMS: 1957 N	RMS: 2000 N
Slope climbing	Max: 2090 N	Max: 2785 N	Max: 2860 N
	RMS: 913 N	RMS: 898 N	RMS: 1162 N

6 Limitations and discussion

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

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