Parametric Study of the Influence of Nonlinear Elastic Characteristics of Rail Pads on Wheel–Rail Vibrations

The rail pad is the elastic element between the rail and the sleeper that has the role of absorbing the mechanical stresses from the rail and reducing the vibrations and shocks generated by wheel–rail interactions. In this paper, the problem of the influence of the variability of the nonlinear load-deformation characteristic of rail pads (resulting from the manufacturing process) on wheel–rail vibrations is investigated. The limit load-deformation characteristics of a manufactured rail pad and the medium load-deformation characteristic resulting as the arithmetic mean of the two are considered. The nonlinear load-deformation characteristic of the ballast is also considered. All these characteristics are approximated with the help of the bilinear function and are implemented in a track model consisting of an infinite Euler-Bernoulli beam placed on a two-elastic layer continuous foundation with inertial insertion, resulting in a model with an inhomogeneous foundation. The parameters of the inhomogeneous foundation are established from the equilibrium condition under a static load. Wheel–rail vibrations are studied in terms of the contact force and the acceleration of the rail and wheel. The influence of the variability of the elastic characteristics of the rail pad manifests itself in the field of medium frequencies, which amplify or attenuate the vibration levels in certain bands of one-third of an octave.


Introduction
The rail pad is a piece of elastic material comprising part of the fastening system that includes other components that provide a structural link between rail and sleeper ( Figure 1 based on [1]). Depending on the type of fastening system, namely, a direct (Figure 1a) or indirect system (Figure 1b), the rail pad is inserted between the rail foot and the sleeper or includes a steel baseplate between the rail pad and sleeper which is fixed to the sleeper using other fasteners than those used to fix the rail by the steel baseplate.
The main roles of the rail pad are to elastically absorb the force from the rail by transmitting it to the sleeper and to damp the vibrations and shocks caused by traffic. To this end, rail pads are typically fabricated from resilient materials like rubber, ethylene propylene diene monomer (EPDM), thermoplastic polyester elastomer (TPE), ethylenevinyl acetate (EVA) or high-density polyethylene (HDPE) [2,3]. Some of these rail pads are presented in Figure 2.
The load-deformation of the rail pad has a specific shape dictated by the way it must behave on the track under the action of the force coming from the rail. The rail pad stiffness must be small at the beginning so that the deformation under the action of the clamping force is large enough to ensure intimate contact with the rail foot, regardless of the vertical movements of the rail. When a train wheel is above the rail pad, its stiffness must be high to prevent large movements that could lead to a weakening of the grip. Such an  It follows from the above that rail pads have nonlinear, quasi-static elastic characteristics, as described in several works [3,5,6]. This nonlinear characteristic of the rail pad, together with the non-linear elastic characteristic of the ballast and subgrade [5,7,8], determine the specific shape of the load-deformation curve of the track [9][10][11][12]. Additionally, the dynamic behaviour of the rail pad exhibits nonlinear aspects in terms of the frequency response function [3,[13][14][15][16][17][18].
As a component of the track, the fastening system containing the rail pad is a critical element which needs to be included in models when studying problems of practical interest, such as those related to rolling noise, the mitigation of the vibrations in soil, wear of rolling surfaces, etc.
Mechanical representations of tracks are based on the continuous medium theory, whose equations can be resolved by analytical, semi-analytical or numerical methods. In the case of analytical and semi-analytical methods, the rail can be modelled using Euler-Bernoulli beam theory [19,20] or Timoshenko beam theory [21,22]. There are two approaches: either the beam is finite, in which case the modal analysis method is applied to obtain time-domain simulations [23], or the beam is infinite, which has the advantage of eliminating the effect of the waves which are reflected from the edges of the model [24].
The main numerical methods applied for track modelling are the finite element method, the boundary element method and the discrete element method. The finite element method is recommended for track modelling in the domain of high frequencies, because this method offers more possibilities in terms of reproducing constructive details.    It follows from the above that rail pads have nonlinear, quasi-static elastic characteristics, as described in several works [3,5,6]. This nonlinear characteristic of the rail pad, together with the non-linear elastic characteristic of the ballast and subgrade [5,7,8], determine the specific shape of the load-deformation curve of the track [9][10][11][12]. Additionally, the dynamic behaviour of the rail pad exhibits nonlinear aspects in terms of the frequency response function [3,[13][14][15][16][17][18].
As a component of the track, the fastening system containing the rail pad is a critical element which needs to be included in models when studying problems of practical interest, such as those related to rolling noise, the mitigation of the vibrations in soil, wear of rolling surfaces, etc.
Mechanical representations of tracks are based on the continuous medium theory, whose equations can be resolved by analytical, semi-analytical or numerical methods. In the case of analytical and semi-analytical methods, the rail can be modelled using Euler-Bernoulli beam theory [19,20] or Timoshenko beam theory [21,22]. There are two approaches: either the beam is finite, in which case the modal analysis method is applied to obtain time-domain simulations [23], or the beam is infinite, which has the advantage of eliminating the effect of the waves which are reflected from the edges of the model [24].
The main numerical methods applied for track modelling are the finite element method, the boundary element method and the discrete element method. The finite element method is recommended for track modelling in the domain of high frequencies, because this method offers more possibilities in terms of reproducing constructive details. It follows from the above that rail pads have nonlinear, quasi-static elastic characteristics, as described in several works [3,5,6]. This nonlinear characteristic of the rail pad, together with the non-linear elastic characteristic of the ballast and subgrade [5,7,8], determine the specific shape of the load-deformation curve of the track [9][10][11][12]. Additionally, the dynamic behaviour of the rail pad exhibits nonlinear aspects in terms of the frequency response function [3,[13][14][15][16][17][18].
As a component of the track, the fastening system containing the rail pad is a critical element which needs to be included in models when studying problems of practical interest, such as those related to rolling noise, the mitigation of the vibrations in soil, wear of rolling surfaces, etc.
Mechanical representations of tracks are based on the continuous medium theory, whose equations can be resolved by analytical, semi-analytical or numerical methods. In the case of analytical and semi-analytical methods, the rail can be modelled using Euler-Bernoulli beam theory [19,20] or Timoshenko beam theory [21,22]. There are two approaches: either the beam is finite, in which case the modal analysis method is applied to obtain time-domain simulations [23], or the beam is infinite, which has the advantage of eliminating the effect of the waves which are reflected from the edges of the model [24].
The main numerical methods applied for track modelling are the finite element method, the boundary element method and the discrete element method. The finite element method is recommended for track modelling in the domain of high frequencies, because this method offers more possibilities in terms of reproducing constructive details. When the finite element method is used, the length of the model is always finite, and for this reason, the method exhibits the disadvantages associated with the finite length model mentioned above. The calculation time is considerably higher because the models have a very large number of finite elements and nodes [25,26]. A direction of research in which the finite element method is applied is that of rolling noise. In this case, the finite element method is supplemented with the boundary element method to calculate the sound pressure level [27,28].
The discrete element method subsumes any numerical method that allows the calculation of the motion of a large number of small bodies that are in contact or collide with each other or with a surface. The method has many applications, especially in problems involving granular materials such as the ballast. The use of this method is recommended to solve the problems concerning the specific phenomena of the ballast that lead to track unevenness, i.e., the migration of the ballast, resulting in hanging sleepers, and ballast settlement [29,30].
Regarding rail pad modelling, the simplest representations are linear models with concentrated parameters of the Kelvin-Voigt [25,31] or Poynting-Thompson [14] types, or with distributed parameters, i.e., of the Winkler type with viscous damping [23]. In [32,33], the rail and the sleepers are modelled using the finite element method, and the rail pad is modelled with the help of several discrete Kelvin-Voight systems that are regularly distributed on a line segment or on a rectangular surface. Obviously, these models do not reflect the nonlinearity of the load-deformation characteristic of the rail pad. Further, the loss factor is proportional to the frequency, which limits the frequency up to which they can be used. Better results from this point of view are obtained when hysteretic damping is introduced instead of viscous damping [34].
Other techniques with which to study the viscoelastic features of rail pads include the use of the fractional derivative model [35].
The nonlinear load-deformation characteristics of rail pads can be modelled using polynomial functions [36] or the bilinear function [37,38].
In this paper, the impact of the load-deformation characteristic of the rail pad on the wheel-rail vibration behaviour is treated. Specifically, we aim to highlight the impact of the variability of this characteristic, resulting from the manufacturing of the pads. For instance, Figure 3 shows the load-deformation characteristic for a certain rubber rail pad used at CFR (Romanian Railway), as well as the two limits, namely, the stiff limit and the soft limit, within which this characteristic must fall [39]; otherwise, the rail pad is rejected. It is observed that there is a wide margin between the two limit characteristics, and the question is asked: What is the potential impact on the wheel-rail vibration behaviour and, in turn, on the functionality of the rail pad? Obviously, narrowing the accepted tolerance margin leads to an increase in the number of rejected rail pads and to an increase in manufacturing costs, while widening this margin can have negative effects on the functionality of the rail pads.
To deal with this problem, a wheel-rail interaction model was developed in which the wheel is represented by a rigid body with one degree of freedom loaded by the vehicle weight on a wheel. This simple wheel model was chosen to limit the number of parameters in order to highlight more clearly the impact of the load-deformation characteristic of the rail pad. Such an approach is frequently encountered [20,22,23,34,37]. As for the track, our model considers the nonlinear load-deformation characteristics of the rail pad and ballast, as represented by bilinear functions [37,38]. The rail equilibrium position under a static load on the wheel was first determined and the elastic load-deformation characteristics of the rail were obtained. The wheel-rail interaction model is that of roughness displacement [22], and the rail receptance was calculated using a model with a nonhomogeneous foundation, whose parameters result from the equilibrium position. The impact of the variability of the rail pad load-deformation characteristic on wheel-rail vibrations in terms of dynamic contact force, wheel speed, and rail speed at the contact point was evaluated.  To deal with this problem, a wheel-rail interaction model was developed in which the wheel is represented by a rigid body with one degree of freedom loaded by the vehicle weight on a wheel. This simple wheel model was chosen to limit the number of parameters in order to highlight more clearly the impact of the load-deformation characteristic of the rail pad. Such an approach is frequently encountered [20,22,23,34,37]. As for the track, our model considers the nonlinear load-deformation characteristics of the rail pad and ballast, as represented by bilinear functions [37,38]. The rail equilibrium position under a static load on the wheel was first determined and the elastic load-deformation characteristics of the rail were obtained. The wheel-rail interaction model is that of roughness displacement [22], and the rail receptance was calculated using a model with a nonhomogeneous foundation, whose parameters result from the equilibrium position. The impact of the variability of the rail pad load-deformation characteristic on wheel-rail vibrations in terms of dynamic contact force, wheel speed, and rail speed at the contact point was evaluated.
To the knowledge of the authors, the influence of the variability of the elastic characteristic of the rail pad on wheel-rail vibrations has not been treated in the past in the manner described in this paper.

Hypotheses and Structure of the Mechanical Model
The modelling of a wheel-rail system to investigate the behaviour of the vertical vibrations is based upon several hypotheses that are widely encountered in the specialized literature. Thus, the vehicle and the running track are considered symmetrical in relation to the median longitudinal plane. Moreover, to operate railway vehicles, they must have balanced wheel loads, both when empty and when loaded. As a result, the two rails of the track are equally stressed during running; a similar observation can be made regarding the wheels of the vehicle. Under these conditions, it is expected that the wear of the running surfaces takes a reasonably similar shape on the two sides of the vehicle-track system. Consequently, the vehicle-track system can be reduced to half, taking the longitudinal vertical plane as the separation plane. In other words, vertical vibration studies consider a 2D mechanical model that represents half of the vehicle and one rail.
The following assumption refers to the frequency domain, in which the natural frequencies of both the vehicle and the wheel-rail system are located. The vehicle's natural frequencies are between 0.5-20 Hz, and the corresponding vehicle vibration modes are excited by irregularities in the track superstructure with wavelengths between 3 and 100-120 m. The wheel-rail natural frequencies are greater than 30-40 Hz, and the vibrations, which can reach 2-3000 Hz, are induced by the irregularities of the rolling surfaces, with To the knowledge of the authors, the influence of the variability of the elastic characteristic of the rail pad on wheel-rail vibrations has not been treated in the past in the manner described in this paper.

Hypotheses and Structure of the Mechanical Model
The modelling of a wheel-rail system to investigate the behaviour of the vertical vibrations is based upon several hypotheses that are widely encountered in the specialized literature. Thus, the vehicle and the running track are considered symmetrical in relation to the median longitudinal plane. Moreover, to operate railway vehicles, they must have balanced wheel loads, both when empty and when loaded. As a result, the two rails of the track are equally stressed during running; a similar observation can be made regarding the wheels of the vehicle. Under these conditions, it is expected that the wear of the running surfaces takes a reasonably similar shape on the two sides of the vehicle-track system. Consequently, the vehicle-track system can be reduced to half, taking the longitudinal vertical plane as the separation plane. In other words, vertical vibration studies consider a 2D mechanical model that represents half of the vehicle and one rail.
The following assumption refers to the frequency domain, in which the natural frequencies of both the vehicle and the wheel-rail system are located. The vehicle's natural frequencies are between 0.5-20 Hz, and the corresponding vehicle vibration modes are excited by irregularities in the track superstructure with wavelengths between 3 and 100-120 m. The wheel-rail natural frequencies are greater than 30-40 Hz, and the vibrations, which can reach 2-3000 Hz, are induced by the irregularities of the rolling surfaces, with wavelengths between 2-4 cm and 3 m. From the two facts presented above, it follows that the vehicle vibration and wheel train vibration are decoupled, and that the wheels vibrate independently of each other, because the coupling resulting from the suspension is unimportant. The vehicle model is thus reduced to a train of loaded wheels running on a rail; the load is actually the static load corresponding to a wheel.
The third assumption is that the rail bending waves generated by the vibration of a wheel-rail pair will have a negligeable effect upon the vibration of other wheel-rail pair. Indeed, depending on the frequency, the rail bending waves have two or one evanescent waves and zero or one propagative wave. Evanescent waves are strongly attenuated by their very nature, and the attenuation of the propagative waves is reduced as the frequency increases. Since the present study deals with low and medium frequencies, it follows from the above that the coupling between the vehicle wheels due to rail bending waves can be neglected.
Adopting the above assumptions, the mechanical model shown in Figure 4 is proposed to analyze the wheel-rail vibrations. The vehicle model is reduced to a wheel which is considered as a rigid body of mass (M w ) loaded with the static load Q, and that of the track is simplified to a rail represented by an infinitely uniform Euler-Bernoulli beam on an elastic two-layer foundation between which an inertial layer is interposed [40]. The elastic layers model the effect of the rail pads and the ballast, while the inertial layer models the influence of the sleepers on the rail vibrations. The frequency domain of this type of model is limited to 6-700 Hz, because the effect of the spacing of the sleepers on the rail bending vibrations (the pinned-pinned mode) is neglected. The track model parameters are the bending stiffness, EI, and the mass per length unit m of the rail, the stiffness per length unit k 1e and k 1r for the rail pad, the mass per length unit m s of the sleepers, and the stiffness per length unit k 2e and k 2r for the ballast. Stiffness k 1e , k 1s , k 2e and k 2s are determined by bilinear approximations of the nonlinear load-displacement characteristics of the rail pad and ballast. Distances l 1 and l 2 are determined from the equilibrium condition under static load Q.
waves and zero or one propagative wave. Evanescent waves are strongly attenuated by their very nature, and the attenuation of the propagative waves is reduced as the frequency increases. Since the present study deals with low and medium frequencies, it follows from the above that the coupling between the vehicle wheels due to rail bending waves can be neglected.
Adopting the above assumptions, the mechanical model shown in Figure 4 is proposed to analyze the wheel-rail vibrations. The vehicle model is reduced to a wheel which is considered as a rigid body of mass (Mw) loaded with the static load Q, and that of the track is simplified to a rail represented by an infinitely uniform Euler-Bernoulli beam on an elastic two-layer foundation between which an inertial layer is interposed [40]. The elastic layers model the effect of the rail pads and the ballast, while the inertial layer models the influence of the sleepers on the rail vibrations. The frequency domain of this type of model is limited to 6-700 Hz, because the effect of the spacing of the sleepers on the rail bending vibrations (the pinned-pinned mode) is neglected. The track model parameters are the bending stiffness, EI, and the mass per length unit m of the rail, the stiffness per length unit k1e and k1r for the rail pad, the mass per length unit ms of the sleepers, and the stiffness per length unit k2e and k2r for the ballast. Stiffness k1e, k1s, k2e and k2s are determined by bilinear approximations of the nonlinear load-displacement characteristics of the rail pad and ballast. Distances l1 and l2 are determined from the equilibrium condition under static load Q.   Figure 5 displays the stiff and soft limits of the load-displacement characteristics of the rail pad and the mean characteristic calculated as the arithmetic mean of the limit characteristics. Typically, these three characteristics are associated with stiff, medium, and soft rail pads. Additionally, the nonlinear load-displacement characteristic of the ballast per semi-sleeper is presented. This characteristic obeys Hertz's law [5]

Bilinear Approximation of the Load-Displacement Characteristics of the Rail Pad and Ballast
where Q is the load, u b is ballast deformation, and C H is the Hertz constant. Each nonlinear characteristic is approximated by a bilinear function which has elastic and rigid parts ( Figure 6) where k rr and k re are the rigid and elastic parts of the bilinear function of the rail pad, k br and k be are the rigid and elastic parts of the bilinear function of the ballast, u r and u b are the rail pad and ballast deformation due to loads Q(u r ) and Q(u b ), respectively, u ro and u bo are the rail pad and ballast deformation when it is passing from the elastic to the rigid part, and u rm and u bm are the rail pad and ballast deformation for the maximum value of the load. Applying the least mean squares method, we obtain the parameters of the bilinear function.  Figure 5 displays the stiff and soft limits of the load-displacement characteristics of the rail pad and the mean characteristic calculated as the arithmetic mean of the limit characteristics. Typically, these three characteristics are associated with stiff, medium, and soft rail pads. Additionally, the nonlinear load-displacement characteristic of the ballast per semi-sleeper is presented. This characteristic obeys Hertz's law [5] 3 ,

Bilinear Approximation of the Load-Displacement Characteristics of the Rail Pad and Ballast
where Q is the load, ub is ballast deformation, and CH is the Hertz constant. Figure 5. Nonlinear characteristic of the rail pad (based on [40]) and ballast (based on [5]).
Each nonlinear characteristic is approximated by a bilinear function which has elastic and rigid parts ( Figure 6) where krr and kre are the rigid and elastic parts of the bilinear function of the rail pad, kbr and kbe are the rigid and elastic parts of the bilinear function of the ballast, ur and ub are the rail pad and ballast deformation due to loads Q(ur) and Q(ub), respectively, uro and ubo are the rail pad and ballast deformation when it is passing from the elastic to the rigid part, and urm and ubm are the rail pad and ballast deformation for the maximum value of the load. Applying the least mean squares method, we obtain the parameters of the bilinear function.   Table 1 presents the values of the parameters of the bilinear functions.   Table 1 presents the values of the parameters of the bilinear functions. Table 1 also contains the values of the stiffnesses of the two elastic layers calculated depending on a sleeper bay of 0.577 m:   Table 1 presents the values of the parameters of the bilinear functions.
where k1r and k1e are the stiffness of the rigid/elastic part of the bilinear function of the first elastic layer (rail pad), k2r and k2e are the stiffness of the rigid/elastic part of the bilinear function of the second elastic layer (ballast), and d is the sleeper bay. There are very large differences in the stiffness of the elastic part of the bilinear function associated with the load-deformation characteristic of the rail pad and small variations in the stiffness of its rigid part. To evaluate the error introduced by the bilinear function and the range of applicability, the track equivalent characteristics must be calculated. This represents the uniform distributed load, q, depending on the rail displacement ( Figure 8).  To evaluate the error introduced by the bilinear function and the range of applicability, the track equivalent characteristics must be calculated. This represents the uniform distributed load, q, depending on the rail displacement ( Figure 8).
where w is the rail displacement and the equivalent stiffnesses are where k 1,2,3 are the equivalent stiffnesses corresponding to the elastic, moderate, and stiff portions of the uniform distributed load-rail displacement characteristic, and the rail displacement for the transition points, w 12 and w 23 , are given by: with w 12 < w 23 , because this considers that (see Figures 6 and 7): By imposing the maximum admissible error, we obtain the maximum rail displacement. For instance, considering a maximum error of 5%, Figure 9 presents the distributed load versus the rail displacement for both the theoretical nonlinear characteristics and a bilinear approximation for the three characteristics of the rail pad. The main results are summarized in Table 2. The Rms errors are acceptable, i.e., smaller than 1.91%.
By imposing the maximum admissible error, we obtain the maximum rail displacement. For instance, considering a maximum error of 5%, Figure 9 presents the distributed load versus the rail displacement for both the theoretical nonlinear characteristics and a bilinear approximation for the three characteristics of the rail pad. The main results are summarized in Table 2. The Rms errors are acceptable, i.e., smaller than 1.91%.

Equilibrium Position of the Rail
In this section, the rail equilibrium position under a static load is determined. There are two reasons for this: (a) to check the applicability domain of the bilinear approximation in the sense that the maximum displacement of the rail is not higher than the limit of the applicability determined in the previous section (Table 2); and (b) to identify the transition sections of the elastic foundation, i.e., to calculate the distances l1 and l2, respectively. Figure 10 presents the calculation model derived from the one presented in Figure 4; the inertial layer of the sleepers is not displayed for reasons of simplicity, because its position is not relevant here.

Equilibrium Position of the Rail
In this section, the rail equilibrium position under a static load is determined. There are two reasons for this: (a) to check the applicability domain of the bilinear approximation in the sense that the maximum displacement of the rail is not higher than the limit of the applicability determined in the previous section (Table 2); and (b) to identify the transition sections of the elastic foundation, i.e., to calculate the distances l 1 and l 2 , respectively.  Figure 4; the inertial layer of the sleepers is not displayed for reasons of simplicity, because its position is not relevant here.

Equilibrium Position of the Rail
In this section, the rail equilibrium position under a static load is determined. There are two reasons for this: (a) to check the applicability domain of the bilinear approximation in the sense that the maximum displacement of the rail is not higher than the limit of the applicability determined in the previous section (Table 2); and (b) to identify the transition sections of the elastic foundation, i.e., to calculate the distances l1 and l2, respectively. Figure 10 presents the calculation model derived from the one presented in Figure 4; the inertial layer of the sleepers is not displayed for reasons of simplicity, because its position is not relevant here. The equilibrium equation of the rail under a static load can be refined from the two equilibrium equations of the track by eliminating the inertial layer displacement as: The equilibrium equation of the rail under a static load can be refined from the two equilibrium equations of the track by eliminating the inertial layer displacement as: where w(x) is the rail displacement, k(x) is the equivalent stiffness of the two elastic layers depending on coordinate x along the rail, and q(x) is the distributed load supported by the rail with q(x) = Qδ(x), where δ(.) is the Dirac delta function.
The equivalent stiffness has the following form: where ±l 1 and ±(l 1 + l 2 ) are the coordinates along the rail corresponding to the transition sections of the uniform distributed load-rail displacement characteristic. The rail displacement at the transition section, i.e., x = ±l 1 and x = ± (l 1 + l 2 ), takes the values: w(±l 1 ) = w 12 , w(±(l 1 + l 2 )) = w 23 .
The boundary condition associated with Equation (9) is: Considering only a half rail due to symmetry and applying the direct method, which implies that the rail to be segmented is in the concentrated load section and in the transition sections, the equations of equilibrium are obtained: where displacement w i (x i ) is the displacement of the i rail segment in section x i , which is related to rail displacement w(x) as follows: The boundary conditions associated with Equation (13) are as follows: at x 3 = 0, the rail slope is null and the shear force equals -Q/2 dw 3 (0) at x 3 = l 1 and x 2 = 0, the rail displacement is w 23 and the first three derivatives of the rail displacement are continuous functions: at x 2 = l 2 and x 1 = 0, the rail displacement is w 12 and the continuity conditions must be met as above: The solutions to differential Equation (13) are: where A n with n =1÷9 are constants to be determined, and We observe that the shape of w 1 (x 1 ) accomplishes the boundary condition at x 1 = 0 and for x 1 → ∞ .
Equation (20) is inserted in the boundary conditions (15)- (19), resulting 11 non-linear equations of unknowns A n , n = 1 ÷ 9 and l 1 and l 2 . Applying the Newton-Raphson algorithm, we obtain the unknown values, including the distances l 1 and l 2 . Finally, the load-rail displacement characteristic can be drawn, as seen in Figure 11. Under a static load of 100 kN, the rail displacement is 1.465 mm for the soft rail pad, 1.352 mm for the medium rail pad, and 1.239 mm when the rail pad is stiff; in all cases, the limit of applicability of the bilinear approximation is not exceeded.
and the continuity conditions must be met as above: The solutions to differential Equation (13) where An with n =1÷9 are constants to be determined, and We observe that the shape of w1(x1) accomplishes the boundary condition at x1 = 0 and for x1 →∞.
Equation (20) is inserted in the boundary conditions (15)- (19), resulting 11 non-linear equations of unknowns An, n = 1 ÷ 9 and l1 and l2. Applying the Newton-Raphson algorithm, we obtain the unknown values, including the distances l1 and l2. Finally, the loadrail displacement characteristic can be drawn, as seen in Figure 11. Under a static load of 100 kN, the rail displacement is 1.465 mm for the soft rail pad, 1.352 mm for the medium rail pad, and 1.239 mm when the rail pad is stiff; in all cases, the limit of applicability of the bilinear approximation is not exceeded. The results regarding distances l1 and l2 are listed in Table 3. The results regarding distances l 1 and l 2 are listed in Table 3.

Equations of the Dynamic Behavior
Wheel-rail vibrations develop around the equilibrium position and may be described by the displacement of the wheel, ∆u(t), the displacement of the rail, ∆w(x,t), and the displacement of the sleepers, ∆z(x,t).
The equations of motion are as follows: for the wheel for the track where k I,II (x) is the stiffness of the i elastic layer (i = 1,2) k 2r , |x| ≤ l 1 ; k 2r , l 1 < |x| ≤ l 2 ; k 2e , l 1 + l 2 < |x|, (24) and ∆Q(t) is the dynamic component of the contact force.
The following contact equation should be considered: where k H is the contact stiffness, and r(t) is the irregularity of the rolling surfaces.
Considering the harmonic steady-state behavior, the equations of motion can be written as: and the contact equation as: where ∆u,∆w(x),∆z(x), ∆Q, and r are the complex amplitudes of the wheel, rail, sleeper, contact force, and roughness, and where η I,II is the loss factor of the first and second elastic layer, and η H is the loss factor of the elastic contact, and i 2 = −1. Equation (27) is equivalent to where In fact, the complex parameter K(x, ω) has a similar shape to those in Equation (24): Following a similar method as in the equilibrium case, the following differential equations are considered: The characteristic equations associated to the above differential equations are: that gives the eigenvalues λ 1,3 = ±(a n − ib n ) λ 2,4 = ±(b n + ia n ) (37) where a n and b n are positive real numbers. The complex amplitude of the rail is: where W n with n = 1÷10 are constants to be determined. The last form meets the boundary condition for x → ∞ .
The following boundary conditions hold: at x = 0, the rail slope is null due to the symmetry and the shear force amplitude equalling −∆Q/2 at x = l 1 , the continuity conditions are imposed: where the left parts are calculated from Equation (38) and the right parts are calculated from Equation (39).
at x = l 1 + l 2 , the continuity conditions are imposed where the left parts are calculated from Equation (39) and the right parts from Equation (40). Finally, the following set of linear algebraic equations are obtained from Equations (38)-(43) where l = l 1 + l 2 .
The solution to the above system can be obtained numerically, depending on contact force complex amplitude ∆Q. Then, the amplitude of the rail is calculated by inserting W 1÷10 into Equations (38)- (40).
Rail receptance is the ratio between the rail amplitude and contact force amplitude: Taking ∆Q = 1, the rail receptance may be determined using Equations (44)-(54). Similarly, the wheel receptance is the ratio between the wheel amplitude and the force amplitude where the force and the wheel displacement have the same orientation. Considering the above, the wheel/rail vibrations are described by the following equations: where G H (ω) = 1/k H is the wheel/rail contact receptance. Generally, the receptance is a complex number depending on the angular frequency, and because of this, all receptances have been notated corresponding to these features. However, the wheel receptance is a real number and the wheel/rail contact receptance does not depend on angular frequency according to the model.
The solution to Equations (56)-(58) is: where are the frequency response functions.
When it is assumed that the roughness of the rolling surfaces is a stationary and ergodic process described by the power spectral density S r (Ω), where Ω is the wavenumber, then the power spectral density depending on the angular frequency ω = VΩ is: The power spectral density of the output quantity p (∆u, ∆w and ∆Q) that describes the vibration behavior is: where H p (ω) is the frequency response function associated with the quantity p, and the rms value calculated within an angular frequency interval is: where ω i and ω s are the interval limits and ω c the central angular frequency.

Numerical Application
In this section, the wheel-rail model with a nonhomogeneous foundation is used to identify the impact of the variability of the rail pad elastic characteristics upon the wheel-rail vibrations, employing the MATLAB environment. Table 3 includes the values of the wheel-rail model parameters. The track is fitted with UIC 60 rail type, concrete sleepers of 268 kg, and the sleeper spacing is 577 mm. Stiffness values of the two elastic layers are the result of approximating the characteristics of the rail pad and the ballast with the help of the bilinear function, as shown in Section 2.1. The wheel mass and static load correspond to an electric locomotive. Figure 12 shows the receptance modulus of the rail at the active point, i.e., where the harmonic force acts, as calculated for the homogenous two-layer model (Figure 13), where the stiffness of the elastic layers corresponds to the medium characteristic case (rigid part-k 1r = 443.464 MN/m 2 , k 2r = 96.681 MN/m 2 ). This stiffness configuration is typically for the zone around the wheel-rail contact point, as shown in Figure 4. Considering the undamped case, the rail receptance exhibits two resonance frequencies, a low resonance at 91 Hz and a high one at 482 Hz. Between them, an antiresonance frequency at 240 Hz occurs due to the dynamic absorber effect of the sleepers. Conventionally, there are three ranges of frequency: the low frequency range, lower than the low resonance frequency, the medium frequency range between the two resonance frequencies, and the high frequency range located at frequencies higher than the high resonance frequency. The response of the rail is the effect of the dynamic balance of the forces acting on it, i.e., the elastic force, inertial force, and the excitation force. Thus, the rail receptance is almost constant at low frequency, because the elastic force is much greater than the inertial force, while at high frequency, the rail receptance decreases with frequency as the inertial force becomes dominant. Figure 14 presents the rail receptance calculated for the nonhomogeneous two-layer track model for the three cases considered: soft, medium, and stiff characteristic. Additionally, the undamped case was considered. The rail receptance diagrams exhibit similar peaks, corresponding to the resonance frequencies indicated in the case of the homogeneous two-layer foundation, with small differences due to the differences between the stiffness values of the foundation around the active point in the three cases. However, these two peaks are flattened as an effect of the direct/reflected wave interference. In the range of low and high frequencies, the rail receptance shows similar characteristics to those highlighted in Figure 12. In the range of medium frequencies, the appearance of the rail receptance is different, being marked by abrupt variations due to the standing waves that are formed because of the overlapping of the direct bending waves with the reflected ones caused by the presence of the sections in which the foundation stiffness changes. the medium frequency range between the two resonance frequencies, and the high frequency range located at frequencies higher than the high resonance frequency. The response of the rail is the effect of the dynamic balance of the forces acting on it, i.e., the elastic force, inertial force, and the excitation force. Thus, the rail receptance is almost constant at low frequency, because the elastic force is much greater than the inertial force, while at high frequency, the rail receptance decreases with frequency as the inertial force becomes dominant.   the medium frequency range between the two resonance frequencies, and the high frequency range located at frequencies higher than the high resonance frequency. The response of the rail is the effect of the dynamic balance of the forces acting on it, i.e., the elastic force, inertial force, and the excitation force. Thus, the rail receptance is almost constant at low frequency, because the elastic force is much greater than the inertial force, while at high frequency, the rail receptance decreases with frequency as the inertial force becomes dominant.    Figure 14 presents the rail receptance calculated for the nonhomogeneous two-laye track model for the three cases considered: soft, medium, and stiff characteristic. Addi tionally, the undamped case was considered. The rail receptance diagrams exhibit simila peaks, corresponding to the resonance frequencies indicated in the case of the homogene ous two-layer foundation, with small differences due to the differences between the stiff ness values of the foundation around the active point in the three cases. However, thes two peaks are flattened as an effect of the direct/reflected wave interference. In the rang of low and high frequencies, the rail receptance shows similar characteristics to thos highlighted in Figure 12. In the range of medium frequencies, the appearance of the ra receptance is different, being marked by abrupt variations due to the standing waves tha are formed because of the overlapping of the direct bending waves with the reflected one caused by the presence of the sections in which the foundation stiffness changes.  Figure 15 shows the impact of hysteretic damping on rail receptance. All abrupt var iations of the rail receptance in the medium frequency range are removed due to the hys teretic damping. The two flattened peaks corresponding to the low and high resonanc frequency and the depth characteristic of the anti-resonant behavior remain the same. Rel  Figure 15 shows the impact of hysteretic damping on rail receptance. All abrupt variations of the rail receptance in the medium frequency range are removed due to the hysteretic damping. The two flattened peaks corresponding to the low and high resonance frequency and the depth characteristic of the anti-resonant behavior remain the same. Relative significant differences are maintained between the three cases considered in the field of low and average frequencies.  Figure 16 displays the frequency response function (F.R.F.) for the wheel displacement, rail displacement, and the contact force. At low frequency, the wheel takes the displacement imposed by the roughness of the rolling surfaces because the rail behaves as a rigid element and, due to this fact, the value of the F.R.F. module is around 1. The maximum value recorded at approx. 44 Hz is the effect of the resonance of the wheel on the track. At medium and high frequencies, the wheel displacement decreases because its receptance decreases as the frequency increases (see Equation (55)), and the rail takes the displacement imposed by the irregularity of the rolling surface, as can see in Figure 16b, where the F.R.F. modulus of the rail displacement is around 1. Regarding the response function of contact forces, this has a general increasing tendency depending on the frequency. It presents maxima relative to the wheel/track resonance frequency and to the anti-resonance frequency of the rail. At the two resonance frequencies of the rail, the contact force shows local minima. The impact of the rail pad elastic characteristic manifests itself mainly at medium frequencies and affects the wheel vibration and the magnitude of the contact force. Meanwhile, the rail vibration is less influenced by the variability of the elastic characteristic of the rail pad.  Figure 16 displays the frequency response function (F.R.F.) for the wheel displacement, rail displacement, and the contact force. At low frequency, the wheel takes the displacement imposed by the roughness of the rolling surfaces because the rail behaves as a rigid element and, due to this fact, the value of the F.R.F. module is around 1. The maximum value recorded at approx. 44 Hz is the effect of the resonance of the wheel on the track. At medium and high frequencies, the wheel displacement decreases because its receptance decreases as the frequency increases (see Equation (55)), and the rail takes the displacement imposed by the irregularity of the rolling surface, as can see in Figure 16b, where the F.R.F. modulus of the rail displacement is around 1. Regarding the response function of contact forces, this has a general increasing tendency depending on the frequency. It presents maxima relative to the wheel/track resonance frequency and to the anti-resonance frequency of the rail. At the two resonance frequencies of the rail, the contact force shows local minima. The impact of the rail pad elastic characteristic manifests itself mainly at medium frequencies and affects the wheel vibration and the magnitude of the contact force. Meanwhile, the rail vibration is less influenced by the variability of the elastic characteristic of the rail pad. Figure 17 shows the spectrum of the rms value calculated for the wheel speed and rail speed at the active contact and for the contact force when the wheel velocity is 160 km/h. The spectrum is calculated in 1/3 octave intervals. The vibration speed was chosen as the parameter of interest because the acoustic power of the source depends on it, and the rolling noise is an important effect of wheel-rail vibrations. The data in Figure 17 are supplemented with those in Figure 18, where the ratio of the spectral components is presented.
It is observed that the variability of the elastic characteristic of the rail pad influences the vibration of the wheel and the contact force in almost the entire considered frequency range. In most cases, the vibration of the wheel intensifies or decreases if the characteristic of the rail pad is stiff or soft. In general, the variations are smaller by a few percent, but there are also situations, especially with medium frequencies, where larger differences are observed. As for the rail, at very low and high frequencies, its vibrations are more intense with a soft rail pad. Otherwise, the rail vibration increases with a stiffer rail pad. In any case, the impact of the variability of the rail pad characteristic is less visible in the case of the rail. (c) contact force; color code as in Figure 15. Figure 17 shows the spectrum of the rms value calculated for the wheel speed and rail speed at the active contact and for the contact force when the wheel velocity is 160 km/h. The spectrum is calculated in 1/3 octave intervals. The vibration speed was chosen as the parameter of interest because the acoustic power of the source depends on it, and the rolling noise is an important effect of wheel-rail vibrations. The data in Figure 17 are supplemented with those in Figure 18, where the ratio of the spectral components is presented.
It is observed that the variability of the elastic characteristic of the rail pad influences the vibration of the wheel and the contact force in almost the entire considered frequency range. In most cases, the vibration of the wheel intensifies or decreases if the characteristic of the rail pad is stiff or soft. In general, the variations are smaller by a few percent, but there are also situations, especially with medium frequencies, where larger differences are observed. As for the rail, at very low and high frequencies, its vibrations are more intense with a soft rail pad. Otherwise, the rail vibration increases with a stiffer rail pad. In any case, the impact of the variability of the rail pad characteristic is less visible in the case of the rail. (c) contact force; ᵒ, stiff characteristic, ᵒ, medium characteristic; ᵒ, soft characteristic; color code as in Figure 15. Figure 17. Spectrum of rms value: (a) wheel speed; (b) rail speed; (c) contact force; • , stiff characteristic, • , medium characteristic; • , soft characteristic; color code as in Figure 15. Figure 17. Spectrum of rms value: (a) wheel speed; (b) rail speed; (c) contact force; ᵒ, stiff characteristic, ᵒ, medium characteristic; ᵒ, soft characteristic; color code as in Figure 15.

Conclusions
Although at first glance, it seems to be a simple small piece of elastic material, the rail pad is an essential component of the track superstructure due to the role it plays in terms of vibration damping and reducing the shocks transmitted to the sleepers and ballast.
The functionality of the rail pad depends on its most important property, namely, its load-deformation characteristics, which have a specific, nonlinear shape, and which present significant variability in mass production.
In this paper, the impact of the variability of the load-deformation characteristics of the rail support on the wheel-rail vibration behavior is studied. To this end, the nonlinear loaddeformation characteristics have been approximated with the help of the bilinear function and implemented in a track model with a nonhomogeneous two-layer elastic foundation.
The advantage of this method is that the model of the track with a nonhomogeneous foundation with continuous variation due the nonlinear characteristic of the rail pad and ballast is transformed into a model with nonhomogeneous foundation with two-step variation. In other words, the track model has three distinct portions with a homogeneous foundation, which simplifies the mathematical solution of the equations that describe the equilibrium state and the dynamic behavior.
The load-deformation characteristic of the rail pad shows great variability in terms of the stiffness of the elastic part, while the stiffness of the rigid part is much more stable.
The rail pad works under the static load on the rigid side of the bilinear loaddeformation characteristic, which has small deviations; this explains the relatively low impact of the variability of the rail pad characteristic on the wheel-rail vibrations.
The results show that the impact of the variability of the load-deformation characteristic of the rail pad is most pronounced on the wheel-rail contact force and on the wheel vibrations, especially under the resonance frequencies of the wheel on the track and in the mid-frequency range. In contrast, rail vibrations have less influence.
Future research will focus on expanding the frequency domain of the investigation by adapting the bilinear characteristic of the rail pad to the track model with discrete supports.