1 Introduction

Rail corrugation has become one of the major problems in the field of railway engineering. The strong impact that this pathology has on the noise and vibration generation phenomenon affects users’ comfort and represents one of the main transport externalities.

In order to study in depth this pathology, several authors as Grassie and Kalousek [1] and Grassie [2, 3] analyzed from a theoretical point of view the generation of corrugation phenomenon and its possible treatments, concluding that all types of corrugation known were associated to the resonant behavior of the vehicle/track system. Other authors like Suda et al. [4] conducted an experimental study of the corrugation generation phenomenon on the top surface of a rail in a tight curved track. They concluded that at steady state, the longitudinal slip determines the corrugation type and its wavelength.

Further investigations, as the one developed by Jin et al. [5], studied the influence of rail corrugation in the vehicle–track dynamics through a calculation model of a half-passenger car coupled with a curved track. Conclusions showed the great influence of corrugation morphology in the wheel–rail dynamic response. Thus, the deeper the corrugation depth and the shorter the corrugation wavelength were, the greater the influence of corrugation was. Other authors, as Zhao et al. [6] or Torstensson and Nielsen [7], continued with this aim and presented simplified multi-body finite element (FE) models of damaged rails to study the wheel–rail interaction in the time domain. Thus, [6] studied the stress and strain states of the rails through a model composed of the primary suspension of vehicle, half locomotive wheelset, one rail, rail pads, and ballast demonstrating the influence of rail surface on the oscillations of contact forces and stresses. Likewise, [7] simulated a single bogie circulating through a corrugated curved track, demonstrating the influence of wheelset eigenmodes in the wheel–rail dynamic response.

The high loads generated in the wheel–rail contact depend on the vehicle speed, as Hawari and Murray [8] demonstrated through a multi-body interaction model in NUCARS software able to reproduce accurately the dynamic performance of the vehicle. This not only affects the severity of track damage, but also affects the track elements, as Ling et al. [9] demonstrated through a multi-body FE simulation. In this study, some of the main causes of the fracture of clips were performed, and several mitigation measures to reduce track vibrations were suggested. In fact, the great adaptability of numerical models has led to incorporate preventive and corrective maintenance measures in the calculations. Regarding the preventive measures, Collette et al. [10] developed a multi-body FE model capable to predict the effectiveness of a dynamic vibration absorber located at the vehicle axle for mitigating rutting corrugation. They concluded that this phenomenon was amplified when the vertical resonance of the vehicle track and the torsional resonance of the wheelset matched. Hence, acting on the torsional modes of wheelsets could reduce this phenomenon appearing. However, from the corrective measures’ point of view, Egaña et al. [11] studied the effectiveness of a liquid capable of modifying friction of a corrugated rail surface. They established that by applying this liquid, the accelerations generated in the wheel–rail contact and the corrugation growth were reduced.

The aim of this paper is to continue in this research field by analyzing the influence of rail grinding on the vibration generation phenomenon by a time domain feedback between a multi-body model of the vehicle and a FE model of the track. For this purpose, the dynamic forces resulting from the passage of the train through the track (defined as dynamic overloads) are calculated with a multi-body model and then incorporated to the FE model. For calculating of the dynamic overload, two cases are presented: the first one (case A) includes a track where the inner rail is corrugated, while the second one (case B) consists in a track without defects. Then, the influence of rail grinding will be evaluated by analyzing the different vibratory response of the FE model in both cases.

2 Data-gathering campaign

The studied stretch (Fig. 1) is a ballasted track of gauge S = 1.668 m, radius R = 300 m, and cant h = 0.027 m. A rail corrugation on the inner rail has been detected of wavelength λ c = 0.12 m and amplitude A c = 400 μm.

Fig. 1
figure 1

Studied stretch

Full size image

The existing superstructure is formed by an UIC-54 rail resting on a rail pad of 0.007-m thick. The sleepers are made of prestressed concrete of 0.3-m width separated every 0.6 m. The ballast layer is 0.5-m thick and it is resting on a large layer of green clays. The main characteristics of the materials that form the section are shown in Table 1, where E represents the Young Modulus, ν represents the Poisson’s coefficient, and ρ represents the density of each material.

Table 1 Main mechanical properties of the materials composing section
Full size table

The registered vehicle was a CAF S-599 (Fig. 2). This is a symmetric vehicle composed of three wagons and six bogies whose main characteristics are shown in Table 2.

Fig. 2
figure 2

CAF S-599

Full size image
Table 2 S-599 main characteristics
Full size table

For the data gathering, three PCB 354C02 accelerometers were disposed in the studied stretch. Their main characteristics are shown in Table 3.

Table 3 Main characteristics of the PCB 354C02 accelerometers
Full size table

With the data collected by the accelerometers, it was possible to characterize the real vibratory response of the track. Moreover, the vehicle speed was calculated using Eq. (1):

$$v = \frac{{d_{\text{e}} }}{{d_{\text{p}} }},$$
(1)

where v is the circulation speed, d e represents the distance between two consecutive axles, and d p represents the gap between peaks of the accelerogram. Thus, the train speed is set to 27 km/h.

3 Calculation method

In order to obtain the railway-induced vibrations, two different models are developed. In the first one (model 1), the track is developed by a numerical FE model in the time domain, while in the second one (model 2), the vehicle is developed with a multi-body model. Then, the vibratory response of the track is obtained by matching both models.

However, this interaction has a very significant drawback. One of the inputs necessary for its development is the equivalent stiffness of the track (\(K_{{{\text{eq}}_{n} }}\)), which cannot be calculated because some of the main properties of materials of the tack are undefined (Table 1). To solve this problem, the following calculation procedure is proposed (Fig. 3):

Fig. 3
figure 3

Calculation method

Full size image

Firstly, as the dynamic overloads caused by the rail corrugation are unknown, a simplified numerical model that disregards the presence of such pathology is developed (model 1a). Next, a first approximation to the dynamic overloads is obtained by experimental formulation and it is then applied to the numerical model. In this moment, a process of pre-calibration is carried out on the high rail (iteration n = 0), obtaining a first approximation of the equivalent stiffness of the track (\(K_{{{\text{eq}}_{0} }}\)).

Subsequently, a multi-body model is developed by VAMPIRE commercial software. This model takes into account the equivalent stiffness calculated above, the vehicle characteristics, and the mechanical and geometrical characteristics of the superstructure. The result is a dynamic overloads vector F n generated in the wheel–rail contact of the high rail (with no defects) and low rail (corrugated). Hence, this vector will depend on the calculated stiffness \(K_{{{\text{eq}}_{0} }}\).

Secondly, a more realistic numerical model is developed (model 1b) and the dynamic overloads vector F n is applied. From this point on, a new calibration process in both rails begins, and a new stiffness value (\(K_{{{\text{eq}}_{n} }}\)) is obtained. It is expected that, once the calibration process ends, the stiffness obtained in the model 1b will be different to that obtained in model 1a due to the new material properties.

The new stiffness value is used to recalculate the dynamic overloads, obtaining a new dynamic overloads vector (F n+1). This vector is applied on model 1b and a new vibrational response of the track is obtained. Then it is compared with the one previously calculated. If there is a relevant difference, the calibration process should be repeated until \(K_{{{\text{eq}}_{n - 1} }} \approx K_{{{\text{eq}}_{n} }}\). In that moment, the model will be calibrated and validated for a \(K_{{{\text{eq}}_{n} }}\) stiffness and a F n dynamic forces vector.

3.1 Numerical model

The FE model is developed using ANSYS LS-DYNA V14 software. It aims to solve the following equilibrium equation on all nodes of the model:

$${\varvec {M}}{\ddot{\varvec u}} + {\varvec{C}}{\dot{\varvec u}} + \varvec{Ku} = \varvec{f}^{\alpha } \left( t \right),$$
(2)

where \(\varvec{M}\) is the mass matrix, \(\varvec{C}\) is the damping matrix, \(\varvec{K}\) is the stiffness matrix, \(\varvec{u}\) is the displacement vector, \(\dot{\varvec{u}}\) is the velocity vector, and \(\user2{\ddot{\varvec u}}\) is the acceleration vector. Combining these elements, the external force vector f α (t) is then obtained.

For calculating the damping matrix C, the Rayleigh’s damping theory is applied as

$$\varvec{C} = \alpha \varvec{M} + \beta \varvec{K},$$
(3)

where α and β represent the Rayleigh’s damping coefficients.

As the predominant frequency in the numerical model corresponds to the loadsteps application [Eq. (4)], where v is the velocity of the vehicle and d s is the distance between points of load application (estimated in \(\lambda_{\text{c}} /2\)), the contribution of the Rayleigh coefficient α can be neglected according to Real et al. [12].

$$f = \frac{v}{{d_{\text{s}} }} = \frac{7.5}{0.06} = 125\;{\text{Hz}} = 785\;{\text{rad/s}} .$$
(4)

Hence, Eq. (5) is obtained by combining Eqs. (2) and (3):

$${\varvec {M}}{\ddot{\varvec u}} + \beta {\varvec {K}}{\dot{\varvec u}} + \varvec{Ku} = \varvec{f}^{\alpha } \left( t \right).$$
(5)

3.1.1 Meshing

According to [12], the huge calculation time needed to develop numerical models requires the adoption of certain simplifications which do not involve a loss of accuracy. Thus, it is assumed that all materials are homogenous and isotropic and have a linear elastic behavior.

Furthermore, the minimum frequency range to study is set from 2 to 50 Hz in order to include the frequency interval relevant to the whole human body perception. Nevertheless, this range is amplified in the model 1b in order to reproduce the excitation frequency associated to the rail corrugation calculated in Eq. (4). Hence, the frequency range studied is set between 2 and 125 Hz.

According to [13], the maximum wavelength to be studied in the model 1a is λ max = 60 m, while the minimum is λ min = 1.5 m. Likewise, the maximum wavelength to be studied in the model 1b is λ max = 60 m, while the minimum is λ min = 0.96 m.

With regard to the dimensions of the model and its elements, according to [14], the minimum wavelength must be developed in 6 nodes at least. It must also be taken into account the fact that the model length must be sufficient to allow an accurate study of the maximum wavelength. Thus, the minimum length of the model is set to 60 m, while the maximum size of the mesh elements (d max) is set to 0.48 m (model 1a) and 0.192 m (model 1b). However, in order to represent the sinusoidal loads caused by rail corrugation, a maximum distance of separation between sections is set in the full model of d s = 0.06 m (Fig. 4), coinciding with \(\lambda_{\text{c}} /2\). Furthermore, in order to size elements as regularly as possible, the thickness of the rail pad is increased from 0.007 to 0.05 m maintaining its equivalent stiffness.

Fig. 4
figure 4

Model mesh

Full size image

Although all the elements of the model have their movements coupled, this does not happen in the ballast-sleeper contact. The existing friction between these two materials requires decoupling movements on their interface. Therefore, it is established that the only movement coupled between them is the perpendicular to their planes, leaving the other movements free.

Furthermore, in order to represent the confinement of the ground, the degrees of freedom in the boundaries are restricted. Thus, in the boundaries of the model, the transversal displacements are set to zero.

3.1.2 Load application

Loads applied by the vehicle movement are caused by two main components: one is the quasi-static component (Q E) due to the weight of the vehicle and another is the dynamic component, resulting from imperfections in the wheel–rail contact, geometric changes or stiffness variations.

The quasi-static component has been estimated as Q E = 142.25 kN/axis in Table 2, while the dynamic component has been calculated with experimental equations during the pre-calibration process and with a dynamic model during the calibration process.

In this sense, the Eisenmann’s equation is adopted to consider the vertical dynamic loads as Eq. (6), where Q dv is the dynamic load transmitted to the track; t is the corresponding standard deviation; \(\bar{s}\) is a coefficient that evaluates the quality of the track; and φ is a factor which depends on the train speed obtained in Eq. (7), wherein v corresponds to the speed of the vehicle in km/h.

$$Q_{\text{dv}} = Q_{\text{E}} \left( {1 + t\,\bar{s}\,\varphi } \right).$$
(6)
$$\varphi = 1 + \frac{v - 60}{380}.$$
(7)

Similarly, lateral dynamic overloads are calculated according to [15] as

$$H = H_{\text{C}} + H_{\text{A}} ,$$
(8)

where \(H\) represents the lateral dynamic overloads component, H C represents the acceleration component due to the unbalanced acceleration obtained in Eq. (9), and H A represents the component related to possible defects existing in the wheel–rail contact obtained in Eq. (10).

$$H_{\text{C}} = {\frac{{P\,{v^{2}} }}{Rg} - \frac{P\,h}{S}}.$$
(9)
$$H_{\text{A}} = {\frac{P\,v}{1,000}}.$$
(10)

In these last two equations, P represents the load transmitted by vehicle axle and g is the gravity acceleration.

Then, these data are introduced into the numerical model 1a and the pre-calibration process begins. The dynamic overloads applied to the complete model 1b are developed in Sect. 3.2.

3.1.3 Model pre-calibration

For the first approach, results are studied only for the high rail (which presented no corrugation). The accelerations registered in the rail web are used to calibrate the model (Fig. 5a), while the accelerations registered in the sleeper are used to validate the model (Fig. 5b). The numerical model’s results (blue) show an acceptable adjustment to the ones obtained in the data-gathering campaign (black). Anyway, the differences between both are accepted due to the indicative nature of this process.

Fig. 5
figure 5

Vertical accelerations obtained in the pre-calibration in the rail web (a) and sleeper (b)

Full size image

Once mechanical properties of the track materials are approximated, the equivalent stiffness of the track is calculated by running a static analysis in the numerical model. In this analysis, a unit load Q a = 1 N is applied and the deflection δ is obtained. Subsequently, the equivalent stiffness \(K_{{{\text{eq}}_{n} }}\) is obtained using Eq. (11):

$$K_{{{\text{eq}}_{n} }} = \frac{{Q_{\text{a}} }}{\delta } .$$
(11)

This equation, applied to the pre-calibrated model, provides a vertical equivalent stiffness of \(K_{{{\text{ev}}_{0} }} = 9 5 {\text{ MN/m}}\), while lateral equivalent stiffness is \(K_{{{\text{et}}_{0} }} = 30{\text{ MN/m}}\).

3.2 Multi-body model

Based on the studies carried out by several authors as [16], the non-negligible influence of the vehicle model in the evaluation of the train–track interaction loads allows adopting significant simplifications in vehicle–track model interaction. Although a greater level of detail would provide more accurate results, the great mathematical simplification of the two mass models against complete models allows us to omit certain suspension parameters that are unknown.

The dynamic model has been developed in the VAMPIRE PRO software. The three model inputs are the main features of the vehicle included in Table 2, the geometrical characteristics, and the equivalent stiffness of the track. According to previous works, as those carried out by Jin et al. [5], Ferrara et al. [17], or Uzzal et al. [18], the vehicle is simplified to a single wagon of two bogies and four axles. Only the primary suspension is considered.

For the wheel–rail contact, a full non-linear analysis based on the Kalker’s contact theory has been adopted. Thus, the creepages and creep forces are calculated from the instantaneous axle load and contact area generated at each timestep t. Hence, the wheel profile S1002 has been adopted and, according to [19], a friction coefficient μ = 0.35 is set for the high rail and μ = 0.4 for low rail. Likewise, in order to represent the contamination on the rail surface, the Kalker coefficient is reduced to 0.8. The results of the dynamic overloads obtained in each axis are shown in Fig. 6.

Fig. 6
figure 6

Dynamic overloads obtained in the wheel–rail contact for the 4 axles of the vehicle. a Lateral overloads. b Vertical overloads

Full size image

The graphs in Fig. 6 show the great importance of corrugation on the loads generated in the wheel–rail contact. It can be seen that the dynamic vertical load generated in the low rail presents an oscillatory behavior with higher amplitude than in the high rail (no corrugated).

In all cases, a sinusoidal behavior is observed, whose wavelength matches the wavelength of rail corrugation. However, the forces generated at the wheel–rail contact of both axles of the same bogie show different behaviors. The main reason for this phenomenon lies in the different angle in which each wheelset attacks the track. Thus, different oscillations take place, causing creepages and variations in the contact forces and in the wheel–rail contact point. This leads to different dynamical behavior of the wheelsets. Consequently, the larger the yaw angles is, the greater creep forces are. Nevertheless, in order to evaluate the influence of vehicle suspension parameters on the dynamic behavior of wheel–rail contact forces, the stiffness and damping values of the vehicle suspension have been modified in a sensitive analysis according to Table 4.

Table 4 Sensitivity analysis of the vehicle suspensions
Full size table

Once the sensitivity analysis was carried out, both amplitude and average value of the dynamic overloads were analyzed. As an example, in order to evaluate these differences, the variations of the vertical dynamic overloads generated in the wheel–rail contact of the second axle of the vehicle, when damping coefficient is reduced to 10 kN s/m (Fig. 7b), are compared to the original case (Fig. 7a).

Fig. 7
figure 7

Comparison between the vertical dynamic overloads in axle 2 due to variations of suspension damping coefficient. a Original case. b Case I

Full size image

We can see that a new sinusoidal behavior induces some differences in the amplitude of dynamic overloads. Nevertheless, in all the cases, only minor differences were obtained.

4 Calibration and validation of the FE model

The closest approximation to reality that involves the calculation of dynamic overloads using multi-body models requires dismissing those forces obtained by experimental formulations [Eqs. (6) and (8)]. Hence, a new calibration and validation process is carried out in order to achieve an accurate value of equivalent track stiffness.

In this way, the overloads generated in the wheel–rail contact shown in Fig. 6 are applied in every section of the model 1b. For this purpose, an auxiliary vector (\(\varvec{q}_{i}\)) is created for each wheel representing the sum of the quasi-static component (Q E) and the dynamic component (F n ) obtained according to Eq. (12):

$$\varvec{q}_{i} = Q_{\text{E}} + \varvec{F}_{n} .$$
(12)

Thus, the passage of a complete wagon (4 axles, 8 wheels) is applied to the FE model. Then, assuming linear elastic material behavior, the superposition principle is applied to the three wagons of the train. In order to achieve the ultimate objective of this paper, the iterative calibration process is performed on both rails of the track. In particular, the registers obtained in the low rail are used for the calibration, while the ones obtained in the high rail are used in the validation process. The results of both vertical and lateral dynamics obtained in the numerical model (red) are compared with the results of the data-gathering campaign (black) in Figs. 8 and 9.

Fig. 8
figure 8

Results of the calibration process (a) and validation (b) of rail vertical accelerations

Full size image
Fig. 9
figure 9

Results of the calibration (a) and validation process (b) of rail lateral accelerations

Full size image

In this case, a correct performance of the model is achieved for vertical and lateral accelerations in both rails. Hence, the model is calibrated and validated, and the mechanical properties resulting from the calibration process are shown in Table 5.

Table 5 Main properties of the materials calibrated
Full size table

For these values, the vertical equivalent track stiffness has been reduced a −10.5 % over \(K_{{{\text{ev}}_{0} }}\), while the lateral ones has been also reduced a −6.6 % over \(K_{{{\text{et}}_{0} }}\). The dynamic overloads obtained with the new stiffness are applied to the numerical model and the results are shown in Figs. 10 and 11.

Fig. 10
figure 10

Vertical accelerations obtained for both stiffness and both rail. a Low rail. b High rail

Full size image
Fig. 11
figure 11

Lateral accelerations obtained for both stiffness and both rail. a Low rail. b High rail

Full size image

As expected, due to the track stiffness variation, the vertical and lateral accelerations have been modified. However, this difference is lower than 9 % in all cases. These variations are considered perfectly acceptable, and therefore, the model is finally calibrated and validated.

5 Influence of rail grinding

In order to reproduce the effect of rail grinding, a recalculation in the dynamic model is carried out. For this case, rail corrugation of the inner rail is deleted and the rest of the mechanical properties of the track are kept. Furthermore, as a result of the rail grinding, friction coefficients of both rails are set to μ = 0.35 and Kalker’s coefficient is increased to 1. The new overloads obtained are shown in Fig. 12.

Fig. 12
figure 12

Dynamic overloads calculated in the wheel–rail contact for the 4 axles of the vehicle in a grinded track. a Lateral overloads. b Vertical overloads

Full size image

It can be observed that the amplitude of the dynamic overloads disappears in every case. However, the sign of force is kept in 15 of the 16 cases (Tables 6 and 7). Moreover, although the relative difference of the average value of these forces is higher than 10 % in 4 cases, this relative difference is only remarkable in the lateral forces resulting from the passage of the wheel of the second axle through the high rail (Figs. 13 and 14).

Table 6 Average value of the vertical dynamic overloads
Full size table
Table 7 Average value of the lateral dynamic overloads
Full size table
Fig. 13
figure 13

Average value of the vertical dynamic overloads for the 4 axles of the vehicle

Full size image
Fig. 14
figure 14

Average value of the lateral dynamic overloads for the 4 axles of the vehicle

Full size image

After analyzing the results, it can be concluded that the amplitude of the dynamic overloads obtained is mainly caused by rail corrugation, while the average value is mainly caused by geometric variations in the track. Therefore, in the absence of any other pathology, rail grinding can counteract the amplitude variation of the wheel–rail contact forces.

From the vibrational point of view, due to the variation of the excitation forces generated in the wheel–rail contact, a decrease can be expected in the accelerations registered. This variation is calculated in the FE model according to Eq. (5). Thus, since the mass, damping and stiffness matrices, as well as β Rayleigh’s coefficient remain constant throughout the simulation, any variation in the forces vector should compel a variation in the displacement, speed, and acceleration vectors.

Moreover, since all nodes are coupled, any variation in the excitation forces not only will generate a displacement of the node application, but also will displace the ones located in its vicinity. Thus, the velocity and acceleration component will also be altered, generating a wave which is propagated through the ground in all directions. Hence, the greater the amplitude of the excitation forces is, the higher the vibrational response is.

To estimate this variation, the new calculated overloads are introduced into the calibrated and validated numerical model and the vibrational behavior of the track is shown as Figs. 15 and 16.

Fig. 15
figure 15

Vertical vibration response in the low rail (a) and high rail (b)

Full size image
Fig. 16
figure 16

Lateral vibration response in the low rail (a) and high rail (b)

Full size image

Indeed, a severe reduction of the vibratory response of the track has been reached. Regarding to vertical dynamics, results show a 50 % reduction over the maximum peaks of the low rail. This reduction becomes less noticeable on the high rail, where difference is set to 28 %. Regarding lateral dynamics, the differences have been increased. Here, an 87 % reduction in low rail peaks and 61 % reduction in the high rail are observed.

This demonstrates the importance of rail corrugation in the vibration generation phenomenon. Results have shown that the vibrational response of a corrugated track is much higher than that with no defects; this condition is especially noticeable in lateral dynamics.

6 Conclusions

The proposed method, which consists of a time domain feedback between a multi-body model of the vehicle and a finite element numerical model of the track, is an interesting tool to study railway vibrations in time domain if some of the ground properties are unknown.

With regard to the dynamic overloads obtained by the dynamic model, it has been observed that the amplitude in the low rail is much higher than the one obtained in the high rail. The reason for this is that the corrugation phenomenon has only been depicted on the inner rail, while the response of high rail is the result of the interaction between the two wheels on the same axle. Furthermore, rail corrugation wavelength has been perfectly reproduced in dynamic overloads, obtaining relative maximums and minimums separated each rail corrugation wavelength.

In addition, the dynamic overloads obtained for the four wheels of the same bogie have presented different behaviors as a result of the curved stretch. For this reason, the eight wheels of the wagon have been implemented in the numerical model instead of considering the superposition principle. This has considerably increased calculation times, but an improvement in the accuracy of the results has been obtained.

Moreover, a noticeable reduction of dynamic overloads after rail grinding has been achieved. Although it is observed that the sign and average value of the overload have remained virtually unchanged in 15 of the 16 cases after the grinding process, the amplitude of the forces of wheel–rail interaction has disappeared. Therefore, in the absence of any other pathology, rail grinding will allow counteracting the amplitude variation of the wheel–rail contact forces.

Finally, from a vibrational point of view, it has been proved that the numerical model has adequately responded to variations in the dynamic overload calculated. Once the loads corresponding to a grinded track are introduced, the vibrational response is reduced considerably. Consequently, it can be concluded that the grinding process of a corrugated rail surface can significantly reduce the vibrations generated by passage of a vehicle, especially in lateral dynamics. Nevertheless, the fact of introducing loads as discrete forces at different timesteps has induced a predominant frequency in the numerical model corresponding to the loadsteps application. For further investigations, the authors will work to overcome this handicap.