Metamodelling of wheel–rail normal contact in railway crossings with elasto-plastic material behaviour

Abstract

A metamodel considering material plasticity is presented for computationally efficient prediction of wheel–rail normal contact in railway switches and crossings (S&C). The metamodel is inspired by the contact theory of Hertz, and for a given material, it computes the size of the contact patch and the maximum contact pressure as a function of the normal force and the local curvatures of the bodies in contact. The model is calibrated based on finite element (FE) simulations with an elasto-plastic material model and is demonstrated for rail steel grade R350HT. The error of simplifying the contact geometry is discussed and quantified. For a moderate difference in contact curvature between wheel and rail, the metamodel is able to accurately predict the size of the contact patch and the maximum contact pressure. The accuracy is worse when there is a small difference in contact curvature, where the influence of variation in curvature within the contact patch becomes more significant. However, it is shown that such conditions lead to contact stresses that contribute less to accumulated plastic deformation. The metamodel allows for a vast reduction of computational effort compared to the original FE model as it is given in analytical form.

1 Introduction

Railway turnouts (switches & crossings, S&C) are subjected to a severe load environment, in particular in the switch and crossing panels, see, e.g., [8, 18, 23]. High maintenance costs are generated because of the needs for repair and replacement of switch rails and crossings due to damage in the forms of wear, plastic deformation and breaking out of material caused by surface or subsurface initiated rolling contact fatigue cracking.

A wheel travelling over a crossing generally induces an impact load on the crossing rail (for traffic in the facing move: from switch panel to crossing panel) or on the wing rail (for traffic in the trailing move: from crossing panel to switch panel), see Fig. 1. The impact load is caused by the downwards–upwards motion experienced by the wheel as it rolls along the wing rail and over to the crossing rail or vice versa. On the wing rail, the vertical motion is caused by the conicity of the wheel and the significant change in lateral wheel–rail contact position that occurs due to the lateral deviation of the wing rail in the running direction. On the crossing rail, the lateral contact position is relatively constant but the crossing rail has a vertical inclination.

Fig. 1

Components of a turnout (from [10])

Field tests have been reported, where impact loads at crossings have been measured by instrumented wheelsets [9, 20] or by rail-mounted accelerometers [13]. Depending on vehicle speed and wheel–rail contact geometry, the impact load may be considerably higher than the nominal static wheel load [1]. Increasing vehicle speeds and axle loads, and wheels with severely worn profiles, induce contact conditions that generate higher magnitudes of wheel–rail contact forces and slip. Accelerated rail profile degradation and damage occur if the rail profiles are not corrected in time, since the deteriorated rail profiles induce contact conditions that further magnify the dynamic loads.

Optimisation of rail profiles [19, 31] and support stiffness [3] are two means to reduce the wheel–rail contact forces and damage in the crossing panel. Another approach is to use a higher strength material for the crossing. Common materials used in crossings include austenitic steels (such as the explosively hardened manganese steel Mn13), bainitic steels (B360) and high strength pearlitic rail steels (such as the fine-pearlitic rail grade R350HT). Based on field experience, each of the materials has different advantages in terms of resistance to the various types of damage. For example, the manganese steel has a higher rate of plastic deformation and hardening during the initial load cycles but at the same time a better adaptability of rail profile to meet the variation of worn wheel profiles in traffic. In addition, due to its work-hardening, the manganese steel has good toughness and wear resistance under severe impact loading conditions (see, e.g., [4]).

A cross-disciplinary simulation methodology, integrating several numerical tools, for predicting the damage in rail profiles of S&C by an iterative approach has been proposed by Johansson et al. [5]. In each iteration step, the methodology consists of four parts: (I) simulation of dynamic vehicle–track interaction to predict wheel–rail contact forces, creepages and contact positions; (II) analysis of wheel–rail contacts, using a three-dimensional finite element (FE) model with an elasto-plastic material model, to determine the size of the contact patches and stresses in the material due to the normal contact; (III) prediction of accumulated damage including plastic deformation, wear, and rolling contact fatigue; and (IV) updating of rail profiles, which are then used as input in the next iteration step. The methodology has been validated by comparing predicted rail profiles with those measured in the field at Härad in Sweden [5] and Haste in Germany [16]. A similar method based on an explicit FE model has been presented in [34].

The long-term objective of the described methodology is to enable efficient and accurate predictions of the influence of material selection on long-term S&C rail profile degradation to design better S&C. However, to represent the variation in traffic that a given crossing is subjected to, a distribution of load parameters such as wheel profile, train speed, axle load and wheel–rail friction coefficient needs to be accounted for. This in turn leads to a wide range of contact load scenarios with different magnitudes of wheel–rail contact load occurring at various positions along the crossing panel. Each of these scenarios needs to be accounted for by an FE simulation in part II of the methodology, which, due to the required combination of high mesh resolution in three dimensions and non-linear material model, is computationally expensive. Thus, there is a need to improve the methodology in terms of computational efficiency. Substituting all of the required FE simulations in each iteration step (part II above) with low-cost calculations using a calibrated analytical metamodel, with the distribution of output data calculated in step I as input, would reduce the computational cost significantly. A procedure to generate such a metamodel with limited loss in accuracy compared to the full elasto-plastic FE model is described in the present paper.

This paper is organized in the following way: Sect. 2 describes tools available for the modelling of wheel–rail contact; the FE model used in this study is presented in Sect. 3; two alternative metamodels that are used to replace the FE simulations are described and compared in Sect. 4; Sect. 5 investigates the influence of constant curvature assumed for the contact geometry on the accuracy of the metamodels; Sect. 6 presents a demonstration of the suggested approach, and finally, concluding remarks are stated in Sect. 7.

The following notation convention has been adopted in this paper: small Latin characters in bold (e.g., $\varvec{a}$) denote vectors, while capital characters (e.g., $\varvec{A}$) or small Greek characters (e.g., $\varvec{\sigma }$) in bold denote tensors of the second order, unless stated differently by the context.

2 Wheel–rail contact in the crossing panel

One key aspect in the simulation of rail profile degradation is to determine the correct contact pressure distribution in the wheel–rail contact. An overview of available contact models can be found in the literature (see, e.g., [26]). Several tools exist for solving the normal contact problem, for example, solutions based on (1) the theory of Hertz [12]; (2) Kalker’s variational method [7]; and (3) FE simulation. The first two approaches are applicable for elastic material response only. In addition, they rely on the assumption that the bodies in contact are large compared with the dimensions of the contact patch, such that they can be considered as infinite half-spaces. The same holds for most fast methods for vehicle dynamics simulations (see, e.g., [25]). Furthermore, the theory of Hertz assumes that the geometry of each contact surface can be approximated by an elliptic paraboloid. Kalker’s method imposes no restriction on the wheel–rail contact geometry, while the FE simulation approach allows for both arbitrary geometry and inelastic material response.

Typically, the size of the contact patch is in the order of 1 cm${}^{2}$, which is small in comparison with the dimensions of the wheel and rail [11]. Furthermore, the curvatures of wheel and rail are constant over multiple regions of the profiles, see Fig. 2. However, there are cases when at least one of the assumptions in Hertzian theory is violated. In particular, such situations can occur (especially with worn profiles) when the contact is at the gauge corner of the rail or at the flange of the wheel. For such situations, it has been reported [11] that the results of the Hertzian solution are in vast disagreement with FE simulations. In the presence of plastic deformation, the FE method (FEM) is the only solution that is available, although at a large computational cost. The local coordinate systems used to describe the wheel and rail profiles, the contact positions on wheel and rail, and the relative lateral wheel–rail displacement $\varDelta y$ are marked in Fig. 3. Contact points positions on the wheel and the rail cross section are plotted in Fig. 4 for different $\varDelta y$.

Fig. 2

Calculated curvatures within the potential contact area

Fig. 3

Example of wheel and rail profiles in the crossing panel. Local coordinate systems $(y_{\text {rail}}, z_{\text {rail}})$ and $(y_{\text {wheel}}, z_{\text {wheel}})$ with origins at the centre of the rail and at the nominal running circle, respectively. Relative wheel–rail displacement $\varDelta y$ (positive direction when wheel flange is approaching the rail gauge corner). Contact positions on wheel $y_{\text {wheel}}^{\text {contact}}\Bigr |_{\varDelta y}$ and rail $y_{\text {rail}}^{\text {contact}}\Bigr |_{\varDelta y}$

Fig. 4

Positions of the contact points on the wheel and rail for different relative lateral wheel–rail displacements $\varDelta y$

A comparison between calculated contact conditions for elastic material when applying an in-house implementation of Kalker’s variational method [21] or an FE model, using either the nominal geometry (accounting for the variation in curvature within the contact patch) or the simplified geometry (constant radius of curvature), will be presented in Sect. 5. It will be shown that there is excellent agreement between the two methods when the nominal geometry is used. However, the aim of the procedure presented in this paper is to account for plasticity due to wheel–rail contact in a parameterised and numerically efficient way. If plasticity is overlooked and the contact pressure is computed based on the elastic material response, the stress in the rail will, at some contact locations, be so high that it will (in part III of the simulation methodology) lead to a drastic overestimation of the accumulated plastic strain or to an absence of the equilibrium solution when the estimated stress exceeds the ultimate stress. Therefore, a plasticity model together with a simplified geometry of the wheel–rail contact will be proposed in this paper.

3 FE model of wheel–rail contact

3.1 FE mesh

The FE model of the wheel–rail contact represents two cut-out pieces of wheel and rail that are pressed against each other in the normal direction. To allow for a straightforward parameterisation of the contact geometry, the contact surfaces of the two profiles are approximated based on the assumption of constant curvature within the contact patch (same as in the theory of Hertz). Figure  5 illustrates how a circle is fitted to each profile to determine the wheel $R_{x}^{\text {w}}$ and rail $R_{x}^{\text {r}}$ radii at the point of contact ($yz$ plane). The circle is fitted on the region that spans 1 mm along the $y$-axis on either side of the contact point, which was found sufficient to assure that the radii that define a nominal 60E1 profile are recovered. The lateral positions of the contact points are denoted $y_{\text {wheel}}^{\text {contact}}\Bigr |_{\Delta y}$ and $y_{\text {rail}}^{\text {contact}}\Bigr |_{\Delta y}$, see Fig. 3. Three rail cross sections labelled A, B, and C (located at 754, 1008 and 1166 mm from the theoretical crossing point in a 60E1-R760-1:15 turnout, respectively), and a nominal S1002 wheel profile is used to demonstrate the procedure in this study. For a given wheel–rail pair, the locations of the contact points on the wheel and rail are determined by a simulation in the KPF (contact point function) module of GENSYS [2], which is a commercial software for rail vehicle dynamics (see also the methodology described in [5]). Only a quarter of each cut-out piece needs to be simulated due to the assumption of symmetry. The FE model is generated by means of Python scripts through Abaqus Scripting Interface [29]. An example of how to create 3D contact jobs in Abaqus via its API (Application Programming Interface) is described in [30].

Fig. 5

Approximation of wheel and rail profile geometry for the nominal ($\varDelta y=0$) position of the S1002 wheel with respect to three rail cross sections in the crossing panel. Circles are fitted to determine wheel radius $R_{x}^{\text {w}}$ and rail radius $R_{x}^{\text {r}}$

The modelled quarter of the cut-out piece of the rail is shown in Fig. 6. The upper side of the block, which is pressed against the wheel, has constant curvature $1 / R_{x}^{\text {r}}$ in the $yz$ plane and zero curvature in $xz$ plane (the rail has a constant inclination in the running direction). Its bottom side is fixed, while the two front sides in the $xz$ and $yz$ planes are the planes of symmetry, i.e. the out-of-plane deformation is set to zero. The back sides in the $xz$ and $yz$ planes are set to be free boundaries, because the identification of the stiffness of linear springs that would correspond to physical boundary conditions was considered superfluous, since Fig. 7 shows negligible effect stemming from the free boundaries. The edge length $L$ was chosen to assure that the area $A$ of the contact patch and the maximum contact pressure $p_{0}$ for the elastic wheel and rail are in close agreement with those provided by the contact theory of Hertz. In the following, the FE model with constant curvature illustrated in Fig. 6 will be referred to as “simplified FE model”.

Fig. 6

FE mesh of the modelled piece of rail. Eight of the edges have length $L$. The remaining four edges are different from $L$ due to rail curvature $R_{x}^{\text {r}}$

The influence of the edge length $L$ on the key responses is illustrated in Fig. 7 for the case of nominal ($\varDelta y=0$) wheel position with respect to rail section C and 22.5 tonnes axle load. The responses, obtained with FEM using the penalty method (achieved by using the “surface-to-surface” type of contact pair and the “hard” pressure–overclosure relationship in Abaqus, see [28] ), are normalised by the responses from the Hertzian solution. In these simulations, the wheel and the rail are both assumed to be linear elastic with Young’s modulus $E = 183$ GPa and Poisson’s ratio $\nu = 0.3$. It is observed that the responses deviate by less than 5% from the corresponding Hertzian solutions for all studied edge lengths. The reason for this small deviation is the spatial discretisation of the FE model as the accuracy of measuring the length of the semi-axes is limited by the edge length of the elements in the contact patch. Since it had been found (not shown here) that the length of one of the semi-axes can be greater than 20 mm for certain contact scenarios, the edge length of 30 mm has been adopted in this study. An analogous consideration was undertaken when choosing the mesh size. Figure 8 shows how the responses vary relative to the Hertzian solution for different element sizes in the potential contact zone. Furthermore, two mesh densities were considered on the outer boundary of the model. The 15 and 5 mm element edge lengths $l$ used in the outer mesh are referred to as “coarse mesh” and “fine mesh”, respectively. The difference in computational effort with respect to the finest mesh resolution (0.2 mm in the contact region and 5 mm in the outer region) is illustrated in Fig. 9. To conclude, the mesh adopted in the following of this study has tetrahedral elements with 0.25 mm edge length in the potential contact zone (fine inner mesh) and 15 mm on the outer boundary (coarse outer mesh).

Fig. 7

Influence of block size (edge length $L$) of the FE model on calculated area of contact patch $A$ and maximum contact pressure $p_{0}$ in relation to the Hertzian solution. The values of semi-axes $a$ and $b$ vary around 7.2 and 7.6 mm, respectively

Fig. 8

h-convergence study for FE mesh. Normalisation with respect to Hertzian solution. Solid lines denote coarse outer mesh, dashed–dotted lines—fine outer mesh

Fig. 9

Computational time gain for different discretisation levels of the FE mesh

3.2 Material modelling of R350HT

Since accumulated plastic deformation has been identified as an important damage mechanism in railway turnout components [5], an appropriate material model needs to be utilised. Switch and crossing components may be manufactured out of different materials such as manganese steel, bainitic steel and pearlitic steel. In this paper, the pearlitic steel R350HT is chosen as the prototype material. Cyclic uniaxial tests of R350HT have been conducted by Schilke [24]. The data from these tests have been used to calibrate a cyclic plasticity model. The adopted cyclic plasticity model was proposed by Ohno and Wang, see [17], and it is summarized in Appendix A. The model parameters collected in the vector $\varvec{p}_{\text {r}}$:

$$\begin{aligned} \varvec{p}_{\text {r}} = \left[ E \quad \nu \quad \sigma _{y} \quad C_{1} \quad \gamma _{1} \quad m_{1} \quad C_{2} \quad \gamma _{2} \quad m_{2} \right] \end{aligned}$$

(1)

have been identified from a least-square’s fit against the cyclic uniaxial tests, and the numerical values are given in Table 1. The model response and the test results are compared in Fig. 10. Good agreement between the simulated and measured material response is observed for two different load conditions. To achieve this agreement, it was necessary to use two kinematic hardening stresses in the Ohno–Wang model. For the wheel, linear isotropic elasticity is assumed with $E_{\text {w}} = 183$ GPa and $\nu _{\text {w}} = 0.3$. These parameters define the vector:

$$\begin{aligned} \varvec{p}_{\text {w}} = \left[ E_{\text {w}} \quad \nu _{\text {w}} \right] . \end{aligned}$$

(2)

The effect of assuming elastic material behaviour for the wheel is examined in Appendix C.

Table 1 Identified material parameter values (Poisson’s ratio $\nu$ is not calibrated) for the Ohno–Wang model for rail grade R350HT

Fig. 10

Experiments and simulations with cyclic plasticity model for rail grade R350HT

4 Metamodelling

A metamodel (surrogate model) can be described as a model of a model, i.e. an approximation of another model, and metamodelling is the process of generating such metamodels. Metamodels find vast applications in different fields, such as engineering design optimisation, software, and systems engineering, and are mostly being employed to reduce the computational cost. Reviews of different approaches used in metamodelling can be found in the literature, see e.g. [27, 32]. It is concluded that metamodelling involves three basic steps: (1) choosing an experimental design (sampling); (2) choosing the type of metamodel; and (3) fitting the metamodel to the data (calibration). Prior to that, the interface to the metamodel needs to be identified, i.e., the input and output variables. The interface to the metamodel developed in the present study is summarized as

$$\begin{aligned} \hat{\varvec{g}} = \hat{\varvec{g}} \left( R_{x}^{\text {w}}, R_{y}^{\text {w}}, R_{x}^{\text {r}}, F_{n}; \varvec{\beta } \left( \varvec{p}_{\text {w}}, \varvec{p}_{\text {r}} \right) \right) \end{aligned}$$

(3)

where the output variables $\hat{\varvec{g}}$ are the semi-axes $\hat{a}$ and $\hat{b}$ of the contact patch and the maximum contact pressure $\hat{p}_{0}$, i.e. $\hat{\varvec{g}} = \left[ \hat{a}, \hat{b}, \hat{p_{0}} \right] ^{T}$. Furthermore, the input variables are the wheel radius $R_{x}^{\text {w}}$ and rail radius $R_{x}^{\text {r}}$ in the plane of the wheel–rail cross section, wheel radius $R_{y}^{\text {w}}$ in the circumferential direction (see also Sect. 3) and the magnitude of the normal wheel–rail contact force $F_{n}$. The metamodel parameters are collected in the vector $\varvec{\beta }$ which are determined from calibration. Since these are calibrated for a certain combination of wheel and rail material, they are functions of $\varvec{p}_{\text {w}}$ and $\varvec{p}_{\text {r}}$.

4.1 Sampling

The simulation data that need to be sampled are the results of the FE simulations in terms of maximum contact pressure $p_{0i}^\mathrm{s}$ and semi-axes $a_i^\mathrm{s},\, b_i^\mathrm{s}$ (superscript s for sampling) of the wheel–rail contact described in Sect. 3. These data will be used to calibrate the metamodel. The number of samples per output is assumed to be M, i.e. $i=1,\ldots ,M$. For a given pair of wheel and rail profiles, and for a given relative lateral displacement $\varDelta y$ between wheel and rail (recall Fig. 3), the contact positions on wheel and rail are solved in the KPF module of GENSYS. The range of relative lateral wheel–rail displacements $\varDelta y$ in the crossing panel is constrained in one direction by the contact between the wheel flange and the gauge side of the crossing rail, and in the opposite direction by the contact between the back of the wheel flange (on the other wheel in the wheelset) and the check rail. The range of relative wheel–rail displacements considered in this study is set to $\varDelta y \in \left[ -5,5 \right]$ mm. Note that a positive value of $\varDelta y$ means that the wheel flange has moved from the nominal wheel–rail contact position towards the rail. For each predicted contact position on wheel and rail, the curvature is calculated by fitting two circles (see Sect. 3) and thereby obtaining $R_{x}^{\text {w}}, R_{x}^{\text {r}}$. The other wheel radius $R_{y}^{\text {w}}$ is calculated by knowing the location of the contact point and the wheel radius at the running circle.

Based on the combination of a nominal S1002 wheel profile and rail cross section C, the influence of relative lateral wheel–rail displacement $\varDelta y$ on the radii of curvature on wheel and rail is illustrated in Fig. 11. It is observed that the magnitude of the wheel radius is decreasing with increasing values of $\varDelta y$ as the contact position is approaching the flange root. For $\varDelta y \in \left[ -1,3.5 \right]$ mm, the rail radius is in the order of 200 mm as contact is taking place close to the centre of the crossing rail. Note that the values of the wheel radius $R_{x}^{\text {w}}$ about the $x$-axis are negative. This means the wheel tread has a concave surface that together with the opposite sign of the rail radius $R_{x}^{\text {r}}$ leads to a certain degree of conformity within the selected range of $\varDelta y$.

Fig. 11

Relation between wheel/rail radii $R_{x}^{\text {w}}$ and $R_{x}^{\text {r}}$ (negative sign denotes a concave surface) and relative lateral wheel–rail displacement $\varDelta y$

Conventional sampling techniques, e.g., factorial and central composite designs, see, e.g., [27], assume that all combinations of the variables within the design space are possible. Consider the case denoted by the square in Fig. 12. This could correspond to one of the sampling points produced by such a sampling technique, since it is a combination of the minimum value of $R_{x}^{\text {w}}$ and the maximum value of $R_{x}^{\text {r}}$. However, as can be concluded from Fig. 11, the corresponding scenario cannot occur. Hence, such a sampling point should be avoided. Instead, a hand-picked technique based on the contact geometry and pre-calculations of contact positions in GENSYS is applied. For the chosen range of $\varDelta y$ with step length 0.5 mm, all integer values of $\varDelta y$ are selected as sampling points, while all the fractional values are used as test points to compare the results from the calibrated metamodel and the FE simulations. At this point, it is worth mentioning that not all the values of $\varDelta y$ lead to permissible contact scenarios due to the local constant curvature approximation. For instance, the value of $\varDelta y = 4$ mm corresponds to $|R_{x}^{\text {w}} |< R_{x}^{\text {r}}$, which leads to an impossible contact scenario. Such cases are ignored in the sampling. In addition, note that the radius of curvature of the rail profile along the running direction is constant, i.e., $R_{y}^{\text {r}}=\infty$.

Fig. 12

Taxonomy of contact scenarios for $\varDelta y \in [-5, 5]$ mm

4.2 Metamodels

4.2.1 Polynomial metamodel

One particular metamodelling approach that is very common and does not require many sampling data points is Response Surface Methodology (RSM), see, e.g., [14]. In RSM, the model is a polynomial of a certain order. The selected order of the polynomial depends on the anticipated curvature of the response surface within the given design space. Typically, a first-order polynomial is used for flat response surfaces, while a second-order polynomial is employed when the curvature is strong. Here, the quadratic metamodel is written as (superscript Q for quadratic)

$$\begin{aligned} \hat{\varvec{g}}^{\text {Q}}&= \varvec{\beta }_{0}^{\text {Q}} + \varvec{\beta }_{1}^{\text {Q}} R_{x}^{\text {w}} + \varvec{\beta }_{2}^{\text {Q}} R_{y}^{\text {w}} + \varvec{\beta }_{3}^{\text {Q}} R_{x}^{\text {r}} + \varvec{\beta }_{4}^{\text {Q}} F_{n} \nonumber \\&\quad + \varvec{\beta }_{5}^{\text {Q}} \left( R_{x}^{\text {w}}\right) ^{2} + \varvec{\beta }_{6}^{\text {Q}} \left( R_{y}^{\text {w}}\right) ^{2} + \varvec{\beta }_{7}^{\text {Q}} \left( R_{x}^{\text {r}}\right) ^{2} + \varvec{\beta }_{8}^{\text {Q}} \left( F_{n}\right) ^{2} \nonumber \\&\quad + \varvec{\beta }_{9}^{\text {Q}} R_{x}^{\text {w}} R_{y}^{\text {w}} + \varvec{\beta }_{10}^{\text {Q}} R_{x}^{\text {w}} R_{x}^{\text {r}} + \varvec{\beta }_{11}^{\text {Q}} R_{x}^{\text {w}} F_{n} \\&\quad + \varvec{\beta }_{12}^{\text {Q}} R_{y}^{\text {w}} R_{x}^{\text {r}} + \varvec{\beta }_{13}^{\text {Q}} R_{y}^{\text {w}} F_{n} + \varvec{\beta }_{14}^{\text {Q}} R_{x}^{\text {r}} F_{n} \nonumber \end{aligned}$$

(4)

with $\hat{\varvec{g}}^{\text {Q}}$ defined as $\hat{\varvec{g}}^{\text {Q}} = \left[ \hat{a}^{\text {Q}} \quad \hat{b}^{\text {Q}} \quad \hat{p}_{0}^{\text {Q}} \right] ^{T}$ and $\varvec{\beta }_{i}^{\text {Q}} = \left[ \beta _{ia}^{\text {Q}} \quad \beta _{ib}^{\text {Q}} \quad \beta _{ip_{0}}^{\text {Q}} \right] ^{T}$ for $i=0,\ldots ,14$. The $\hat{(\cdot )}$ denotes an estimate of $(\cdot )$. For instance, $\hat{a}$ is the semi-axis within the elliptic contact patch along the $x$-axis (in the direction of travel) given by the metamodel; its counterpart $a$ is the reference value of the semi-axis obtained by the experiment (FE simulation).

4.2.2 Hertzian-based metamodel

In this paper, a Hertzian-based metamodel is proposed with the main assumption that the geometry of the contact patch is Hertzian, i.e. elliptic. Yet the size of the contact patch might differ from the Hertz’s solution since plasticity of the rail is accounted for. According to the theory of Hertz, the ratio $\alpha = a / b$ of the semi-axes of the contact patch ellipse can be calculated, see [12], by solving the equation:

$$\begin{aligned} f(\alpha ;R_{x}^{\text {w}}, R_{y}^{\text {w}}, R_{x}^{\text {r}} )&= \left| \left( \frac{1}{R_{x}^{\text {w}}} - \frac{1}{R_{y}^{\text {w}}} \right) +\frac{1}{R_{x}^{\text {r}}} \right| \nonumber \\&\quad -\frac{(\alpha ^{2}+1)E_{2} - 2 E_{1}}{(\alpha ^{2} - 1)E_{2}} \sum \frac{1}{R} = 0 \end{aligned}$$

(5)

where $E_{1}$ and $E_{2}$ are the complete elliptic integrals of the first and second kinds which are functions of $\alpha$ via the modulus $k^{2} = 1 - 1/\alpha ^{2}$, and $\sum \frac{1}{R}= \frac{1}{R_{x}^{\text {w}}} + \frac{1}{R_{y}^{\text {w}}} + \frac{1}{R_{x}^{\text {r}}}$. It is worth noting that Eq. (5) only contains geometrical variables. Hence, the following metamodel is proposed (superscript H for Hertz):

$$\begin{aligned} \alpha&= {\text {arg}} \left( f \left( \alpha ; R_{x}^{\text {w}}, R_{y}^{\text {w}}, R_{x}^{\text {r}} \right) = 0 \right) \end{aligned}$$

(6)

$$\begin{aligned} \hat{a}&= \beta _{0a}^{\text {H}} \left[ \frac{\alpha ^{2} E_{2} F_{n} }{\sum \frac{1}{R}} \right] ^{\beta _{1a}^{\text {H}}} \end{aligned}$$

(7)

$$\begin{aligned} \hat{b}&= \frac{\hat{a}}{\alpha } \end{aligned}$$

(8)

$$\begin{aligned} \hat{p}_{0}&= \beta _{0p_{0}}^{\text {H}} \left[ \frac{F_{n}}{\pi \hat{a} \hat{b}} \right] ^{\beta _{1p_{0}}^{\text {H}}}. \end{aligned}$$

(9)

The expressions for the semi-axis $\hat{a}$ and the maximum normal contact pressure $\hat{p}_0$ are also inspired by the Hertzian solution. Note that for this metamodel the number of parameters is only 4, which can be compared with the quadratic metamodel in Eq. (4) with 45 parameters.

4.3 Calibration of metamodels

The metamodel parameters $\varvec{\beta }_{i}$ are identified using a non-linear regression analysis. The analysis is non-linear because for the Hertzian case the model $\hat{\varvec{g}}$ is a non-linear function of the parameters $\varvec{\beta }_{i}$. The parameters are identified by defining an objective function h and performing a minimisation of this function. The objective function is chosen in a staggered fashion with a relative error for each of the responses:

$$\begin{aligned} h = \frac{1}{M} \sum _{e=1}^{M} \left[ \frac{ \hat{g}_{e} - g_{e}^{\text {s}} }{g_{e}^{\text {s}}} \right] ^{2} \end{aligned}$$

(10)

where $M$ is the number of sampling results (FE simulations) and superscript s stands for sampling. The minimisation of h for the polynomial metamodel is performed using the “Powell” algorithm, see [22], and for the Hertzian-based metamodel using the “Nelder-Mead” algorithm, see [15].

The metamodels were implemented in Python with an extensive use of algorithms and routines for scientific computing provided by the SciPy library [6].

5 Assumption of constant local curvature

For nominal wheel and rail profiles, elastic material response and $\varDelta y = 0$, the effect of the simplification of wheel–rail contact geometry applied in Sect. 3 (often made elsewhere in the literature, see, e.g. [33]) is illustrated in Fig. 13. The maximum contact pressure predicted by Kalker’s variational method is considerably higher when accounting for the variation of curvature within the contact zone. Further, it is observed that Kalker’s solution is not symmetric with respect to $y = 0$. It is also shown that the elastic FE solution based on a model of the nominal surface geometry (extended FE model not accounting for the symmetric boundary conditions) is in excellent agreement with Kalker’s solution.

Fig. 13

Influence of contact geometry on the distribution of normal contact pressure $p$ for $\varDelta y = 0$, $F_{n} = 111.5$ kN and elastic material model. Along the $x$-axis, only half of the plot is shown due to symmetry

For the elasto-plastic material model, a comparison of the maximum contact pressure and maximum von Mises stress when using either the FE model considering the nominal contact geometry or the simplified FE model with constant local curvature is given in Table 2. For some $\varDelta y$, it is observed that with the simplified FE model, there might be an error of up to 33% in maximum pressure and up to 14% in maximum von Mises stress. Furthermore, the influence of $\varDelta y$ on the wheel and rail curvatures of the simplified FE model is summarized in Table 3. By examining Tables 2 and 3, it is found that the accuracy of the simplified models drops significantly for the contacts with smaller (than about 7 m$^{-1}$) difference in the contact curvature.

Table 2 Influence of contact geometry on maximum contact pressure $p_{0}$ and maximum von Mises stress $\sigma ^{\text {Mises}}$ in case of elasto-plastic model for rail grade R350HT

Table 3 Computed rail and wheel curvatures at the points of contact for different wheel–rail lateral displacements

Figure 14 illustrates two cases with either a large or a small difference in contact curvature, which correspond to good or poor agreement with the nominal geometry, respectively. The constant curvatures tuned by the procedure described in Sect. 3 are marked by dots in the figure. The difference between the nominal geometry when there is a variation of curvature within the contact patch (shaded regions in the plots), and the approximated geometry when the curvature is taken as constant and equivalent with the one at the point of contact is shown. It is concluded that when there is a small difference in contact curvature with significant variation in curvature within the contact patch, the accuracy in prediction of maximum contact pressure is reduced. However, in Table 2, it is also observed that the von Mises stress is decreasing with increasing conformity. Since the metamodel is a part of the simulation methodology to predict plastic deformation of S&C components, the cases with larger difference in contact curvature and higher von Mises stress are of main importance. From Table 2 it can also be observed that the maximum von Mises stress is predicted well for the severe load cases, which further strengthens the utility of the simplified FE model as part of the simulation methodology.

Fig. 14

Illustration of a low and b high conformity of wheel and rail profiles at the point of contact. Contact points are marked by dots. Shaded areas denote regions in contact

6 Results

An evaluation of the adopted approach and results of application of the calibrated metamodels are presented in this section. The metamodels were calibrated using the simplified (constant curvature) FE model with the elasto-plastic material model and with the sampling points (integer values of $\varDelta y$) on cross sections A and C. Furthermore, the sampling was performed for four different normal contact loads $F_{n} \in \left\{ 0.5 N, N, 1.5N, 2 N \right\}$, where $N = 111.5$ kN is the nominal static wheel load. The values of the calibrated model parameters are given in Appendix B. The metamodel will be evaluated for section B.

The effect of plastic deformation is shown in Fig. 15, where the pressure distribution within the contact patch calculated with the simplified FE model is plotted together with the distribution when assuming Hertzian theory:

$$\begin{aligned} p (x,y) = p_{0} \sqrt{1 - \left( \frac{x}{a} \right) ^{2} - \left( \frac{y}{b} \right) ^{2}}. \end{aligned}$$

(11)

The distribution denoted as “Hertzian-based metamodel” also uses the Hertzian expression but replaces the parameters $p_{0}, a$ and $b$ with the values obtained by the calibration. As expected, it is observed that the size of the contact patch is increased and the maximum contact pressure is reduced in the presence of plasticity compared with the Hertzian solution that only accounts for the elastic material response.

As described in Sect. 1, the goal of the ongoing research is to study the degradation of rails in S&C. A simplification that has been done is that the wheels are assumed to have been used in service such that they are work-hardened and can be assumed to be linear elastic.

Fig. 15

Elastic vs inelastic material response for $\varDelta y = -5$ mm, $F_{n} = 223$ kN and rail grade R350HT

6.1 Quadratic metamodel

For different positions of contact point on cross section C, a comparison of the responses $a, b$ and $p_{0}$ predicted by the quadratic metamodel and the simplified FE model is shown in Fig. 16. An analogous comparison for cross section B, which the metamodel was not calibrated for, is shown in Fig. 17. The figures reveal a significant spread in contact patch size and maximum contact pressure depending on the lateral position of the point of contact. In addition, it can be seen how non-uniformly the contact points on the rail are located when shifting the wheel by a constant value, since the data points are not uniformly spaced. Furthermore, it is evident that the polynomial metamodel may lose precision for certain contact point locations. This phenomenon can be explained by the fact that the response surface at these locations cannot be accurately approximated by a quadratic function. Another important observation is that the metamodel is able to predict responses for a cross section it was not calibrated for (cross section B), see Fig. 17.

In summary, based on Figs. 16 and 17, it is concluded that the performance of the quadratic metamodel is relatively poor. Table 4 summarizes the error estimates for the quadratic metamodel for the full load range $F_{n} \in \left[ 0.5N, 2N \right]$.

Fig. 16

Calculated a, b semi-axes and c maximum contact pressure based on quadratic metamodel calibrated for integer values of $\varDelta y$ on sections A and C. The evaluation was performed for fractional values of $\varDelta y$ on section C, $F_{n} = 223$ kN

Fig. 17

Calculated a, b semi-axes and c maximum contact pressure based on quadratic metamodel calibrated for integer values of $\varDelta y$ on sections A and C. The evaluation was performed for fractional values of $\varDelta y$ on section B, $F_{n} = 223$ kN

Table 4 Errors in quadratic metamodel

6.2 Hertzian-based metamodel

Analogous results from the Hertzian-based metamodel are given in Figs. 18 and 19 and in Table 5. It can be seen that this metamodel is superior to the quadratic metamodel as it successfully predicts the responses for all considered contact scenarios. The maximum error in the estimated responses for the considered variation in load is up to 6%. On average, the error for $a, b$ and $p_{0}$ is 2%.

Fig. 18

Calculated a, b semi-axes and c maximum contact pressure from the Hertzian-based metamodel calibrated for integer values of $\varDelta y$ on sections A and C. Evaluation was performed for fractional values of $\varDelta y$ on section C, $F_{n} = 223$ kN

Fig. 19

Calculated a, b semi-axes and c maximum contact pressure from the Hertzian-based metamodel calibrated for integer values of $\varDelta y$ on sections A and C. Evaluation was performed for fractional values of $\varDelta y$ on section B, $F_{n} = 223$ kN

Table 5 Errors in Hertzian-based metamodel

6.3 Coefficient of determination

The coefficient of determination $R^{2}$ is used to estimate the accuracy of the metamodels. It is defined as

$$\begin{aligned} R_{g}^{2} = 1 - \frac{ \sum _{i=1}^{M} \left( g_{i} - \hat{g}_{i} \right) ^{2} }{ \sum _{i=1}^{M} \left( g_{i} - \overline{g} \right) ^{2} } \end{aligned}$$

(12)

where $g \in \left\{ a, b, p_{0} \right\}$ and $\overline{g} = \frac{1}{M} \sum _{i=1}^{M} g_{i}$. A value of $R^{2}$ close to 1 indicates a good performance of the metamodel.

Figure 20 shows the plots of estimates of the key responses by the metamodel versus the reference values obtained from the simplified FE models. The closer the markers are to the line, the better is the agreement. It is observed that the Hertzian-based metamodel shows an excellent agreement with the FE simulations. A comparison of the two presented models is summarized in Table 6. The ability of the quadratic metamodel to predict $b$ is poor. Based on Tables 4, 5 and 6, it is concluded that the Hertzian-based metamodel is more accurate than its polynomial counterpart, even though it has much fewer model coefficients.

Fig. 20

Comparison of a, b semi-axes and c maximum contact pressure computed using Hertzian-based metamodel (calibrated for integer values of $\varDelta y$ on sections A and C, denoted by $\hat{()}$) and simplified FE model. Evaluation was performed for fractional values of $\varDelta y$ on section B, $F_{n} = \left[ 55.75, 223 \right]$ kN

Table 6 $R^{2}$ values for different models (closer to one is better)

The major benefit of using the metamodel can be seen when replacing the FE contact simulations in the simulation methodology briefly described in Sect. 1. Each FE simulation requires approximately 5 min on a computer with a 6-core “haswell” CPU and 128 GB of RAM. This leads to a significant computational effort when considering that thousands of wheel–rail contact scenarios need to be simulated in each iteration step.

7 Conclusions

Metamodels for predicting the size of the contact patch and the maximum contact pressure in simulation of wheel–rail contact in railway crossings have been presented. The first type of model is inspired by response surface methodology and is, therefore, represented by a polynomial. The second type is based on the Hertzian contact theory for two bodies in contact. The performance of both types of metamodels has been quantified and it was found that the Hertzian-based metamodel is more accurate.

The metamodel can take into account plastic deformations and, therefore, is an enhancement of the Hertzian solution. The error stemming from the assumption made in the theory of Hertz regarding the geometry of objects in contact was quantified by comparing with the results from an FE model with the true contact geometry. For the cases considered, the discrepancy ranged from 1 to 33% in maximum contact pressure and from 2 to 13% in maximum von Mises stress. The accuracy was found to be poor for the cases with a smaller difference in contact curvature. However, more importantly from the material degradation point of view is that the cases with a larger difference in contact curvature are captured with a higher precision, i.e. smaller errors were observed for cases with higher values of von Mises stress. This enables the use of the simplified geometry in the simulation methodology [5]. In future work, the presented metamodel will be incorporated in the methodology aiming at robust and efficient predictions of long-term material deterioration in S&C components. The metamodel will facilitate attaining the goal of the methodology since a considerable amount of computational effort is saved by replacing the FE contact simulations by a closed-form solution. In this context, a closed-form solution is a solution expressed in terms of input parameters and simple mathematical operations. Furthermore, the simplification of the geometry allows for the straightforward parameterisation needed for creating a metamodel. There is an obvious trade-off between the accuracy of the single contact simulation and the robustness and efficiency of the overall approach. A further investigation is needed to draw conclusions about the impact of the made assumption.

Appendices

Ohno–Wang plasticity model

In the current section, the adopted plasticity model is briefly summarized. For further details, see [17]. The model is formulated in the small strain framework and the strain $\varvec{\epsilon }$ is assumed to be additively decomposed into an elastic $\varvec{\epsilon }^\mathrm{e}$ and a plastic $\varvec{\epsilon }^\mathrm{p}$ part:

$$\begin{aligned} \varvec{\epsilon }=\varvec{\epsilon }^\mathrm{e}+\varvec{\epsilon }^\mathrm{p}. \end{aligned}$$

(13)

Linear isotropic elasticity is assumed whereby the following relations for the volumetric stress $\sigma _\mathrm{vol}=\varvec{I}: \varvec{\sigma }$ and the deviatoric stress $\varvec{\sigma }_\mathrm{dev}=\varvec{\sigma }- \sigma _\mathrm{vol}/3 \, \varvec{I}$ hold in terms of the elastic strain:

$$\begin{aligned} \sigma _\mathrm{vol}=3 \, K_\mathrm{b} \, \epsilon ^\mathrm{e}_\mathrm{vol}, \quad \varvec{\sigma }_\mathrm{dev}=2 \, G \, \varvec{\epsilon }^\mathrm{e}_\mathrm{dev} \end{aligned}$$

(14)

where G is the shear modulus and $K_\mathrm{b}$ is the bulk modulus. The von Mises yield function $\varPhi$ is adopted which is defined as

$$\begin{aligned} \varPhi =\sqrt{\frac{3}{2} (\varvec{\sigma }_\mathrm{dev} - \varvec{X}):(\varvec{\sigma }_\mathrm{dev} - \varvec{X})}-\sigma _\mathrm{y}=\sqrt{\frac{3}{2}} | \varvec{\sigma }_\mathrm{dev} - \varvec{X} |- \sigma _\mathrm{y} \end{aligned}$$

(15)

where $\varvec{X}$ is the kinematic hardening (back-stress) and $\sigma _\mathrm{y}$ is the initial yield stress. The evolution of plastic strain is assumed to be of associative type:

$$\begin{aligned} \dot{\varvec{\epsilon }}^\mathrm{p}=\dot{\lambda } \, \frac{\partial f}{\partial \varvec{\sigma }}= \dot{\lambda } \, \sqrt{\frac{3}{2}} \frac{\varvec{\sigma }_\mathrm{dev} - \varvec{X}}{| \varvec{\sigma }_\mathrm{dev} - \varvec{X} |} \end{aligned}$$

(16)

where the plastic multiplier $\dot{\lambda }$ is determined from the loading/unloading conditions:

$$\begin{aligned} \varPhi \le 0, \quad \dot{\lambda } \ge 0, \quad \dot{\lambda }\, \varPhi =0. \end{aligned}$$

(17)

The kinematic hardening $\varvec{X}$ is assumed to be additively decomposed into several hardening stresses $\varvec{X}_i$, i.e., $\varvec{X}=\sum _{i=1}^{n_X} \varvec{X}_i$. In the current paper, $n_X=2$ was found to be sufficient to capture the experimental data for rail grade R350HT. Following the proposal in [17], the evolution equation for the kinematic hardening is given by

$$\dot{\varvec{{X}}}_{i} = \frac{2}{3}C_{i} \dot{\epsilon }^{p} - \frac{{\gamma _{i}^{{m_{i} + 1}} }}{{C_{i}^{{m_{i} }} }}\left( {\sqrt {\frac{3}{2}} |{\varvec{X}}_{i} |} \right)^{{m_{i} }} \left\langle {\dot{\epsilon}^{p} :{\varvec{X}}_{i} } \right\rangle \sqrt {\frac{2}{3}} {\mkern 1mu} \frac{{{\varvec{X}}_{i} }}{{|{\varvec{X}}_{i} |}}$$

(18)

To increase the identifiability, the values of the parameters: $C_i$, $m_i$ and $p_i=\gamma _i^{m_i+1}/C_i^{m_i}$ are determined in the calibration against the experimental data.

Metamodel coefficients

Values of the metamodel coefficients for the elasto-plastic material model, calibrated using rail cross sections A and C, are listed in Table 7.

Table 7 Model coefficients $\varvec{\beta }$

Effect of elastic wheel assumption

The effect of assuming the wheel to be elastic, contrary to elasto-plastic behaviour, is illustrated in Fig. 21. The difference in the maximum pressure is less than 4.5%. The observed absence of difference in contact patch area is explained by a too small change in contact patch size in each direction, less than one element, to be resolved by the used mesh.

Fig. 21

Influence of elastic wheel on the distribution of normal contact pressure $p$ for $\varDelta y = 0$, $F_{n} = 111.5$ kN. Along the $x$-axis, only half of the plot is shown due to symmetry