Frequency Analysis of the Railway Track under Loads Caused by the Hunting Phenomenon Frequency Analysis of the Railway Track under Loads Caused by the Hunting Phenomenon

: Hunting is a potentially dangerous phenomenon related to the lateral oscillation of the wheels that impacts the rails and causes the wear of the infrastructure’s components. Therefore, the analysis and timely detection of hunting can lead to the application of corrective maintenance tasks, reducing damages, and costs and even derailments as a result. In this work, the vibration response of a ﬁnite element model of a rail with hunting-induced loads by a single wheel is analyzed in three directions: longitudinal, lateral, and vertical. The contact patch forces are calculated by means of Simpack ® using the Kalker linear theory and the contact Hertz theory. The system is solved by using the Newmark- β approach. The results of the deﬂection and vibration analysis, following the proposed methodology, show how the different characteristics of the loads impact the rail. Abstract: Hunting is a potentially dangerous phenomenon related to the lateral oscillation of the wheels that impacts the rails and causes the wear of the infrastructure’s components. Therefore, the analysis and timely detection of hunting can lead to the application of corrective maintenance tasks, reducing damages, and costs and even derailments as a result. In this work, the vibration response of a finite element model of a rail with hunting-induced loads by a single wheel is analyzed in three directions: longitudinal, lateral, and vertical. The contact patch forces are calculated by means of Simpack ® using the Kalker linear theory and the contact Hertz theory. The system is solved by using the Newmark- β approach. The results of the deflection and vibration analysis, following the proposed methodology, show how the different characteristics of the loads impact the rail.


Introduction
Hunting is a wheel lateral oscillation that occurs when the wheelset reaches its critical speed. It is well kNown that the main factor causing this phenomenon is the semi-conical shape of the threads, as shown in Figure 1. Nonetheless, the conical shape is needed to keep trains on the tracks and able to turn on curves.

Introduction
Hunting is a wheel lateral oscillation that occurs when the wheelset reaches its critical speed. It is well known that the main factor causing this phenomenon is the semiconical shape of the threads, as shown in Figure 1. Nonetheless, the conical shape is needed to keep trains on the tracks and able to turn on curves. Hunting limits the stable and safe maximum speed of any railway vehicle. Over the past few decades, numerous authors have made many attempts to understand the Hunting limits the stable and safe maximum speed of any railway vehicle. Over the past few decades, numerous authors have made many attempts to understand the mechanism of this phenomenon, aiming to increase the top-speed of trains [1]. One the biggest risks in this violent oscillation is the wear of the railway infrastructure [2]. As the vertical and lateral loads between the wheel and the rail increase, the system becomes unstable which may result in vehicle derailment [3].
Hunting frequencies vary between 2 to 8 Hz [4], depending on the characteristics of the suspension components. In this regard, it is highly helpful to understand how the rail will respond when hunting occurs. During this oscillatory phenomenon, the body of the rails vibrates in a very specific way [5], which can be easily seen with a specialized simulation tool.
Simulators are utilized in the analysis of many important issues in the railway field, such as condition monitoring [6][7][8][9], rail noise [10,11], natural frequencies and vibration modes in bridges [12][13][14][15], and the analysis of track sleeper voids [16,17]. To simulate the vibrations of the rails, many authors use numerical methods, such as the finite element method (FEM). The FEM consists of discretizing the bodies into smaller elements to be analyzed. Such elements are interconnected through nodes with at least one degree of freedom. An example of a railway wheel moving on a rail track at a speed v is shown in Figure 2a. The forces F , F , F arise on the wheel-rail interaction contact. Figure 2b shows a section of the rail modeled with FEM using 12-degree-of-freedom beam elements. The stiffness of the rail and the sleeper are represented by k , k , k , and k , k , k , accordingly. The Figure 2c represents a single beam element between nodes and . In [18], the author obtains a coupling vibration model to investigate the impact when a vehicle runs over the connection between a floating-plate track and a ballasted track. There is a difference in rail stiffness, which makes impact between the rail and the wheel likely to take place. In [19], a vertical track model based on Timoshenko elements was developed. The author implemented improvements that now make it suitable to perform an accurate time domain. As for the frequency domain, the dynamic local system approach is capable of representing the track frequency content by the use of a single FEM element. In [20], the effect of sleeper pitch on rail corrugation on a tangent track is investigated using a numerical method in vehicle hunting. If the hunting of the vehicle running on the tangent track does not occur, rail corrugation, caused by the wear mechanism, does not form easily. In [5], the effect of different models of track flexibility to simulate the running dynamics of railway vehicles on straight and curved tracks is investigated. It is shown that neglecting the effect of track flexibility results in an overestimation of critical speed by more than 10%.
The above literature analysis depicts a range of studies on rail models and their applications. However, none of them analyze the frequency response of a 3D rail model In [18], the author obtains a coupling vibration model to investigate the impact when a vehicle runs over the connection between a floating-plate track and a ballasted track. There is a difference in rail stiffness, which makes impact between the rail and the wheel likely to take place. In [19], a vertical track model based on Timoshenko elements was developed. The author implemented improvements that now make it suitable to perform an accurate time domain. As for the frequency domain, the dynamic local system approach is capable of representing the track frequency content by the use of a single FEM element. In [20], the effect of sleeper pitch on rail corrugation on a tangent track is investigated using a numerical method in vehicle hunting. If the hunting of the vehicle running on the tangent track does not occur, rail corrugation, caused by the wear mechanism, does not form easily. In [5], the effect of different models of track flexibility to simulate the running dynamics of railway vehicles on straight and curved tracks is investigated. It is shown that neglecting the effect of track flexibility results in an overestimation of critical speed by more than 10%.
The above literature analysis depicts a range of studies on rail models and their applications. However, none of them analyze the frequency response of a 3D rail model when the hunting is presented, which is useful in distinguishing the different components of frequency in all directions in the impact zones, in detail. The rail model can be useful to test a track design or to understand how some phenomena, such as hunting and rail corrugations, develop. In this work, as the main contribution, the vibration response of a Mathematics 2022, 10, 0 3 of 17 rail FEM model with hunting-induced loads by the passage of a single wheel is analyzed in the longitudinal, lateral, and vertical directions. The contact patch forces are calculated with Simpack using the Kalker linear theory and the contact Hertz theory. The rail dynamic response is solved using the Newmark-β approach. Some sections of the rail, subjected to high forces, are compared with sections subjected to smaller forces. Additionally, exploiting the potential of new and powerful signal processing techniques, the high frequencies of the rail response are exposed using the Empirical mode decomposition (EMD), and analyzed with the continuous wavelet transform using fast Fourier transform (CWTFT) to show the frequencies at selected sections of the track. The response of the rail model is validated with real measurements. The calculations involving the contact friction and tangential forces are solved by the Simpack software. It is demonstrated that the combination of the multibody simulation software and a simple track model is realistic and saves computational resources.

Mathematical Model
In this section, the development of the proposed rail model is presented.

Rail Model
The finite element model of the rail is built on beam elements connected in a series and is supported by springs that simulate the sleeper stiffness, as illustrated in Figure 3. The hunting loads, pre-calculated with the commercial software Simpack, vary in time and position. As the wheel moves forward, it interacts with each of the nodes through the wheel-rail forces in three directions. The model contains 1020 beam elements connected through nodes n 1 , n 2 , n 3 . . . n 1021 . The rail is mounted on 18 sleepers (sleeper D1, sleeper D2, . . . sleeper D18) with 60 elements between them. The number of elements between sleepers was chosen to be 60 to expand the contact patch size of the wheel and the rail (11 mm) to the beam element size [21]. A two-dimensional FEM representation of the first rail span is illustrated in Figure 3 on the x-z plane. The interaction forces are represented by the red arrows. The nodes in this rail span are n 1 , n 2 , n 3 . . . n 61 , and k SZ and k SθY are the vertical and rotational stiffness of the sleepers, accordingly.  Each node features six degrees of freedom, u , u , u , u , u and u , meant to capture all the forces and characteristics of the hunting behavior of the wheelset. Each sleeper is positioned on one single node. Thus, the stiffness of those nodes increases in all directions and rotations, i.e., the total stiffness of the nodes of the sleepers is k k , in the vertical direction, and k k , in the rotation around y.
The 12 degrees of freedom of each element are represented by the vector (see Figure   Figure 3. The FEM model of the rail. The wheel forces are applied to one-by-one nodes at each step of the simulation. View 1 illustrates in detail how the interaction forces are applied to the nodes.
Each node features six degrees of freedom, u x , u y , u z , u θx , u θy and u θz , meant to capture all the forces and characteristics of the hunting behavior of the wheelset. Each sleeper is positioned on one single node. Thus, the stiffness of those nodes increases in all directions and rotations, i.e., the total stiffness of the nodes of the sleepers is k Z + k SZ , in the vertical direction, and k θY + k SθY , in the rotation around y.
The 12 degrees of freedom of each element are represented by the vector (see Figure 4): {u i } = u xi u yi u zi u θxi u θyi u θzi u xi+1 u yi+1 u zi+1 u θxi+1 u θyi+1 u θzi+1 (1) Figure 3. The FEM model of the rail. The wheel forces are applied to one-by-one nodes at each step of the simulation. View 1 illustrates in detail how the interaction forces are applied to the nodes.
Each node features six degrees of freedom, u , u , u , u , u and u , meant to capture all the forces and characteristics of the hunting behavior of the wheelset. Each sleeper is positioned on one single node. Thus, the stiffness of those nodes increases in all directions and rotations, i.e., the total stiffness of the nodes of the sleepers is k k , in the vertical direction, and k k , in the rotation around y.
The 12 degrees of freedom of each element are represented by the vector (see Figure  4): The dynamical model equation that rules this model is: Figure 4. A single 3D beam element featuring the 12 degrees of freedom of the rail model. The forces F xi , F yi , F zi and moments M xi , M yi , M zi are the wheel forces applied to the node n i . At the next step, the forces F xi+1 , F yi+1 , and F zi+1 , and the moments M xi+1 , M yi+1 , and M zi+1 are applied to the node n i+1 , and so on.
The dynamical model equation that rules this model is: ..
The vectors {u t } and ..   For each beam element, there is a 12 × 12 stiffness matrix as follows [22]: where k 1 , k 2 , k 3 , and k 4 are the stiffness element matrices of the Euler-Bernoulli beam. The stiffness depends on the cross-sectional area A, Young's modulus E, and the moments of inertia I y and I t , which are characteristics of the rail. The mass matrices are of the same dimension as the stiffness matrix, having the form:

Wheel-Rail Interaction Forces
The forces applied to the node n i of the rail model are F i = Fx i Fy i Fz i Mθ xi Mθ yi Mθ zi (5) where i = 1, 2, 3, . . . 1021 is the index of the FEM nodes. The elements of vector {F} change over time, as depicted in Table 1. When the wheel is at the first position, t = 0 s, the wheel is applying the forces at node n 1 . At t = ∆t s, the wheel is over node n 2 . At t = 2∆t s, the wheel is over n 2 , and so on.
In this work, the wheel-rail interaction forces {F} were pre-calculated using Simpack. These forces are applied to the rail during the simulation. Many authors [23][24][25][26] use Simpack to simulate rail vehicle dynamics and wheel-rail interaction forces. Simpack computes the forces that are developed on the contact patch using an elastic approach with spring/damper elements [27]. The contact is calculated using the Hertz method, see Figure 5. The semi-axes a and b of the equivalent contact ellipse are calculated using this theory.
Mathematics 2022, 10, x FOR PEER REVIEW 6 of 18 Figure 5. According to the Hertz contact theory, the contact patch between the rail and the wheel is an ellipse with semi-axis and . View 2 illustrates in detail the forces on the contact patch. The semi-axes of the contact patch depend on the magnitude of the force P. The forces applied by the wheel F , F , and F , and the longitudinal velocity cause the appearance of tangential forces , and a normal force at the contact patch. The reference frame of the contact patch xy'z' change during the passage because of the conical shape and the oscillation of the wheel, and so it is slightly different from the orthogonal reference frame xyz (see Figure 6).
If the pressure distribution x, y is known, then the total normal load can be found by integration: The maximum pressure depends on the combination of Young's modulus * , the Poisson ratio , and the applied force over the bodies [21]. The transformation from the xy'z' frame to the orthogonal xyz frame is illustrated in Figure 6. The normal force is in the z' direction. The tangential force is in the y direction. The transformation matrix is Figure 5. According to the Hertz contact theory, the contact patch between the rail and the wheel is an ellipse with semi-axis a and b. View 2 illustrates in detail the forces on the contact patch. The semi-axes of the contact patch depend on the magnitude of the force P. The forces applied by the wheel F x , F y , and F z , and the longitudinal velocity v cause the appearance of tangential forces T x , T y and a normal force N at the contact patch. The reference frame of the contact patch xy z change during the passage because of the conical shape and the oscillation of the wheel, and so it is slightly different from the orthogonal reference frame xyz (see Figure 6).
If the pressure distribution p(x, y) is kNown, then the total normal load can be found by integration: The maximum pressure p 0 depends on the combination of Young's modulus E * , the Poisson ratio ν, and the applied force P over the bodies [21]. The transformation from the xy z frame to the orthogonal xyz frame is illustrated in Figure 6. The normal force N is in the z direction. The tangential force T y is in the y direction. The transformation matrix is given by: during the passage because of the conical shape and the oscillation of the wheel, and so it is slightly different from the orthogonal reference frame xyz (see Figure 6).
If the pressure distribution x, y is known, then the total normal load can be found by integration: The maximum pressure depends on the combination of Young's modulus * , the Poisson ratio , and the applied force over the bodies [21]. The transformation from the xy'z' frame to the orthogonal xyz frame is illustrated in Figure 6. The normal force is in the z' direction. The tangential force is in the y direction. The transformation matrix is given by: According to Kalker's linear theory, the tangential forces developed in the wheel-rail contact area are defined as [28]: According to Kalker's linear theory, the tangential forces developed in the wheel-rail contact area are defined as [28]: in which f 11 , f 12 , and f 33 are the creep coefficients, in function of the modulus of rigidity and the semi-axis of the contact ellipse. These coefficients were studied by Kalker and have been used for decades in the railway field to calculate the tangential forces on the contact patch. The terms ξ x , ξ y , and ξ sp are the longitudinal, lateral, and spin creepages, respectively. Creepages, the main factor that triggers the hunting, refer to the difference in tangential speeds between the wheel and the rail on the contact patch. Therefore, the tangential forces on the contact area rely on the normal force, the friction coefficient, as well as the contact patch geometry and the pressure distribution throughout the patch. The model built in Simpack involves a single axis running at the critical speed of 50 m/s in a straight line. The wheel type S1002 and the rail profile 54E1 are the default components in Simpack. To induce hunting oscillation, the axis is moved initially to y = 0.0055 m. The axle has no primary suspension and has a total vertical load of −130 kN so the oscillation never stops, and all the characteristics of the hunting behavior can be easily obtained. The maximum pressure p 0 between wheel and rail computed by Simpack is displayed in Figure 7. in which , , and are the creep coefficients, in function of the modulus of rigidity and the semi-axis of the contact ellipse. These coefficients were studied by Kalker and have been used for decades in the railway field to calculate the tangential forces on the contact patch. The terms , , and are the longitudinal, lateral, and spin creepages, respectively. Creepages, the main factor that triggers the hunting, refer to the difference in tangential speeds between the wheel and the rail on the contact patch. Therefore, the tangential forces on the contact area rely on the normal force, the friction coefficient, as well as the contact patch geometry and the pressure distribution throughout the patch.
The model built in Simpack involves a single axis running at the critical speed of 50 m/s in a straight line. The wheel type S1002 and the rail profile 54E1 are the default components in Simpack. To induce hunting oscillation, the axis is moved initially to = 0.0055 m. The axle has no primary suspension and has a total vertical load of −130 kN so the oscillation never stops, and all the characteristics of the hunting behavior can be easily obtained. The maximum pressure between wheel and rail computed by Simpack is displayed in Figure 7. The results of the Simpack simulation are the forces , and . The , forces are transformed into the orthogonal forces and , respectively, which are taken as the input forces for the rail model. According to the wheel position on the track, the overall forces are rotated and aligned from the reference frame xy'z' to the reference system xyz, see The results of the Simpack simulation are the forces N, T y and T x . The N, T y forces are transformed into the orthogonal forces Y and Q, respectively, which are taken as the input forces for the rail model. According to the wheel position on the track, the overall forces are rotated and aligned from the reference frame xy z to the reference system xyz, see Figure 5; thus: where [R] is the rotation matrix (7). The tangential force T x does not need a transformation and can be taken as it is. Figure 8 shows the forces resulting from the Simpack model. The time interval is 0.238 s, which corresponds to exactly 17 sleepers at 50 m/s and covers approximately two hunting cycles. The magnitude and the frequencies of the forces do change during this interval. The hunting main-oscillation frequency is 7.8 Hz, yet higher frequencies emerge due to the friction of the interaction surfaces. Note that the Y and T x forces oscillate around the equilibrium point, while the Q force is the highest and is almost flat. The vertical force Q is slightly different from the normal force from Figure 7a; N is in the z direction, whereas Q is in the z direction. The low-amplitude force occurs around t = 0.045 s and the peaks of high amplitude at t = 0.126 s and t = 0.184 s. The high peak around t = 0.126 s is caused by the lateral movement of the wheel during the hunting. These points are of interest for the analysis of the rail response. The forces of Figure 8 are sampled as they contain the same number of elements as the total number of nodes in the finite element model. In this way, each of these forces takes the vector form: The forces of Figure 8 are sampled as they contain the same number of elements as the total number of nodes in the finite element model. In this way, each of these forces takes the vector form: Mathematics 2022, 10, 0 F x (n) stands for the samples of T x (nt) at uniform intervals of ∆t seconds, in which n is the total number of nodes in the rail. The same procedure is applied to F y (n) and F z (n). Therefore, as the wheel moves forward, each of the elements of these vectors excites the nodes one by one until the simulation finishes.
Another output of Simpack is the location of the contact point U, which is important so as to apply moments to the degrees u θx , u θy , and u θz of the nodes in the rail. Consequently, the moments applied to the nodes are: where U y and U z are the coordinates of the contact point in the y and z axis, respectively. These coordinates are measured with ground on the orthogonal reference frame of the Simpack model.

Method of Solution
To solve this system, the Newmark-β method [28] was used. The vectors . u t and .. u t at time t + ∆t are calculated based upon the following expressions: .. .
Within this method, an effective mass matrix [m] and an effective force vector F t+∆t are defined as follows: [m] = 1 Thus, it is possible to calculate the displacement of the nodes {u t+∆t }:

Methodology
In general, once the FEM model is obtained and its dynamics are computed by the solver, the frequency information from the time response of the rail induced by the wheel forces can be obtained and analyzed. The flowchart is shown in Figure 9, whereas the stages for Dynamic Models and Solver have been described in previous sections of this paper.

Methodology
In general, once the FEM model is obtained and its dynamics are computed by the solver, the frequency information from the time response of the rail induced by the wheel forces can be obtained and analyzed. The flowchart is shown in Figure 9, whereas the stages for Dynamic Models and Solver have been described in previous sections of this paper. Figure 9. Flowchart for the frequency analysis of the rail response.
The first stage of the flowchart describes the wheel and rail interaction models. In these blocks, the general system properties and forces of the wheel and rail, listed in Table  2, are entered. The wheel forces are computed from the dynamic model created in Simpack's software. The rail model is described by using FEM. The parameters used in this model are summarized in Table 2 [18,29]. Table 2. Principal parameters of the rail model (data from [18,29] Figure 9. Flowchart for the frequency analysis of the rail response. The first stage of the flowchart describes the wheel and rail interaction models. In these blocks, the general system properties and forces of the wheel and rail, listed in Table 2, are entered. The wheel forces are computed from the dynamic model created in Simpack's software. The rail model is described by using FEM. The parameters used in this model are summarized in Table 2 [18,29]. Table 2. Principal parameters of the rail model (data from [18,29] The rail model and forces are arranged in a matrix form, Equation (2), before entering the solver. The displacements of the nodes resulting from the wheel forces interaction are obtained following the Newmark-β method discussed in Section 2.3, Equations (18)- (22). The integration time step is 10 −4 s. Now, for the frequency analysis, the vibrations of the nodes are analyzed with the EMD method. In this part, the high-frequency vibration components are preferred, so that small changes in the wheel-rail contact area can be exposed in detail.
The high-frequency changes are analyzed with the Matlab's CWTFT function. The output of this analysis is a spectrogram showing the different frequencies that appear at every point of the rail during the simulation. In general, these techniques are chosen because of their advantages in applications on nonlinear mechanics [30] and computational time [31], respectively. The results of rail frequencies are explained in the next section.

Rail Deflections
By solving the system of Equation (2), the linear and angular displacements of the rail nodes over time and position are obtained. The largest deflections take place in the vertical direction. In Figure 10, a series of curves showing the deflections of the rail in the vertical direction when the wheel passes from sleeper D4 to D7 are plotted. Since the vertical forces seem to be almost steady (no oscillation) and in only one direction, the behavior of the rail deflections in this direction is rather simple. This pattern repeats over and over until the end of the track. The sleepers are represented by red dots. When the wheel is in close proximity to the sleepers, the deflection is small, see Figure 10b,f. In the midspan, the rail stiffness is softer and, consequently, the deflection is bigger, see Figure 10c-e. vertical direction. In Figure 10, a series of curves showing the deflections of the rail in the vertical direction when the wheel passes from sleeper D4 to D7 are plotted. Since the vertical forces seem to be almost steady (no oscillation) and in only one direction, the behavior of the rail deflections in this direction is rather simple. This pattern repeats over and over until the end of the track. The sleepers are represented by red dots. When the whee is in close proximity to the sleepers, the deflection is small, see Figure 10b,f. In the midspan, the rail stiffness is softer and, consequently, the deflection is bigger, see Figure 10ce. In Figure 11, the lateral deflections of the rail during the wheel passing over sleepers D4 to D7 are displayed. The pattern is different in comparison with Figure 10, owing to the fact that the force oscillates around the equilibrium point, as shown in Figure 8b However, the amplitude of these lateral deflections is smaller and oscillates from one side In Figure 11, the lateral deflections of the rail during the wheel passing over sleepers D4 to D7 are displayed. The pattern is different in comparison with Figure 10, owing to the fact that the force Y oscillates around the equilibrium point, as shown in Figure 8b. However, the amplitude of these lateral deflections is smaller and oscillates from one side to the other. The pattern is different in every step of the track, depending on the amplitude and direction of the lateral force Y.  Additionally, a convergence test of the FEM model is shown in Figure 12, varying the number of elements from 3 to 100 (x-axis) between sleepers. The sleeper deflection was selected in this test because it is the most stable displacement. Results presented in Figure  12 show that 40 elements in the x-axis are enough to model the response of the rail since a relatively stable convergence is obtained; however, an approach of 60 elements in the xaxis is used in order to reach the size of the contact patch as mentioned in Section 2.1, indicating that the high frequencies of the hunting forces have a roll on the rail response. Finally, it is observed that a larger number of elements in the x-axis does not provide more information about the rail deflection. Additionally, a convergence test of the FEM model is shown in Figure 12, varying the number of elements from 3 to 100 (x-axis) between sleepers. The sleeper deflection was selected in this test because it is the most stable displacement. Results presented in Figure 12 show that 40 elements in the x-axis are enough to model the response of the rail since a relatively stable convergence is obtained; however, an approach of 60 elements in the x-axis is used in order to reach the size of the contact patch as mentioned in Section 2.1, indicating that the high frequencies of the hunting forces have a roll on the rail response. Finally, it is observed that a larger number of elements in the x-axis does not provide more information about the rail deflection. selected in this test because it is the most stable displacement. Results presented in Figure  12 show that 40 elements in the x-axis are enough to model the response of the rail since a relatively stable convergence is obtained; however, an approach of 60 elements in the xaxis is used in order to reach the size of the contact patch as mentioned in Section 2.1, indicating that the high frequencies of the hunting forces have a roll on the rail response. Finally, it is observed that a larger number of elements in the x-axis does not provide more information about the rail deflection.

EMD Decomposition of the Time Response of the Nodes
The time response of the nodes is a very small-amplitude vibration and has a particular shape, similar to the rail deflection, see Figure 13b. The largest amplitude of every node occurs when the wheel passes on (in the figure, around sample 170) and is limited by the sleepers. Before and after the node's midspan, the vibration is minimum and only a single large trough appears on the entire signal, consequently. By observing Figure 13b, it can be seen that the vibration starts with an up-slope. This effect is caused by the wheel exerting the vertical forces on the previous midspan, causing the node to move upward before entering the node's midspan. Once the wheel enters the node's midspan, the slope goes down until the wheel is positioned right over the node.

EMD Decomposition of the Time Response of the Nodes
The time response of the nodes is a very small-amplitude vibration and has a particular shape, similar to the rail deflection, see Figure 13b. The largest amplitude of every node occurs when the wheel passes on (in the figure, around sample 170) and is limited by the sleepers. Before and after the node's midspan, the vibration is minimum and only a single large trough appears on the entire signal, consequently. By observing Figure 13b, it can be seen that the vibration starts with an up-slope. This effect is caused by the wheel exerting the vertical forces on the previous midspan, causing the node to move upward before entering the node's midspan. Once the wheel enters the node's midspan, the slope goes down until the wheel is positioned right over the node. As a comparison, a vibration measurement captured by a strain gauge installed in the midspan of a rail at a 1594-Hz sampling rate is shown in Figure 13a. In this case, the strain gauge is used to measure vibrations by means of the surface deformations of the rail. The vibration measurement was captured in an experiment using a 1/10-scale railcar running on an inclined surface. More details of this experiment can be seen in [32].
Both signals have similar characteristics, so they are compared to prove the reliability of the model. The experimental signal has four peaks due to the four wheels of the railcar, whereas the time response of the nodes does not because of the single wheel vertical force. The application of the EMD demonstrates that these signals have similar features, despite the experimental signal peaks. The EMD algorithm obtains the oscillating components of a signal, called intrinsic mode functions (IMF), by finding the local extrema in the test data. Every IMF contains signal information that is often difficult to study if the full wave is taken. One characteristic of this method is that the summation of the IMFs returns to the original signal. Some of the IMF can be discarded to obtain a selected-feature signal. The EMD method is suitable for nonlinear and nonstationary processes [33] since the decomposition is based on the local characteristic of the data. The EMD was applied to the experimental and simulation response signals, illustrated in Figure 13, and their IMFs are shown in Figure 14. As a comparison, a vibration measurement captured by a strain gauge installed in the midspan of a rail at a 1594-Hz sampling rate is shown in Figure 13a. In this case, the strain gauge is used to measure vibrations by means of the surface deformations of the rail. The vibration measurement was captured in an experiment using a 1/10-scale railcar running on an inclined surface. More details of this experiment can be seen in [32].
Both signals have similar characteristics, so they are compared to prove the reliability of the model. The experimental signal has four peaks due to the four wheels of the railcar, whereas the time response of the nodes does not because of the single wheel vertical force. The application of the EMD demonstrates that these signals have similar features, despite the experimental signal peaks.
The EMD algorithm obtains the oscillating components of a signal, called intrinsic mode functions (IMF), by finding the local extrema in the test data. Every IMF contains signal information that is often difficult to study if the full wave is taken. One characteristic of this method is that the summation of the IMFs returns to the original signal. Some of the IMF can be discarded to obtain a selected-feature signal. The EMD method is suitable for nonlinear and nonstationary processes [33] since the decomposition is based on the local characteristic of the data. The EMD was applied to the experimental and simulation response signals, illustrated in Figure 13, and their IMFs are shown in Figure 14. The summation of the experimental IMFs is the wave of Figure 15a, whereas the summation of the simulation IMFs is the wave of Figure 15b. It should be noted that the 1, 5, 6, and 7 IMFs of the experimental signal were discarded in order to remove the peaks. A normalized cross-correlation analysis, by means of the xcorr Matlab's function, indicates a 3% error. The similarity of the signals shows that the vibrations of the proposed rail model are similar to real measurements and therefore the model is reliable and suitable for the simulation of the rail response under hunting loads. The summation of the experimental IMFs is the wave of Figure 15a, whereas the summation of the simulation IMFs is the wave of Figure 15b. It should be noted that the 1, 5, 6, and 7 IMFs of the experimental signal were discarded in order to remove the peaks. A normalized cross-correlation analysis, by means of the xcorr Matlab's function, indicates a 3% error. The summation of the experimental IMFs is the wave of Figure 15a, whereas the summation of the simulation IMFs is the wave of Figure 15b. It should be noted that the 1, 5, 6, and 7 IMFs of the experimental signal were discarded in order to remove the peaks. A normalized cross-correlation analysis, by means of the xcorr Matlab's function, indicates a 3% error. The similarity of the signals shows that the vibrations of the proposed rail model are similar to real measurements and therefore the model is reliable and suitable for the simulation of the rail response under hunting loads. The similarity of the signals shows that the vibrations of the proposed rail model are similar to real measurements and therefore the model is reliable and suitable for the simulation of the rail response under hunting loads.

Frequency Domain Analysis during Hunting
The frequency analysis is carried out on the sections of the rail where the highamplitude forces appear. By observing the behavior of the wheel forces in Figure 8, the highlighted spots are the peaks of high amplitude and frequency, around t = 0.126 s and t = 0.184 s. These points are compared with points of small forces, such as the ones at t = 0.045 s to enhance the differences in the rail response. The EMD method is applied to the time domain response of the FEM nodes, prior to making the frequency analysis, to discard the IMFs of lower frequencies caused by the rigid body dynamics. Then, the CWTFT function of Matlab, by using the Morlet wavelet, is applied to the first two IMFs to visualize the different high frequencies in that particular section of the track. The first two IMFs were chosen for the analysis in the longitudinal and vertical directions, whereas only the first one was chosen for the vertical direction. That is because there are fewer high frequencies in the vertical direction. Figure 16 displays the time response and the frequency response of the nodes n 214 , n 542 , and n 790 , which correspond to the forces at times t = 0.045 s, t = 0.126 s, and t = 0.184 s, accordingly, in the longitudinal direction. The left charts are the original time response of the rail at the mentioned times. The charts in the middle are the IMF signals as a result of the applied EMD of those original responses. The right charts are the CWTFT spectrograms of the EMD signals presenting frequencies ranging from 0 to 820 Hz.
The frequency analysis is carried out on the sections of the rail where the high-amplitude forces appear. By observing the behavior of the wheel forces in Figure 8, the highlighted spots are the peaks of high amplitude and frequency, around = 0.126 s and = 0.184 s. These points are compared with points of small forces, such as the ones at = 0.045 s to enhance the differences in the rail response. The EMD method is applied to the time domain response of the FEM nodes, prior to making the frequency analysis, to discard the IMFs of lower frequencies caused by the rigid body dynamics. Then, the CWTFT function of Matlab, by using the Morlet wavelet, is applied to the first two IMFs to visualize the different high frequencies in that particular section of the track. The first two IMFs were chosen for the analysis in the longitudinal and vertical directions, whereas only the first one was chosen for the vertical direction. That is because there are fewer high frequencies in the vertical direction. Figure 16 displays the time response and the frequency response of the nodes , , and , which correspond to the forces at times = 0.045 s, = 0.126 s, and = 0.184 s, accordingly, in the longitudinal direction. The left charts are the original time response of the rail at the mentioned times. The charts in the middle are the IMF signals as a result of the applied EMD of those original responses. The right charts are the CWTFT spectrograms of the EMD signals presenting frequencies ranging from 0 to 820 Hz.
The node , subjected to small forces during = 0.045 s, does not show any interesting behavior as expected, see Figure 16a. However, subjected to bigger forces during = 0.126 s, the time response of the node has a different shape and the CWTFT spectrogram shows higher frequencies only at = 0.126 s, when the wheel passes, see Figure  16b. At = 0.184 s, the node is subjected to smaller forces than node and the spectrogram displays a combination of many frequencies, see Figure 16c. Most of the frequencies in this node are low (around 130 Hz), but at the instant when the wheel passes, a dominant frequency of 250 Hz is presented. The same analysis is applied in the lateral direction at the same locations. The node , to which the force exerts small forces during = 0.045 s, does not show a defined point of impact, thus the frequencies in the spectrogram are not defined, see Figure 17a. The node n 214 , subjected to small forces during t = 0.045 s, does not show any interesting behavior as expected, see Figure 16a. However, subjected to bigger forces during t = 0.126 s, the time response of the node n 542 has a different shape and the CWTFT spectrogram shows higher frequencies only at t = 0.126 s, when the wheel passes, see Figure 16b. At t = 0.184 s, the node n 790 is subjected to smaller forces than node n 542 and the spectrogram displays a combination of many frequencies, see Figure 16c. Most of the frequencies in this node are low (around 130 Hz), but at the instant when the wheel passes, a dominant frequency of 250 Hz is presented.
The same analysis is applied in the lateral direction at the same locations. The node n 214 , to which the force exerts small forces during t = 0.045 s, does not show a defined point of impact, thus the frequencies in the spectrogram are not defined, see Figure 17a. However, the time response of the node n 542 has a more defined shape, due to higher forces. Only a small peak appears at the instant t = 0.126 s and, after that, the node n 542 is pulled to the other side. The CWTFT spectrogram shows the moment of impact at t = 0.128 s, when the high frequencies of the rail stop, and a component of 200 Hz appears.
athematics 2022, 10, x FOR PEER REVIEW 15 of 18 However, the time response of the node has a more defined shape, due to higher forces. Only a small peak appears at the instant = 0.126 s and, after that, the node is pulled to the other side. The CWTFT spectrogram shows the moment of impact at = 0.128 s, when the high frequencies of the rail stop, and a component of 200 Hz appears. At = 0.184 s, the node is subjected to smaller forces than node and the spectrogram displays a combination of many frequencies, see Figure 17c. The dominant frequency is 290 Hz, in this case. The frequency of 120 Hz is displayed due to the wheel arriving at the sleeper. It is important to see that, unlike longitudinal direction, in the lateral direction the nodes and take a moment to react, but this is caused primarily by a larger duration of the peaks in the longitudinal direction, see the differences in Figure  8a,b. High-amplitude forces reduce the appearance of high frequencies; therefore, in order to see the high frequencies in this direction, only the first IMF was selected. Figure 18a displays the behavior of the rail at = 0.045 s, when the force F z is almost constant and delivers low frequencies. The rail behavior is very different compared to the results in Figures 16a and 17a. This is owing to the force F z being very high in comparison to F y and F x , and not having major oscillations. Around = 0.126 s, see Figure 18b, the changes in the frequency of F z cause the rail to have higher frequencies caused by the extreme force at that point. At time = 0.184 s a 200 Hz component is also presented, see Figure 18c. At t = 0.184 s, the node n 790 is subjected to smaller forces than node n 542 and the spectrogram displays a combination of many frequencies, see Figure 17c. The dominant frequency is 290 Hz, in this case. The frequency of 120 Hz is displayed due to the wheel arriving at the sleeper. It is important to see that, unlike longitudinal direction, in the lateral direction the nodes n 542 and n 790 take a moment to react, but this is caused primarily by a larger duration of the peaks in the longitudinal direction, see the differences in Figure 8a,b.
The behavior of the rail in the vertical direction at the locations, t = 0.045 s, t = 0.126 s, and t = 0.184 s, are quite different from the longitudinal and lateral directions. Highamplitude forces reduce the appearance of high frequencies; therefore, in order to see the high frequencies in this direction, only the first IMF was selected. Figure 18a displays the behavior of the rail at t = 0.045 s, when the force F z is almost constant and delivers low frequencies. The rail behavior is very different compared to the results in Figures 16a and 17a. This is owing to the force F z being very high in comparison to F y and F x , and not having major oscillations. Around t = 0.126 s, see Figure 18b, the changes in the frequency of F z cause the rail to have higher frequencies caused by the extreme force at that point. At time t = 0.184 s a 200 Hz component is also presented, see Figure 18c. The amplitude of the beam deflection determines the magnitude of the force, but a frequency analysis can reveal further features in the points of higher impact. The fre quency analysis reveals the different frequencies that appear at the selected points. When the force is high with low frequency, the rail frequency is also low. When the force is high with high frequency, the rail vibration is also high with high amplitude. By observing the CWTFT charts, the wheel impacts are pointed out by abrupt changes in the rail frequency When there are suddenly high interaction forces, a low component appears, indicating the instant when the impact occurs.
Higher frequencies than 800 Hz are due to the way in which the forces are applied during the solution of the system. Thus, these frequencies can be ignored. Despite the high noise, the method is powerful enough to find important features of the wheel-rail interac tion in the rail behavior.

Conclusions
In this work, the vibration response of a rail FEM model with hunting-induced loads by a single wheel is analyzed in longitudinal, lateral, and vertical directions. The resulting hunting forces from a multibody software, based on the Kalker linear theory and the con tact Hertz theory, are applied to the FEM model to induce vibrations to the rail. The high frequencies of the rail during the run are isolated using the EMD method to see the effec of small changes in rail behavior. Then, the filtered signal is analyzed with the CWTFT to reveal where the different frequencies appear in the impact zone. The reliability of the model was proved by comparing the time response of a single node with a signal captured with an experimental railcar.
The analysis of the frequency response was made upon three points of interest in the three directions to look closely at what the response is when the external forces increase The amplitude of the beam deflection determines the magnitude of the force, but a frequency analysis can reveal further features in the points of higher impact. The frequency analysis reveals the different frequencies that appear at the selected points. When the force is high with low frequency, the rail frequency is also low. When the force is high with high frequency, the rail vibration is also high with high amplitude. By observing the CWTFT charts, the wheel impacts are pointed out by abrupt changes in the rail frequency. When there are suddenly high interaction forces, a low component appears, indicating the instant when the impact occurs.
Higher frequencies than 800 Hz are due to the way in which the forces are applied during the solution of the system. Thus, these frequencies can be ignored. Despite the high noise, the method is powerful enough to find important features of the wheel-rail interaction in the rail behavior.

Conclusions
In this work, the vibration response of a rail FEM model with hunting-induced loads by a single wheel is analyzed in longitudinal, lateral, and vertical directions. The resulting hunting forces from a multibody software, based on the Kalker linear theory and the contact Hertz theory, are applied to the FEM model to induce vibrations to the rail. The high frequencies of the rail during the run are isolated using the EMD method to see the effect of small changes in rail behavior. Then, the filtered signal is analyzed with the CWTFT to reveal where the different frequencies appear in the impact zone. The reliability of the model was proved by comparing the time response of a single node with a signal captured with an experimental railcar.
The analysis of the frequency response was made upon three points of interest in the three directions to look closely at what the response is when the external forces increase.
During hunting, the rail oscillates more in the lateral and longitudinal directions. However in the vertical direction, the oscillation is small and the amplitude is high. The instant of higher impact can be easily seen in the charts as low frequencies around 200-300 Hz. Yet, some other frequencies were found in the sections where the loads were bigger and oscillating. Throughout the analysis hereby submitted, it is feasible to spot the moment at which the loads increase since the frequencies change from high to low when there is a sudden increase in the interaction forces.
The hunting frequency did not appear in the analysis because the study only considers the interaction in the zone of bigger deflections. This zone covers approximately 1 m, taking into account the speed of the wheel. At that speed, one cycle of the hunting frequency (7.8 Hz) corresponds to 6.4 m. So, further improvement is needed to see this phenomenon from every point of the rail, i.e., including mechanical physics, such as wave propagation.
Further research related to the actual work should investigate how other mechanical issues of the rail, e.g., fatigue and strain, are affected by the wheel passage using this rail model. Another opportunity of research is to study the rail response considering the ballast damping, which was neglected during the actual simulation.

Data Availability Statement:
The data presented in this study are not publicly available due to privacy issues.