1 Introduction

The study of deterioration mechanisms occurring at the interface between rail head and wheel tread often involves the solution of the associated contact problem. Determining the size and shape of the contact area is obviously a necessary condition to find the pressure distribution that characterizes the solution [1]. Interestingly, if both contact partners are considered as elastic half-spaces of the same material, tangential tractions do not disturb the normal pressure distribution; hence, the shape and size of the contact area depend exclusively on the surface profiles and the normal force [2].

Standard wheel and rail profiles are usually composed of circular arcs with different radii and can also exhibit a continuously variable curvature [3, 4]. The complicated geometry rules out the application of the Hertzian theory to obtain an accurate analytical solution. However, the interface has been thoroughly analyzed by means of numerical procedures with predominant presence of boundary and finite element methods [1, 5]. On the one hand, numerical simulation of wheel-rail contact as a key interface in vehicle dynamics has been strongly influenced by the pioneering work of Kalker, who employed boundary element analysis with a half-space assumption to seek an efficient solution to the rolling contact problem [6,7,8]. Later developments aimed to reduce computational cost requirements, as done for instance by Polach [9], as well as overcome initial concerns related to the poor performance of boundary element analysis in non-Hertzian scenarios [10, 11], resulting in reliable and accurate predictions also in two-point and conformal contact problems [12, 13]. On the other hand, finite element analysis remains the preferred choice as regards general applicability, especially concerning the study of localized structural phenomena, such as crack initiation [14] and crack growth [15].

Despite the fact that the predictive capability of a computational model needs to be validated by comparison between simulation and experimental outcomes [16], the contrast between the highly-developed numerical procedures and the few experimental solutions available to determine the shape and size of the wheel-rail contact area is remarkable. In fact, little research has been done in the field since Labrijn introduced carbon paper between wheel and rail in the fifties [17]. An interesting approach consists of letting air flow through holes drilled on a rail block and deducing the contact area from the pressure variation in the supply line [18]. Another possibility is to hit the interface with ultrasonic waves and detect waves reflected in the contactless domain due to the low acoustic impedance of air [19].

In recent years, monosheet Prescale pressure measurement films [20] have joined the family of experimental techniques successfully applied to wheel-rail contact area measurements [21, 22]. Even though limited to static contact, the nondestructive character of this procedure makes it an ideal complement to other methods. The films are composed of a polyethylene terephthalate (PET) layer and a thin color-developing layer coated with microcapsules that release dye when inserted between contact bodies if a predetermined pressure is exceeded. The three film types, labeled as MS, HS and HHS, have a common PET layer but varying capsule sizes and wall thicknesses that lead to measurement ranges of 10–50 MPa, 50–130 MPa and 130–300 MPa respectively. The stain caused by the released fluid can be digitally post-processed to obtain a faithful depiction of the pressure distribution within said ranges. Maximum pressure values at the wheel-rail interface exceed considerably the measurable ranges, which makes Prescale films unsuitable for pressure distribution analysis. Moreover, contact area measurements would require a lower pressure limit of 0 MPa. However, the markedly steep pressure gradient near the boundary initially suggested that the imprint could be assumed to be in very close agreement with the contact domain.

Said assumption was refuted by the authors in the context of previous research as each film type measured substantially different contact areas with identical experimental setups [23]. In order to quantify the measurement error by comparison with analytical solutions, the three film types were employed to measure the contact area between test bodies that fulfilled the Hertzian assumptions under experimental conditions representative of wheel-rail contact. More precisely, three steel spherical specimens with radii R of 300 mm, 450 mm and 600 mm and smooth surface finish (Ra < 0.1 µm) were pushed one after the other against the flat surface of another steel specimen with normal loads of 25 kN, 50 kN, 75 kN and 100 kN (the last two were skipped for the R300 body in order to avoid plastic strain) leading to a circular contact area [24]. As shown in Table 1, all three types overestimated the radius of the resulting contact patch r in all cases. Interestingly, the error intrinsic to each film type was very similar regardless of geometry and load variations. This and the excellent reproducibility of the experiments led to the determination of film-dependent constant corrections. In addition, the effect of surface roughness was investigated by gradually roughening the originally smooth plane surface up to an approximate Ra value of 1.5 µm so as to consider the typical roughness values present in service [25]. Whereas MS and HS films detected a likewise gradual increase in nominal contact area size, the HHS type was not sensitive to roughness variations within the investigated range.

Table 1 Absolute observational error of contact patch radius r in mm measured with MS, HS and HHS Prescale films as compared to Hertzian solution
Full size table

Subsequent considerations concerning the mechanical alteration caused by the presence of Prescale film between contact surfaces shed light on the dependence of measurement error and roughness sensitivity on the film type. By generating detailed finite element models of the Hertzian experiments including the film, a relation between measurement error and capsule size (and thus film type) could be established [26]. Furthermore, microscopic roughness was implemented in the macroscopic models of the experiments in order to demonstrate that the lack of sensitivity of HHS films to roughness variations is perfectly consistent with the governing contact mechanics [27]. Said studies considered the plastic nature of PET by defining linear Drucker-Prager material behavior [28], which was originally conceived for soil materials but has found wide use in polymer modelling as it considers their characteristic dependence of yield stress on hydrostatic pressure [29]. Even though the film model proved useful to approximate the measurement trends exhibited by real films, it has not been tuned to reproduce measured values accurately and is thus lacking calibration.

After quantifying the contact area deviation originated in the Hertzian experiments and elucidating the mechanical principles underlying the imprint formation on the film surface, the question arose as to how these findings are to be transferred to actual wheel-rail contact. Since contact conditions might differ considerably from the Hertzian scenario, no analytical benchmark for comparison purposes exists and, moreover, the irregular shape of the contact area hinders considerably the determination of correction factors of general applicability that filter out the measurement error. The solution proposed here is to develop an efficient phenomenological model of the film by subjecting it to empirical calibration with the only requirement of matching the contact areas obtained from representative Hertzian experiments (including the error values in Table 1). Then, the model can be incorporated alternately into a wheel-rail contact finite element model in order to estimate the measurement error by comparison between with-film and without-film scenarios.

2 Materials and methods

2.1 Film model calibration

There are compelling reasons not to include the microscopic measurement layer into macroscopic finite element models, namely, its relative thinness and associated aspect ratio considerations [30]. This limitation is especially restrictive for three-dimensional models and becomes prohibitive for element types that require at least a few elements through the thickness to achieve an accurate solution. Therefore, a film model meant for three-dimensional applications such as wheel-rail contact must be computationally inexpensive. This can be achieved with an appropriate choice of element type and material model. The modelling strategy presented here suggests circumventing computational cost restrictions by employing gasket elements with linear elastic material behavior.

Gasket elements are special-purpose finite elements specifically formulated for compressive applications that conveniently allow for a single-element-layer discretization [31]. Since gaskets are highly compressible mechanical seals that deform mainly in the thickness direction, membrane and transverse shear stiffness are sometimes neglected in related material models [32]. Common formulations enable incorporating uncoupled thickness, membrane and transverse behaviors, and derive strains from the relative displacement between the top and bottom surfaces of the elements. Alternatively, a conventional material model can be assigned to the elements by defining a constitutive law [31], as done in this work.

The finite element model employed so far to replicate the Hertzian experiments conducted in the laboratory [26, 27] was accordingly adapted for the material calibration procedure as depicted in Fig. 1. The model was composed of axisymmetric cross-sections of the spherical and flat contact bodies with dimensions of 50 × 50 mm as well as Prescale film interposed in between. Both contact partners remained divided into tied contactless and contact regions meshed with reduced-integration, quadrilateral elements with approximate element sizes of 1.25 mm and 100 μm, respectively. The finely meshed film model used in previous studies, however, was considerably simplified into a single layer of gasket elements. The thickness dimension was set to 90 μm based on experimental observations [26] and a radial dimension of 100 μm was defined for each element so as to match the measurement resolution [33]. The bottom nodes of the flat body were vertically constrained and the normal force was introduced as pressure on the spherical body. Surface-to-surface contact discretization with hard pressure-overclosure behavior implementing augmented Lagrange constraint enforcement and a coefficient of friction of 0.4 were used for both contact interfaces. Additionally, separation after contact in the normal direction was prevented at the flat interface, which increased considerably convergence rates.

Fig. 1
figure 1

Finite element representation of the Hertzian experiments with single-layer gasket element Prescale film for the R300 case

Full size image

Steel bodies were modeled with linear elastic material parameters E = 210,000 MPa and ν = 0.3. As for the film, the efficiency requirement ruled out plasticity models considered so far in favor of linear elasticity. In fact, finite element models of Prescale films have traditionally been linear-elastic, although parameter values are apparently disputed. Hale and Brown determined a Young’s modulus E on the order of 100 MPa and a Poisson’s ratio ν of 0.45 by means of compressive tests [30] and these values have been extensively assigned to finite element representations of the film [34,35,36]. On the one hand, this conflicts with the Young’s modulus listed for PET in the specialized literature, which amounts to 3000 MPa [37]. As a matter of fact, it was concluded from comparative finite element simulations that the Young’s modulus of Prescale films should be near that value [38]. On the other hand, the Poisson’s ratio of PET ranges between 0.37 and 0.44 [39]; however, different experiments suggest that, in the case of ultrathin PET films, it depends considerably on several variables, such as temperature, relative humidity and load, with effective values ranging from 0.29 to 0.37 [40, 41].

In order to investigate the effect of both Young’s modulus and Poisson’s ratio of the film on the contact area, a preliminary computational experiment was designed using the Optimal Latin Hypercube technique as implemented in the design exploration software Isight [42]. The algorithm generated 100 optimally distributed combinations of E and ν ranging approximately from 0 to 5000 MPa and from 0 to 0.5, respectively, resulting in a densely populated factor space. The parameter combinations were successively introduced in the axisymmetric finite element model described above. This was then simulated with the implicit dynamic solver Abaqus/Standard using the Backward Euler integrator, which introduces high numerical damping that provides the energy dissipation necessary for quasi-static applications. As for geometry and load, the two configurations corresponding to the bounds of the design space were considered, namely, the R300 and R600 bodies with normal forces of 25 kN and 100 kN, respectively.

The results of the preliminary study, excluding few failed runs due to an excessively low elastic modulus, are illustrated in Fig. 2. It became evident that modifying ν does not cause significant variation in the contact area registered by the film model excepting a slight alteration near the incompressibility limit. For this reason, its value was fixed to 0.33, which is consistent with own previous models and falls in the range of experimental observations for thin PET films [40, 41].

Fig. 2
figure 2

Effect of Young’s modulus and Poisson’s ratio of a linear-elastic Prescale film model on the simulated sphere-plane contact area for R300 and 25 kN (top) and R600 and 100 kN (bottom)

Full size image

The full study was thus limited to determining whether a value of E can be assigned to any of the film types so that the measurement errors collected in Table 1 are reproduced with high accuracy for all scenarios. To this end, the model described above was tuned to each of the ten combinations of geometry and load that were considered in the experiments. Then, each configuration was subjected to several runs with varying E in order to find the values that lead to a simulated contact area equal to the experimentally determined values. The elastic modulus ranges that fulfilled said requirement are depicted in Fig. 3 classified by film type.

Fig. 3
figure 3

Film stiffness ranges that produce simulated contact areas equal to experimental measurements for different body curvatures and normal loads

Full size image

Regardless of the chosen film, the elastic modulus necessary to reproduce the measured contact area at the lowest load level differed significantly depending on body curvature. For the remaining load levels, however, the experimental values could be replicated with similar film stiffness values for all curvatures. The best fit was achieved for the HS film. In fact, a Young’s modulus of 2800 MPa (dashed line in Fig. 3) produced virtually equal results as compared to experimental outcomes.

To sum up, a linear-elastic film model with material parameters E = 2800 MPa and ν = 0.33 and discretized as a single layer of 100 × 90 μm gasket elements was capable of reproducing exactly the experimental contact area measurements obtained with HS Prescale film under normal loads ranging from 50 to 100 kN. Owing to the purely phenomenological calibration, the application range of the resulting film model is limited to the assessment of contact areas within the parameter space where experimental measurements could be reproduced. However, the load limitation is not particularly restrictive as the lower and upper limits approximately stand for the representative scenarios of half and full axle load.

2.2 Wheel-rail contact

The calibrated HS film model was finally included into a three-dimensional wheel-rail contact model with profiles S1002 [3] and 60E2 [4] as outlined in Fig. 4.

Fig. 4
figure 4

Wheel-rail contact model with interposed film on one side of the symmetry plane

Full size image

On the one hand, the wheel profile extended from the flange tip to the end of the segment with continuously variable curvature. The wheel cutout showed a nominal radius of 460 mm at the tape circle section, which was contained in the xz-plane, and was bounded by a concentric circular surface 12 mm over the tread datum. A circular segment of 2.5º was considered in the circumferential direction, which was effectively doubled by defining a symmetry boundary condition with respect to the yz-plane as shown in Fig. 4. On the other hand, the rail profile extended from the inner gauge corner to the end of the running surface. The rail cutout was bounded by corresponding horizontal and vertical surfaces and extruded to a length of 20 mm in the x-direction, which was also effectively doubled by the symmetry boundary condition. In addition, a rail inclination of 1/40 was considered. The film was modeled with a thickness of 90 μm and conforming to the rail surface in order to facilitate convergence.

The idealization error that originates in wheel-rail contact area simulations when the contact partners are reduced to the described cutouts has already been investigated and is limited to a slight deviation concerning size and position of the contact patch [22]. As was done in said research, both contact bodies were meshed with reduced-integration, hexahedral elements with an approximate element size of 0.75 mm. Following the calibration presented above, the film was discretized as a single layer of gasket elements. However, the mesh was coarser than in the calibrated model as a convergence study confirmed that the element size in the plane of the film can be increased to 0.25 mm for the case under study without significant accuracy loss.

In order to take into consideration the effect of the surrounding material, the front and lateral surfaces of both cutouts were constrained in the x- and y-directions, respectively. Moreover, the lower surface of the rail cutout was fixed in both the y- and z-directions and the upper surface of the wheel part was kinematically connected to a force application point in the yz-plane that was only allowed to undergo vertical translational motion. Contact settings in the normal and tangential directions matched the calibration model at both film-wheel and film-rail interfaces.

The model was simulated for normal loads of 50 kN and 100 kN considering different relative positions of wheel and rail in the y-direction. As depicted in Fig. 5, the transverse offset η was defined as the wheel lateral shift with respect to the position where the tape circle and the longitudinal symmetry plane of the rail intersect at the contact point, and was defined positive when the offset makes the wheel flange approach the rail corner reducing the gap between the profiles. For each of the two load levels, several runs were made starting with a transverse offset of η = -7 mm and increasing it by 1 mm up to η = 2 mm. All runs were then repeated after removing the film from the model and defining hard contact between both steel bodies.

Fig. 5
figure 5

Convention for transverse offset

Full size image

Additionally, the impact of roughness on the measured contact area was simulated. As already demonstrated, the increase in nominal contact area due to roughness, which is indeed detected by HS films, can be reproduced with great precision within the load range considered here if an appropriate interaction characteristic is defined [27]. The hard contact interaction between film and wheel was replaced accordingly by an experimentally obtained pressure-overclosure law that incorporates microscopic asperity interaction between film and roughened steel with an approximate roughness arithmetic average of Ra ≈ 1.5 µm.

All simulations were performed with the implicit dynamic solver Abaqus/Standard using the Backward Euler integrator.

3 Results and discussion

The simulated size of the contact area at the wheel-film interface is plotted in Fig. 6 as a function of normal force and transverse offset for the three analyzed models: without film, with film and with film plus rough wheel surface. These diagrams can be used directly to filter out the systematic deviation caused by the presence of Prescale film of HS type between wheel and rail under smooth surface condition.

Fig. 6
figure 6

Simulated contact area A depending on transverse offset η and normal force F

Full size image

The shape of the contact area varied in a similar way regardless of the load. The resulting shapes for smooth contact with and without film are outlined in Fig. 7 in terms of the transverse offset for a load of 50 kN as an example.

Fig. 7
figure 7

Shape of simulated contact area without (solid line) and with film (dashed line) depending on transverse offset η for a normal force of F = 50 kN

Full size image

In view of the results, it is manifest that Prescale films do not alter significantly the shape of the wheel-rail contact area, allowing thus a straightforward qualitative assessment, but do introduce a non-negligible error when used to determine its size. As expected based on the values already available for related Hertzian geometries and collected in Table 1, little influence of load and transverse offset (thus contact geometry) on measurement error was observed. Therefore, even constant correction values are conceivable depending on the accuracy level required.

However, even a small roughness increase, in this case from smooth surface condition to service roughness of one contact partner, can alter considerably the size of the contact area imprinted on the film as illustrated in Fig. 6. Although much smaller in service, the rail roughness after grinding can amount to Ra = 10 μm [25] and a freshly reprofiled wheel may show values as high as Ra = 30 μm [43]. Taking this into consideration, it is a fundamental requirement to perform a profile roughness analysis of the contact bodies prior to error correction in order to determine whether the curves deduced from smooth contact simulations can yield an acceptable approximation of the actual measurement error. If not, the simulation process presented above should be adapted by implementing the roughness range of interest for both with-film and without-film scenarios and at both contact interfaces if necessary. As a matter of fact, the parametric model calibration proposed here is based on experimental data, so it determines only the fitting ability of the model but not its predictive capability outside the analyzed parameter range [16]. In general, the film model should always be calibrated under consideration of its intended use.

4 Conclusions

The study described above demonstrated, first of all, that it is possible to calibrate a linear-elastic model of HS Prescale film with a single parameter set so that it replicates experimental contact area measurements for dissimilar Hertzian geometries. This is a necessary condition for including it in contact models involving contact partners with variable curvature in different directions as is the case in wheel-rail contact.

Finally, it can be gathered from the results that including the calibrated film model in finite element simulations, as well as characterizing surface roughness properly, are necessary conditions for model validation by comparison with experimental results. Not doing so results in a poor representation of the reality of interest and hence in non-negligible idealization error [44], which lowers the significance of a direct comparison between simulation and experimental outcomes for validation purposes.