65

Contact mechanics analysis of measured wheel–rail

profiles using the finite element method

T Telliskivi¤ and U Olofsson

Department of Machine Design, Royal Institute of Technology, Stockholm, Sweden

Abstract: A tool has been developed for contact mechanics analysis of the wheel–rail contact. Using
measurements of wheel and rail profiles as input, the tool is based on the finite element (FE) code ANSYS.
Traditionally, two methods have been used to investigate the rail–wheel contact, namely Hertz’s analytical
method and Kalker’s software program Contact. Both are based on the half-space assumption as well as on
a linear–elastic material model. The half-space assumption puts geometrical limitations on the contact.
This means that the significant dimensions of the contact area must be small compared with the relative
radii of the curvature of each body. Especially in the gauge corner of the rail profile, the half-space
assumption is questionable since the contact radius here can be as small as 10 mm. By using the FE method
(FEM) the user is not limited by these two assumptions. The profile measurement system Miniprof was
used to measure the wheel and rail profiles that were used as input when generating the FE mesh.
As a test case, a sharp curve (303 m radius) in a unidirectional commuter train track used by X1 and
X10 trains was chosen. The results of two contact cases were compared with the results of the Hertz
analytical method and the program Contact. In the first contact case the wheel was in contact with the rail
gauge corner. In the second case the wheel was in contact with the rail head. In both contact cases Hertz
and Contact presented very similar results for the maximum contact pressure. For the first contact case, a
significant difference was found between the FE method and the Hertz method and the program Contact in
all of output data. The Hertz and Contact methods both presented a maximum contact pressure that was
three times larger (around 3 GPa) than the FE solution. Here, the difference was probably due to the
combination of both the half-space assumption and the elastic–plastic material model. For the second
contact case, there was no significant difference between the maximum contact pressure results of the three
different contact mechanics methods employed.

Keywords: wheel –rail, elastic–plastic contact, finite element method (FEM)

1 INTRODUCTION model is linear–elastic and there is no friction between the

contacting surfaces. Furthermore, the contact is assumed to
Damage mechanisms such as surface cracks, plastic be a half-space. The half-space assumption puts geome-
deformation and wear can reduce the service life of a trical limitations on the contact inasmuch as the significant
railway track and the rolling stock (see Kalousek et al. [1] dimensions of the contact area must be small compared
for a survey of these damage mechanisms). In order to with the relative radii of curvature of each body. Especially
understand these damage mechanisms, the tiny contact area in the gauge corner of the rail profile, the half-space
where steel wheel meets steel rail needs to be investigated. assumption is questionable since the contact radius here
Traditionally, two methods have been used to investigate can be of the same size as the contact zone (10 mm). With
the rail–wheel contact, namely Hertz’s analytical method the Contact software [3] the contact zone is divided into
and Kalker’s software program Contact. The Hertz method cells and a boundary element analysis with a half-space
[2] assumes that the contact surfaces are smooth and can be assumption is used to solve the contact problem. Further-
described by a second-degree polynomial. The material more, a linear material model is used in the analysis.
Results previously presented by Cassidy (Hertz) [4] and
Knothe et al. (Contact) [5] reveal contact stresses up to and
above 3 GPa where the wheel is in contact with the gauge
corner section of the rail. Other examples of tools used in
The MS was received on 2 May 2000 and was accepted after revision for the contact mechanics analysis of the wheel– rail contact
publication on 24 August 2000.
¤ Corresponding author: Department of Machine Design, Royal Institute are Hertzian springs [6] and the Winkler mattress software
of Technology, S-100 44 Stockholm, Sweden. Fastsim by Kalker [3].
F01300 # IMechE 2001 Proc Instn Mech Engrs Vol 215 Part F
66 T TELLISKIVI AND U OLOFSSON

1.1 Rail analysis in curves 1.2 The conditions of finite element analysis
The form change of curves can be large over time [7, 8]. The finite element method (FEM) is not limited by the
Figure 1 shows the form change of a UIC 60 high rail over half-space assumption in contact analysis or by the linear
a period of 2 years in a narrow curve of a commuter train material model used in other structural mechanics methods.
track. As part of an investigation, a commuter train track Despite the fact that contact problems are computationally
was studied over a period of 2 years and the form and the difficult to solve in finite element (FE) codes, great
hardness of the track were characterized as two-dimen- developments in the field of numerical algorithms have
sional profile and surface hardness measurements. New been seen over the past few years [9, 10]. It is the objective
rails of 20 m apiece were inserted in two narrow curves. A of the present research to develop a simulation tool for the
length of the old rail was left in place as the test rail, wheel –rail contact based on the FEM. Since the original
enabling the study of new as well as 3-year-old rail. The wheel and rail profiles differ from those exploited on the
experimental results of the form measurements show that track, the tool should be able to use measured rail and
there was a significant change in the rail profile due to wheel profiles as input data when generating mesh.
wear as well as to plastic deformation and that both
processes influence the form of a rail that has been in use
for more than 5 years. The surface hardness measurement 2 SIMULATION OF THE WHEEL–RAIL
shows that the hardness of the new rail increases, but after CONTACT
2 years’ use the rail has not yet reached the hardness of the
old rail. These experimental results show that plastic
Figure 2 is a schematic illustration of how the wheel–rail
deformation is a necessary part in the wheel–rail contact
surface degradation problem may be solved. The data
analysis.
necessary for the FE analysis come from three sources. The
model is built entirely with the ANSYS solid modelling
tool. Starting from the right in Fig. 2, the material
properties must be decided on and the appropriate
constitutive laws supported by ANSYS must be selected.
This includes the plasticity models with their yielding
criteria and hardening rules. Here, kinematic hardening
criteria and the data on the stress–strain curves provided by
the wheel and rail manufacturers were used (see Table 1
and Fig. 7). The next important input was the prediction of
global forces from the train curving attached to the wheel
centre, usually predicted by multibody system (MBS)
programs such as GENSYS and MEDYNA. All the force
components are dependent on the situation in contact, the
solution algorithms and influence from outside the contact.
The MBS program and FE method results of the so-called
‘creep’ and other force components do not correlate; the
MBS data should therefore be used with care.
Data on the friction coefficient are produced with field
instruments such as the Salient system tribometer. In test
trials, the wheel and rail shape were measured with a
Miniprof instrument. The Miniprof data measured in
Cartesian coordinates are converted to the ANSYS pro-
gramming language, with keypoints created at the wheel–
rail contact surface. Additional keypoints and their connec-
tions by lines (splines, in the case of a curved region) form
the two-dimensional areas (Fig. 3). The three-dimensional
model is extruded rotationally for the wheel and linearly
for the rail. Separate volumes are created for the contact
localities. This helps to define the volumes by different
Fig. 1 Form change of a UIC 60 high rail in a 303 m curve over
meshes and material properties. Appropriate areas on the
a period of 2 years. Top: solid line, 3-year-old rail at test contact locality volumes are defined as ‘superelement
start; dotted line, after 1 year of use; dashed line, after 2 connection areas’ or ‘contact areas’. For the first, the
years of use. Bottom: solid line, new rail at test start; superelement volume shares the area with the contact
dotted line, after 1 year of use; dashed line, after 2 years locality volume. Special computing in ANSYS permits the
of use creation of the superelement of the volume far from the
Proc Instn Mech Engrs Vol 215 Part F F01300 # IMechE 2001
 CONTACT MECHANICS ANALYSIS OF MEASURED WHEEL–RAIL PROFILES USING THE FEM 67

Fig. 2 The wheel–rail analysis simulation tool

Table 1 Forces applied to the centre point of the wheel, and contact closure and meshed by the tetrahedral elements
material data for the wheel and rail providing the linear–elastic material properties with three
translational degrees of freedom. The contact locality is
Force (N) or
Load moment (Nm) meshed by the quadrilateral shape of elements using non-
linear material properties. The rail is completely fixed at
Fx ¡14 962
Fy ¡33 200 the bottom, while the wheel has a so-called ‘steering centre
Fz ¡73 200 node’, which is rigidly fixed to the closest surface nodes at
My 6900 the centre hole of wheel. The length of the contact area in
Mz 600 the rolling direction is 50 mm, permitting rolling if
Yielding point (MPa) Rail: 469 Wheel: 670 necessary. In general, the FE model consists of around
Ultimate stress limit (MPa) Rail: 606 Wheel: 780 40 000 nodes. The contact zone is defined by 5000
elements both for the rail and for the wheel, with 3000
contact elements between them. To reduce the computation
time, a superelement technique is utilized to minimize the
number of degrees of freedom of the model. The computa-
tion time is 48–60 h and the results file takes up consider-
able hard disk space (more than 400 MB). A more detailed
description of the contact model and its build in the
commercial program ANSYS is presented in Telliskivi et
al. [11]. As a test case, a sharp curve (303 m radius) in a
unidirectional commuter train track used by X1 and X10
trains was chosen. The results of two contact cases were
compared with the results of the Hertz analytical method
and the program Contact. In the first contact case the wheel
was in contact with the gauge corner. In the second case the
wheel was in contact with the rail head.

2.1 Contact kinematics

In general, the two steps that have to be followed to set up
the contact geometry are defining the contact and develop-
ing the local kinematic relations including overall solution
strategies. With both new and older rail/wheel shapes, the
location(s) of contact must be found as well as the contact
Fig. 3 Different regions in the FE model locality with the six-degree-of-freedom discretization solid
F01300 # IMechE 2001 Proc Instn Mech Engrs Vol 215 Part F
68 T TELLISKIVI AND U OLOFSSON

elements. Thereafter, the contact elements on the surface Fig. 4). The negative magnitude of the gap, g, indicates the
have to be build. contact and the forces are developed in a direction normal
For the flexible–flexible contact analysis, the three- (n direction) to the target, which will tend to reduce the
dimensional contact element CONTAC49 was used. This penetration to an acceptable numerical level. In addition,
allows large deformations and sliding and enables repre- frictional forces are developed in directions that are tangent
sentation of the general contact of models that are gener- to the target plane. The normal and tangential forces are
ated with arbitrary meshes. Its use is not limited to known referred to by an x– y– z system. The penalty method,
contact or node-to-node configurations. The contact situa- whereby
tion is delineated by various algorithms. Two potential »
contact surfaces are referred to as either the ‘target KN g if g < 0
Fn ˆ (1)
surface’, with nodes I, J, K and L, or the ‘contact surface’, 0 if g . 0
with node M (see Fig. 4). Usually there are several contact
nodes associated with one target surface. Contact occurs enforces compatibility by means of a contact stiffness for
whenever one or more contact nodes penetrate the target the normal forces. F n is the normal load, g is the
surface. The contact location is computed by an iterative penetration and K is the penalty factor. Friction causes the
Newton method based on a normal projection of the contact tangential forces and requires the representation of the
node onto the target plane. At the projected contact point coefficient of friction. Constitutive equations in the contact
the value of the gap, g, is determined by the contact node’s area lead to elastic Coulomb friction.
location with respect to the target plane. Contact penetra- It is necessary to calculate the tangential deformations of
tion is assumed to occur if the value of g is found to be the contact node relative to the target. The deformation is
negative, and the contact locations are found to be in the first separated into x and y components. Next, the
natural space bounds of the target. The wheel is copied to deformation is decomposed into elastic (or ‘sticking’) and
the closure of the rail (approximately 0.1 mm). It turns out sliding (‘inelastic’) components. The actual computation
that knowledge about the exact contact location is not that is performed uses a technique similar to that of the
important if the angular attitude in earth-fixed coordinate non-associative theory of plasticity. For each substep in
system of the wheel is known. This is worked out by means which sliding friction occurs, an elastic predictor is
of the GENSYS multibody dynamic simulation program. computed in contact traction space. A special situation
The linear shift initiates the contact and thereafter stepwise arises when a contact node moves from one target to
loading is performed until the maximum load is reached. another. When this occurs, the contact history is passed
Since contact conditions are represented as inequality from the target that was in contact to the target that is
constraints, the kinematic relations of the motion in the currently subjected to contact.
contact area lead to variational inequalities. Different In the basic Coulomb friction model, two contacting
possibilities exist for the numerical solution of these surfaces can carry shear stresses up to a certain magnitude
problems. Here, the penalty method is used to calculate the across their interface before they start sliding relative to
constraint equations in respect of a normal direction (n each other. This state is known as ‘sticking’. The Coulomb
direction; see Section 2.2). For calculation of the tangential friction defines an equivalent shear stress at which sliding
equations, the relevant constitutive relations are needed, as on the surface begins as a fraction of the contact pressure
discussed in the following section. [12]. Once the shear stress is exceeded, the two surfaces
will slide relative to each other. This state is known as
2.2 Contact forces ‘sliding’. The goal is to determine at what point there will
be a transition from sticking to sliding, and vice versa. The
Contact is indicated when the contact node M penetrates upper estimate for shear stress before sliding will occur is
the target surface defined by target nodes I, J, K and L (see the von Mises yield stress of the material adjacent to the

Fig. 4 Contact identification in FE analysis

Proc Instn Mech Engrs Vol 215 Part F F01300 # IMechE 2001
 CONTACT MECHANICS ANALYSIS OF MEASURED WHEEL–RAIL PROFILES USING THE FEM 69

surface. Spin angular torque redistributes the tangential 3 RESULTS

forces in relation to the pure, normal loading. Every time-
step in rolling determines the proportion of elements that In this section for FE analysis tool is described and the
undergo sliding and those that will stick. In a steady state preliminary results, along with a comparison between the
condition, that is, where no time parameters are included, main methods currently available, are outlined. First, the
the relations between the stick and slip regions do not difference between the various methods is presented (see
change. Table 2 and Fig. 6). The main aim of this work was to
enhance knowledge of the contact pressure and maximum
stresses in bulk material. This should provide an appro-
priate basis for the study of the degradation mechanism and
2.3 Test cases wear simulation.
The results for case 1 could not be compared with the
As a simulation case a sharp curve in a commuter train Hertzian results since the Hertzian method assumes one-
track serving the Stockholm area was chosen. The track point contact. The Hertzian solution showed approximately
carries almost exclusively unidirectional commuter trains twice the contact length in the y direction calculated by the
with an average speed of 75 km=h. Two types of vehicles other methods. The results remained the same with
are used, X1 and X10 trains, both operating in pairs, with variation of the radius of curvature as an important
one powered unit and one trailing unit. This track and the parameter in the range of 0.01–0.015 m for the rail and
rolling stock have been studied previously as part of a 0.02–0.1 m for the wheel. Here, as with Contact, the
national Swedish transport programme in which both rail normal force was split into three parts and the results were
and wheel profiles were measured over a period of 2.5 very similar for the two classic methods developed by
years [13]. Furthermore, the trains used in this study have Hertz and Kalker (see Fig. 6). The difference for the FE
been previously modelled with the train dynamic simula- method was approximately 300 per cent for the contact
tion software GENSYS [14], thus supplying the necessary area and .200 per cent for the maximum pressure.
input data on forces and wheel attitude at optional contact For case 2 the calculated difference in maximum contact
locations. pressure and contact area was small where the minimum
Two load cases were used to study the model. Both were contact radius was large compared with the significant
originated through simulation of the X1 powered unit and dimensions of the contact area. Since the half-space
represent the first and second wheelsets in the leading assumption is valid in this case, the difference of 25–30
bogie. The coordinate system chosen was similar to that of per cent was probably due to the plastic material model
the Deutsche Industrie Normen (DIN), with positive used in the FE analysis.
vertical (z) coordinate upwards, y to the left and x positive Figure 8 illustrates the stress distribution of a fairly
to the train motion direction. Figure 5 presents the contact worn, 2-year-old rail profile. The contact is separated into
point locations on the wheel and rail. Table 1 shows the the two different types of contact. The flange contact is
forces applied to the centre point of the wheel. From a
geometrical perspective, the two load cases represent the
contact points with a large difference in the curvature of
contacting bodies. In load case 1 the minimum contact
Table 2 Comparison between FE and Contact results on maxi-
radius was about 20 mm; in load case 2 it was about mum von Mises equivalent yield stresses
300 mm. New rail profiles and a wheel profile from an X1
train that had been in use for 2 years were used in the two Method Load case 1 (MPa) Load case 2 (MPa)
load cases. Finally, two calculations with more degraded FE with plasticity 606 577
rail profiles were made to demonstrate the tendency of Contact (half-space 3057 715
assumption)
stress distribution changes over time (see Figs 8 and 9).

Fig. 5 Contact point location for the two test cases

F01300 # IMechE 2001 Proc Instn Mech Engrs Vol 215 Part F
70 T TELLISKIVI AND U OLOFSSON

Fig. 6 Comparison, with respect to maximum contact pressure and the contact area, between three different
contact mechanics analysis methods

tion shows that even with the three contact locations, the
contact area is too small and the train forces are too large
to keep the bulk material in the linear–elastic stress region
(see Fig. 9). At the very beginning (case 1), the plastic
work is extremely high. The maximum equivalent von
Mises stress exceeds even the ultimate stress limit, which
was 606 MPa for the actual plastic model (see Table 2).
Here, the additional increase in strain does not result in an
increase in stress and the material behaves perfectly
plastically. The plastic flow hardens the material and makes
the contact more conformal. In the continuing process,
other wear mechanisms will thus become more significant.

Fig. 7 Non-linear material models for the rail (curve A) and

wheel (curve B) portions of the contact region; stress
(MPa) versus strain per cent 4 COMMENTS AND FUTURE WORK

The FE analysis calculates the discretized pressure and

clearly separated from the contact with the rail head. In displacement distribution on the surfaces of contacting
curves, the force increases because of the centrifugal bodies. Simulating the steady curve, the loads and the
acceleration. Because here the rail set is skewed, the locally fields of strain, stress and deformations in the contact
vertical load to the rail also increases. The stress distribu- region do not change with time. The contact region is
Proc Instn Mech Engrs Vol 215 Part F F01300 # IMechE 2001
 CONTACT MECHANICS ANALYSIS OF MEASURED WHEEL–RAIL PROFILES USING THE FEM 71

Fig. 8 Stress distribution (von Mises equivalent stresses) and reaction forces in a wheel run on a rail that has been
in use for 2 years

Fig. 9 Stress distribution (von Mises equivalent stresses) in 3-year-old rail and a 3-year-old wheel

continuously being renewed by the rolling surface because contact passes through the contact field in a manner
the rail and wheels are of the same form in the rolling determined by the relationship between the creep and ‘spin’
direction [15]. After the overrun, the profile is changed as a motions. The separation of normal and tangential forces,
result of plastic flow and wear. This assumes, however, that and the law of friction and its magnitude, affect the contact
during the contact every material point or node that is in balance, which can be conformal, skewed in space and/or
F01300 # IMechE 2001 Proc Instn Mech Engrs Vol 215 Part F
72 T TELLISKIVI AND U OLOFSSON

occur simultaneously at several locations. This causes a ACKNOWLEDGEMENTS

very complex relationship between the motions and all the
forces of the wheel, which are, however, uniquely defined This work was performed as part of the Swedish research
[3]. The results are of value for current and possibly future programme SAMBA. The authors gratefully acknowledge
train dynamic simulations using other contact methods the financial support of the Swedish National Board for
such as Contact and Hertz. For curving analysis and other Industrial and Technical Development (NUTEK), Adtranz
train dynamic studies these contact methods are probably Sweden AB, Stockholm Local Traffic, the Swedish Na-
sufficient, but, if the results from the contact mechanics tional Rail Administration and the Swedish State Railways.
models are used in fracture mechanics, fatigue or wear
studies, more care is necessary in the selection of contact
mechanics methods. This study shows that there can be a REFERENCES
large difference in results between half-space-based meth-
ods (Contact and Hertz) and non-half-space-based methods 1 Kalousek, J., Masel, E. and Grassie S. Perspective on
(FEM). A further study is therefore necessary to investigate metallurgy and contact mechanics. In International Heavy
the limit of the half-space assumption, in particular with Haul Association STS Conference on The Wheel/Rail Inter-
reference to rail wheel applications. face, Moscow, June 1999.
In many situations wear is a significant damage mechan- 2 Johnson, K. L. Contact Mechanics, 1985 (Cambridge
ism. Nilsson [13] showed in an experimental full-scale test University Press, Cambridge).
3 Kalker, J. J. Three-dimensional Elastic Bodies in Rolling
that the losses in rail cross-sectional area remain approxi-
Contact, 1990 (Kluwer, Dordrecht).
mately constant over exploitation time. The wear rate can
4 Cassidy, P. D. Variation of normal contact stresses for
be simulated if the contact state history is known. A wear different wheel/rail profile combinations. European Rail
model includes a wear coefficient, which has a scientific Research Institute, ERRI D 173/DT 336, 1996.
basis; see, for example, Archard [16] and Lim and Ashby 5 Knothe, K., Theiler, A. and Güney, S. Investigation of
[17]. Other important variables in the wear models are the contact stresses on the wheel/rail-system at steady state
contact pressure distribution and the accumulated sliding curving. In 16th IAVSD Conference, Pretoria, South Africa,
distances. Proper estimation of these variables also en- 30 August–3 September 1999.
chances the selection of the right wear regime for each 6 Pascal, J. P. and Sauvage, G. The available methods to
point that has passed the contact zone. Further improve- calculate the wheel/rail forces in non Hertzian contact patches
ment of the calculation tool presented here is therefore and rail damaging. Veh. Syst. Dynamics, 1993, 22, 263–275.
7 Olofsson, U. and Nilsson R. Initial wear of a commuter train
directed to the calculation of the accumulated sliding
track. In Nordtrib’98, Ebeltoft, Denmark, 8–10 June 1998.
distances.
8 Olofsson, U. On the form change due to wear and plastic
deformation in narrow curves on a commuter train track
(submitted).
9 Wriggers, P. Finite element algorithms for contact problems.
5 CONCLUSIONS Arch. Comput. Meth. Engng, 1995, 2, 1–49.
10 Aliabadi, M. H. and Brebbia C. A. Computational Methods
A tool for contact mechanics analysis of the wheel–rail in Contact Mechanics, 1993 (Elsevier, Amsterdam).
contact has been developed that allows easy changing of 11 Telliskivi, T., Olofsson, U., Sellgren, U. and Kruse, P. A tool
and a method for FE analysis of wheel and rail interaction. In
the geometry of the contact. Measured wheel and rail
ANSYS Conference, Pittsburgh, Pennsylvania, 2000.
profiles were used in the model generation. Neither the
12 Refaat, M. H. and Meguid, S. A. On the modelling of
half-space assumption nor a linear material limits the FE frictional contact problem using variational inequalities. Finite
model, while Hertz’s analytical method and the Contact Elements in Analysis and Des., 1995, 19, 89–101.
program, which uses the well-known boundary element 13 Nilsson, R. Wheel and rail wear—measured profiles and
method, are both limited by these assumptions. The results hardness changes during 2.5 years for Stockholm commuter
of two test cases presented show that the difference in traffic. In Hannover EUROMECH 409 Colloquium, Hannover,
maximum contact pressure between the Contact and Hertz Germany, 6–9 March 2000.
methods and the FEM was negligible where the radii of 14 Jendel, T. Prediction of wheel profile wear—comparison with
curvature of two contact bodies at the contact point were field measurements. Presented at the Contact Mechanics and
large compared with the significant dimensions of the Wear of Rail/Wheel Systems Conference, Tokyo, Japan, 25–
27 July 2000.
contact area (in other words, the half-space assumption
15 Rathore, S. K. and Kishore, N. N. Finite element analysis of
was valid). However, in test case 1, where the radius of rolling contact problems using minimum dissipation of energy
curvature of the rail edge was small compared with the principle. J. Appl. Mechanics, December 1997, 64, 271.
significant dimensions of the contact area, the difference 16 Archard, J. F. Contact and rubbing of flat surfaces. J. Appl.
between the model used here and the Contact/Hertz meth- Physics, 1953, 24, 981–988.
ods was as large as 3 GPa, and was probably due to both 17 Lim, S. C. and Ashby, M. F. Wear mechanism maps. Acta
the half-space assumption and the material model. Metall., 1987, 35(1), 1–24.

Proc Instn Mech Engrs Vol 215 Part F F01300 # IMechE 2001