Dynamics of railway track subjected to distributed and local out-of-round wheels

– The dynamic response of a railway track under a moving train in the presence of a defective wheel is studied. The goal of this study is to establish an eﬃcient dynamic model for wheel tread defect analysis. The wheel irregularities which inﬂuence the track vibration characteristics are classiﬁed in 2 major types: periodic out-of-roundness and local defect. While the periodic irregularities are distributed all around the wheel, the local defects such as wheelﬂats have more impact on the dynamic response. The inﬂuence of each defect depends on its geometrical characteristics and the model parameters. The separate eﬀects of diﬀerent irregularities have been the subject of many previous researches. However the contribution of each defect to a coupling of irregularities was not as much considered. For a comprehensive analysis of all types of defects, an eﬃcient global model is presented. It includes 3 sub-models for vehicle, contact force and track. The time-domain responses of the model are evaluated in comparison to those of a ﬁnite element model. It is shown that the semi-analytical model of the track keeps a good precision while being less time consuming. The results of the global model for diﬀerent types of wheel irregularities are presented.


Introduction
One of the most important problems facing the railway maintenance is the analysis of dynamic behavior of tracks under moving loads. In recent years, the speed and load of commercial trains have been increased significantly. The structures are therefore subjected to intense vibrations and dynamic stresses, which in turn are much higher than the corresponding static stresses. In this case, the wheel imperfections highly influence the interaction force which can cause severe vibrations in the railway equipments and so damage both vehicle and track. The main geometrical defects are wheel roughness and profile out-of-roundness (OOR). The OOR defects concern the deviation of the wheel tread geometry from its circular shape. The two major types of OOR are periodically distributed and local. The former is characterized with one or more wavelengths on one turn of the wheel and causes a harmonic excitation to the track, whereas the latter causes an impulse excitation with larger amplitude. The distributed defect is in general a non homogeneous wear gradually created around the circumference of a Corresponding author: antoine.dequidt@univ-valenciennes.fr the wheel [1]. Some experimental measurements on worn wheels have shown that the profile irregularities generally comprise both types of defects. According to Nielsen et al. [2], the wheel OOR is composed of a set of periodic waveforms with corresponding wavelengths. The number of waveforms called OOR order and their wavelengths are related to the train load and speed. The OOR order varies from 1-3 for high-speed trains to 1-11 for freight trains. An eccentric wheel will induce the frequency corresponding to the number of wheel revolutions per second (a wheel of radius R gives an irregularity of wavelength λ = 2πR). An oval wheel will induce twice that frequency (wavelength λ/2), and so on [2].
Local (non-periodic) OOR may appear on the wheel tread, especially if the wheel is equipped with block brakes. The local defect is specifically a wheelflat in 3 stages. Wheelflat is formed when the brakes are applied to a railway wheel with insufficient friction between wheel and rail, so that the wheel locks and slides along the rail. This sliding causes severe wear of the part of the wheel which is in contact with the rail called wheelflat. At first, a small part of the wheel tread is flattened which occurs a newly formed flat (NFF). Soon after the sharp vertices are rounded by wear and a kind of partially round flat (PRF) is formed. The flat's whole length is gradually worn and a completely round flat (CRF) is finally generated. For the same wheelflat with a fixed depth, the CRF length is estimated to be 1.5 times of the NFF length, and the PRF length varies between them. While the initial length is typically between 10 and 50 mm, it can extend to over 100 mm. In wheel maintenance system (WMS) for tourist and freight trains, the wheels containing the flats with up to 40 mm of length are allowed to work. However the impact forces excited by wheelflat change not only with its length, but also with its depth, train speed and static load [3].
In order to investigate the effect of OOR with different forms on track-train dynamic interaction, a proper mathematical model is essential. Such a model can be used for analysis and numerical simulation for different train characteristics and parameters (speed, load . . . ). The aim of the present work is to prepare a valid time-domain model to distinguish the influence of different wheel defects on railroad track. The essential part which determines the accuracy and speed of this model is track. In Section 2.1 a semi-analytical model of discretely-supported Euler-Bernoulli beam with crosslinked ballast is used. Based on modal superposition method, an efficient track model with reduced calculation time is introduced. This approach has been used for a finite track [4][5][6] and the assumed number of modes should correspond to the limited number of spans. To compare and validate the results, a finite element model is used with the same conditions and simulation parameters. The finite element approach allows optimizing different characteristics of the semianalytical model as well.
The defects are introduced on the front wheel profile of a half body vehicle model with 5 degrees-of-freedom presented in Section 2.2. First by introducing a circular wheel profile, the contact force between the wheel and the rail is formulated. Then a non-Hertzian model is applied by using a Winkler bedding method. The results are evaluated by the Hertzian force of interaction for a healthy wheel, i.e. a convenient condition of contact. However in the case of OOR wheels the hypotheses of the Hertz model are not fulfilled. Comparison of dynamic force obtained by the Hertzian contact model with more realistic non-Hertzian models showed that the Hertzian impact force for a wheelflat goes beyond the non-Hertzian ones [3,7,8].
Both types of OOR have been extensively studied form different viewpoints like appearance, mathematical model and dynamic responses. Study of wheels containing distributed polygonal OOR with a few dominating harmonics on high-speed trains have been previously addressed in [9]. Several wheel tread defects were investigated by the peak impact forces and peak rail bending moments and simulation results were validated by measured data [10]. Different authors emphasized more directly the importance of wheelflat in the context of dynamic wheel-rail interaction and track deterioration [7,[11][12][13][14]. The effect of presence of a second wheelflat on impact force of one other wheelflat situated on same wheel or on the other one of a bogie was studied [15]. Then the role of wheelflats length and depth, phase difference and vehicle and track parameters on peak of impact force was given [16].
The third section introduces a novel equation for different types of OOR. Based on some experimental cases of OOR, the corresponding wheel tread geometries are defined. Then the non-Hertzian model of contact is established between rail and wheel. Section four gives the results of the time-domain model when the track is subjected to different defective wheels. The importance of moving speed of wheel is studied subsequently.

Global model of wheel-track
The model of rail and its supports is the most important part which affects the accuracy and speed of the simulation and has attracted much attention. The simplest model which is widely used to represent an infinite rail is Euler-Bernoulli beam [5,[16][17][18]. While continuously supported beam may be used only up to about 1000 Hz, the discretely supported rail model could consider the vertical vibration from about 1000 to 6500 Hz [19]. This type of model is a flexible approach which could be developed in different degrees of freedom for the support. In [20] the three layers of discrete springs and dampers represent the elasticity and damping effects of the rail pads, the ballast, and the subgrade with the crosslinked ballast. This approach is selected for the track model because the parameters of track substructure including mass of sleeper and ballast, stiffness and damping of railpad, ballast and subgarde also shear stiffness and damping of ballast are separately defined. A beam of length L which is pinned up in 2 ends and supported as described in [20] is shown in Figure 1. The separate masses of sleeper, ballast and subgrade with their visco-elastic connection permit a precise identification of model parameters in an experimental platform. The influence of displacement at each point on the other elements of the model could be considered as well.
The equation of motion in matrix form is achieved by using Lagrange's equations and the assumed mode method [6]. In the modal approach, the vertical displacement of a point located through the longitudinal coordinate x at the instant of t is where φ i (x) is the ith vibration mode and q i (t) is the modal coordinates. In assumed mode method, the problem size depends on the number of intermediate points (m) and the number of vibration modes considered in the modal descriptions of these spans (n) [8]. To satisfy the boundary conditions at the two ends of the beam, the assumed function is as following   The kinetic energy of the track (rail and supports), the potential energy of the rail, the potential energy and the dissipation function of m supports used in Lagrange's equations are written in matrix form as where Q andQ are n × 1 vectors containing q i andq i , Z s andŻ s are m × 1 vectors of sleeper vertical positions and speeds, and Z b andŻ b are m × 1 vectors of ballast vertical positions and speeds. The beam mass per unit length, m r , equals to ρ × S. All parameters of the track comprising the masses and the stiffness and the damping coefficients of the elements shown in Figure 1, are given in Table 1, similarly with [20]. The matrices H, M and ψ are introduced in Appendix A. The matrices Δ and K are defined by equations (7) and (8) The ODEs of the system which gives Q, Z s and Z b vectors are d dt where P ψ is the generalized force done by the moving load (P ), L is the Lagrangian function of the track, and is defined by equation (12) If L and D are replaced in equations (9)-(11), the modal coordinates vector of the beam displacement is obtained by the following set of governing equations The numerical integration is performed by the forth order Runge-Kutta method. Afterwards the equation (1) gives z r (x, t) in each point of the rail while the load moves. The simulation uses the material properties of UIC60 as given in Table 1.
In the results shown in Figure 2, it is considered that the beam contains 15 spans, i.e. number of sleepers is 14. The Figure diagrams demonstrate the displacement under a constant load on the whole (a) and five central spans (b) of the beam. The total simulation time could be minimized by limiting the number of assumed modes (n). However insufficient number of modes affects the precision of dynamic response. It is shown that the beam stiffness in a truncated model is overvalued. So there is a difference between deflection results of the beam in this method due to the truncation mode. Some correction methods were proposed to remove the effect of mode-superposition [21]. In order to find the smallest number of modes insuring the convergence of the model, Figure 2 illustrates the effect of increasing number of modes. The results are given for a static load of 250 kN, relating to the maximum load resulting from the permissible axle load, that travels at 30 km.h −1 . Some integer multiple of spans is chosen for mode assumption. Figure 2 shows that in this method three times of the spans number is appropriate associated modes to represent the dynamics of the system and for this track comprising 15 spans, the satisfactory number of modes is 45.

Validation against finite element method
Using finite element method for the discretely supported track model allows improved prediction of dynamic response and offers the potential for refinement by including all conceivable track components as layers. The advantages of finite element modeling is that the dynamic analysis solved by use of numerical time-stepping routine, non-linear components and contacts may be included, using real valued modal analysis and U-transformation, to find the dynamic response of an infinite uniform beam, by exact or approached methods.
The same track with similar condition and parameters is modeled by a finite element method as like in [22] presented in Figure 3a. Each beam section is divided into n elt elements. To evaluate the vertical displacement of the beam, 2 degrees of freedom of each node are considered: vertical translation and planar rotation. To achieve the dynamic response of the beam under moving load, the Newmark's implicit method is used. Figure 3b indicates that for different number of elements/span, the displacement obtained by this method coincides with the one achieved by the semi-analytical method, whereas a reduction in element size (increasing n elt ) leads to results fit of both methods. But a decrease in element length leads to an increase of the same order in calculation time; while the necessary time for semi-analytical method to give the same response is much shorter.

Interaction force of wheel-track contact
The vehicle motion is governed by the wheel-rail contact. A simple way to compute the vertical interaction between vehicle wheels and the track is to treat independently each wheel by using the superposition principle considering the rail vibration as a frequency band average [13]. In most cases a single unsprung mass on a nonlinear Hertzian spring represents the interaction model [11,18,23]. However this model may largely underestimate the dynamic response and a 2D vehicle model containing 2 unsprung masses which forms a 5 degreesof-freedom model is chosen as in [24]. These degrees of freedom comprise two degrees of vertical translation for the wheels, two degrees of translation and rotation for the bogie frame and one degree of vertical translation for the half vehicle body.
The wheels with the same characteristics are connected to a bogie frame via a primary suspension system. The frame is connected beyond a half vehicle body through a secondary suspension and the suspension characteristics are chosen as given in [20].
Despite the popularity of 2D Hertzian contact [19,20,[25][26][27][28], limits of this approach are that geometric requirement in Hertzian model by which the non-deformed surfaces collide, should be elliptic paraboloids. Due to the other assumptions for the contact, the requirements are not satisfied for defective wheels. Thereby some non-Hertzian models of contact are developed in [6-8, 10, 29].
A multipoint contact type used for interaction model between wheel and rail is a distributed bedding which includes some independent springs [19,27,28]. The length of the bedding should be large enough to consider all possible sizes of the interaction zone. It evaluates the normal contact force by the sum of individual forces of the springs. A minimum of spring number is necessary to guarantee a good assumption. The model is two dimensional and the width of contact is supposed constant. Each force takes into account the spring stiffness and the local distance between the interfering surfaces. Assume that the stiffness for all the springs is identical, if n c springs are uniformly distributed on the bedding of 2δ length, as shown in Figure 4a, the common stiffness of the springs is k i = 2δK W /n c where 1 ≤ i ≤ n c and K W is the stiffness per unit length of the springs. For an arbitrary wheel profile and rail surface, the distance of surfaces on each spring position in an instant of time is based on the local wheel and rail displacements ((z p ) i and (z r ) i ).
For ideally round wheel and smooth rail surface, the necessary conditions for the Hertzian model are available. Figure 4b shows that the results of the bedding model for contact force and consequently the beam displacement coincide the Hertzian model. The vehicle moves at a low speed of 30 km.h −1 and the displacement is registered on 5 spans at the middle of the track to receive a complete turn of the wheel.

A new formulation of irregularities on wheel profile
The force of contact is expressed as a function of the relative displacement between wheel and rail at the contact zone. It depends on the un-deformed wheel/rail geometry and the elastic characteristics of the wheel-rail contact. In the present method, a wheel rolls as explained on the track. The moving wheel geometry model is therefore essential to investigate the effects on wheel/rail interaction due to the parametric excitation. The irregularity function is expressed as follows to determine the wheel geometry.

Distributed defects with periodic description
Every distributed defect of a wheel could be expressed by a set of periodic functions based on Fourier Transform. To simplify and confine the number of equations, a limited number of amplitude and frequency (more commonly wavelength) of irregularity is considered for the mathematical model. The number of periodical terms called wheel orders (n w.o. ) is evaluated from 1 to 20 [30]. Each frequency (f i ) is defined as the repetition number of an irregularity and beside its half amplitude (A i ) and phase lag (χ i ), determines the form of a distributed defect. Hereby the periodic irregularity function is related to the wheel rotation angle (ψ) by  where 1 ≤ i ≤ n w.o. . In Figure 5, a typical distributed defect is decomposed into 5 harmonic waves with amplitudes of 0.25, 0.20, 0.15, 0.05 and 0.01 mm respectively. Considering the difference between periodicity and phase lag of each wave, the amplitude of the irregularity varies between 0.1 and 0.5 mm.

Discrete OOR in wheelflat form
The geometry and important role of a shell or a wheelflat as the principal types of local defect on the dynamic interaction was previously detailed. However the coupling effect of different OOR (regardless distributed or local) which is more objective is less investigated. Before studying some coupled cases, progressive formation of a local flat defect in 3 steps is presented as shown in Figure 6a. First appearance of a flat on wheel profile is in the form of NFF. This type of flat is characterized by its sharp edges which are rapidly rounded because of the wear. Hereupon the second type of flat which is PRF is formed. The wheelflat length increases gradually while its depth, d, continues to be constant until the creation of CRF. The flat length in NF, CR and PR forms related to d using could be found respectively by [7] l NFF = 8Rd − 4d 2 (17) where R is the wheel radius. The presence of these 3 types of a wheelflat with depth of 7.4 mm is represented in Figure 6b based on the above formulations for R = 460 mm.
To establish the contact force of each wheelflat, the bedding model uses the interactive surfaces distance, Δz i = (z p ) i − (z r ) i for 1 ≤ i ≤ n c , in terms of a discrete variable ξ i defined by where −δ ≤ ξ i ≤ δ and 1 ≤ i ≤ n c . Each defect is positioned on the wheel tread by its phase lag (χ F ) and center  angle (Φ F ) given in [12]. The wheel angle on ith contact point is evaluated by For ψ i more than 2π, it should be brought back to [0.2π].
In Figure 7a, the length and the depth of the flat are simply represented by l F and d respectively. For all flat types the following statement is proposed to calculate the local irregularity of wheel profile, shown by (LI i ) in Figure 7b relating each point of contact The superposition of distributed and local irregularities beside the vertical position of one of the wheel center gives where 1 ≤ i ≤ n c and z w is the height of the wheel center. Therefore the contact distance Δz i is calculated to achieve the normal force.

Simulation results
The dependence of contact force on wheel tread profile induces to study how each defect influences singularly the dynamic response of the track, which role the characteristic parameters have and do the detail of all defects is interpretable in an output signal of experimental instruments providing track deflection, acceleration or pressing force. The last case will be the subject of another survey. In this part the effect of an irregularity relating its geometry and some combined defects are presented.

Parameters of periodic irregularity
The major parameters of this type of defect are the periodicity of irregularity and its amplitude. For different irregularity periods and a constant amplitude of 0.1 mm, Figure 8a illustrates the dimensionless force applied on the front vehicle wheel which travels at 120 km.h −1 . The dimensionless force is normalized by dividing the contact force to the total static load. The defect is just allocated to the front wheel and its repetition is assumed at 1, 2, 4, 8 and 12 times per round. It is seen that for higher periodicities, the force takes more values and oscillates faster. However a change in amplitude affects severely the interaction force and so the track and wheel displacement. If a runout type OOR with the wavelength of 2πR with variable amplitudes ranging from 0.05 mm to 2 mm degrades the front wheel, the maximum dynamic force changes from 1.02 to 1.35 times of static load, as shown in Figure 8b.

Effect of wheelflat geometry
The depth of a wheelflat holds the main role in equation (22) to obtain the contact force and so the track displacement. Additionally the type of wheelflat and in the case of multiple defects, their arrangement on the wheel is the other important factors. A depth of 3 mm is chosen for 3 types of wheelflat and the length of PRF is taken as the mean of NFF and CRF. In the case of Figure 9a where NF, PR and CR Flats are positioned at the same distance in 90 • , 210 • and 330 • respectively (sequence 1), the dimensionless force diagram shows that the peak of CRF dominates lightly the two other types. In sequence 2, the positions of NFF and CRF are interchanged and the contact force against NFF and PRF remains unchanged, while against CRF experiences falls down for about 10%. It reveals the rather independence in the response under 3 flats with the lengths of 105, 135 and 165 mm with approximately 830 mm of interval.
The main effect of wheelflat on the vehicle parameters concerns the critical speed, defined as the speed in which the wheel separates from the rail [25]. At the separation moment, the normal force tends to zero and just afterwards it rises up to a peak value. The strength of the impact and the frequency of repetition are proportional to the train speed. The wheelflat causes the critical speed to reduce much lower than general conditions. To study the speed effect, a defective wheel as in sequence 1 is used and passed at 3 current speeds of freight trains: 80, 120 and 160 km.h −1 . In Figure 9b, the effect of the wheelflat is shown to aggravate when the wheel speed increases. While the dimensionless force never exceeds 5 at 80 km.h −1 , its maximum for PRF and CRF passes 7 when the speed is twice. Although the separation moment on CRF is about 2 times of the one on NFF, in all of the flat cases this period is slightly larger for higher speeds. The peak of track displacement is lower when the speed increases, because the wheel passes sooner over the flat. However for higher speed and larger flat, the track displacement varies over a longer period of time. Here also the effect of presence of each flat on the others is relatively negligible since the contact forces converge at the beginning of the next flat.

Mixed irregularities cases
The classification of wheel irregularities into distributed, local and a mix of them covers almost all defects including stochastic types. When some different irregularities are coupled, a new wheel tread is resulted which could be a mix of distributed periodic defects, local wheelflats or both of them. In an optimized algorithm, the number of constructing terms should be minimized. Based on equation (16) for n w.o. = 5, an OOR composed of just periodic terms is presented by equation (24). The amplitudes are given in mm and Figure 10a indicates that the irregularity range around the wheel tread varies between 0.3 and 0.55 mm. However this range remains marginally beyond 0.6 mm which is the maximum defect amplitude allowed in the WMS for tourist and freight trains.  Figure 10b illustrates the under front wheel displacement of the track which is caused by this wheel profile geometry. The train travels at 120 km.h −1 and the interaction force is shown to change between 0.95 and 1.07 time of the total static load. It means that for such a defect with marginal tolerance, the contact force and so the track displacement at a usual speed varies below 10% of normal values which happens when a healthy wheel passes on the track at the same speed Another kind of coupled OOR is built up by superposing two or more local defects in one of wheelflat shapes. This type of irregularity is mainly characterized by its abrupt change on the dynamic response. Two major examples are given in Figure 11a, where 2 CRF with the same geometry are combined in two further positions of the front wheel tread. The depth of the flats is considered 1 mm and their length is 95 mm. At the first angular position of 90 • , as the first CRF finishes, the second one starts and an ω-form defect is formed on the wheel profile. At the second wheel position where the angle is 270 • , the second flat is placed inside the first one. Therefore an u-form OOR is created which resembles the concave shell defects. While the former has the same depth with a double length relating the original CRF, the latter is deeper and wider for about 0.4 times.
If no other defect is situated on the wheels and the train travels at the speed of 120 km.h −1 as before, the results of track displacement and dimensionless force will be as presented in Figure 11a. In order to give a comparison idea relating the presence of a one single wheelflat in each position, the relating diagrams are reviewed in this Figure. At the first position, the force diagram of the coupled case follows completely the single wheelflat till the end and after then it responses under the second part of the defect with an obviously higher peak (3.8 versus 3.4). So the residual oscillation of the first part affects clearly the rest of the defect and the separation takes places twice. The second position where the total length less than the previous one, the force curve is the result is achieved through a superposition and the wheel separates the track for one time. It is seen that this force is lower than the force for the single wheelflat.
Lastly the coupling between the mixed distributed and local defects is cited which could result also to stochastic defect generation. On a harmonic defect as presented in Figure 10, two CR flats with 3 different depths are superposed at upper and lower situations of the wheel tread to see the role of each element on the dynamic response. The CRF depths are 1, 0.5 and 0.1 mm and their lengths are respectively 95, 67 and 30 mm. As before, on the distributed part of the coupled defect, the displacement and force follow a round wheel response with a tight variation. On the other hand, at local defect section a sudden change in interaction force is observed even for the smallest wheelflat whose depth is comparable with the difference between upper and lower amplitude of the harmonic defect. Figure 11b demonstrates that for the defective wheel which travels at 120 km.h −1 , the contact force under the flats with the depth of 1, 0.5 and 0.1 mm situated on lower distributed irregularity zone, reaches 3.2, 2.35 and 1.4 times of the total static load respectively. These values for higher distributed irregularity zone reach 3.8, 2.55 and 1.5 times of the static load respectively.
By following the corresponding curves on the track displacement diagram, the same event is interpreted differently. At both positions with upper and lower amplitudes, the flat with 1 mm of depth is separated completely the track for a while and the track displacement tends rapidly toward zero. The return of the wheel applies a high shock to the track and makes it vibrate severely. The contact loss is happened for the flat with 0.5 mm of depth but during a shorter period. Therefore the applied shock and vibration is less considerable. The smallest CRF never separates the wheel from the rail and so the displacement variation continues very slightly after that the flat passes.

Conclusion
This article presented a model for track, train and their contact to study the dynamic response of the system for different types of wheel defect. An assumed-mode model of discretely-supported beam was used and validated by the finite element method. The defective wheel has been positioned at front section of a half vehicle and a Winkler contact model for wheel/track interaction is formulated. Considering the importance of wheel profile in railway issue, large variety of irregularities was categorized in 3 geometrical forms: harmonic distributed, local wheelflat and mix of both.
First the distributed defect has been expressed by a limited number of periodic functions, where the amplitude and frequency of irregularity were determined. Because of the important role of wheelflat in tread OOR, its progressive forming was then introduced by NFF, PRF and CRF, as the principal types of local defect. The coupling effect of different OOR was the other subject of this study and 3 cases with purely periodic irregularities, coupling of two wheelflats and an arbitrary integration of CRF on periodic OOR. By using this new method of defect description, the dynamic responses were represented in the last part of the article.
The results of traveling a wheel containing periodically distributed OOR showed the importance of both amplitude and periodicity. When the defect amplitude exceeds 1 mm, the generated overload could reach more than 10%. However the effect of periodicity on the contact force while it remains beyond four is not much considerable. The presence of multiple flats of different types on further positions of wheel tread shows the independence of caused impact since the influence of this defect disappears soon.
The flat depth and the wheel speed have both remarkable effects on contact force and so track displacement, while the difference in the length following the sharp or rounded sides does not seem to be very important. The most different effect in the results revealed by this work is observed by mixing all the defects in an arbitrary form. The effect of distributed part of the OOR on the contact force differs slightly from the response of a round wheel, whereas local defect even with small size generates a large shock on the track. Therefore the considerable effect of each elementary defect could be detected following the presented base of separation. Lastly a high speed semi-analytical method is obtained to process the passing wheel with geometrical irregularities.
Acknowledgements. This work was developed in the framework of the Track Train System Availability (TTSA) project. This project is supported in part by the European Community (through the FEDER European Funds for Regional Development) and the Région Nord Pas-de-Calais (through i-Trans competitiveness cluster). This work was also supported by International Campus on Safety and Intermodality in Transportation, the Délégation Régionaleà la Recherche et a la Technologie, the Ministère de l'Enseignement supérieur et de la Recherche and the Centre National de la Recherche Scientifique (CNRS). The authors gratefully acknowledge the support of these institutions.

Appendix A
Here the matrices H and M are introduced as follows [5]: where ϕ i (r) = sin iπr L and ϕ i (r) denotes the second derivates of ϕ i (r) with respect to r. ψ is a vector defined as where x is the horizontal position of the load and so ψ should be updated as the load travels on the track.