Multibody simulation benchmark for dynamic vehicle–track interaction in switches and crossings: results and method statements

A Benchmark of railway multibody dynamics software application to switches and crossings (S&C) is presented, comparing all major commercially available software and a few independent codes. Two different representative S&C have been implemented, using the Manchester Benchmark passenger vehicle. The final results show that all software offer a reliable and efficient way to understand the kinematic and dynamics forces between the wheels and the track elements. The highest challenges are found when modelling a combination of multiple rails in simultaneous contact with a wheel (check-stock or switch-stock), large longitudinal variations in rail shape (crossings) and high lateral steering forces (diverging cases in tight radius). In those cases, the codes able to account for the exact relative motion of each wheels with respect to each rails independently are the most apt. The most significant variations between software are found in the contact prediction with an influence on the detailed contact tangential and normal forces. The user variability is found to be very small, with the most time-consuming and error prone being the task of handling the input data for the variable rails definition. All software could benefit from improvements to assist the user and ensure higher reliability and efficiency generally.

of this paper is to present and compare the results submitted by participants in response to the Benchmark task together with an overview of the corresponding modelling and simulation techniques they used to complete it. Full method statements from participants are available online [1]. For the reader interested to perform the exercise after this publication, the task description and data sources required to do so are also available [2,3].

Background and motivation
There have been a number of simulation Benchmarks performed in railway mechanics since the late 1990s on, for example, vehicle dynamics, vehicle-track interaction at high frequencies and more recently on longitudinal train dynamics, contact and vehicle pantograph-catenary interaction [4][5][6][7][8]. A Benchmark on S&C has, however, been missing which led to the present initiative. S&C merit the attention of a tailored Benchmark as they constitute some specific challenges in terms of modelling and simulation of dynamic vehicle-track interaction: 1) In S&C, there are large and sudden changes in rail profile geometry along the track.
This constitutes a challenge in terms of rail surface geometry modelling. Given the stiffness of the wheel/rail (w/r) contact, the slightest distortion in the rail surface description can induce significant shifts in contact conditions and result in large dynamic contributions to the w/r contact forces. 2) Wheels passing through a switch or a crossing panel can make simultaneous contact with multiple rail bodies that can deform and displace relative to one another, i.e. the stock rail and switch rail in the switch panel and the check rail and stock rail in the crossing panel. This calls for more elaborate track and w/r contact modelling compared to plain line. 3) Due to the varying rail and track sections throughout S&C, track properties will vary along the track by design.
For this Benchmark, participants have been given the task to model the switch panel and crossing panel for two different S&C designs and to simulate dynamic vehicle-track interaction in those panels. Rail geometry data have been provided in the form of discrete cross-sections and the track properties are represented using co-running track models with specified properties. Traffic is represented by the passenger vehicle from the Manchester Benchmarks [4].
In doing so, this Benchmark is foremost addressing point (1). This is because the greatest source of results variation between participants and modelling approaches is expected to stem from the rail geometry and how it is represented between the given cross-sections. Point (2) is accounted for in the Benchmark as the track model features individual bodies for each independent rail. The challenge here is to demonstrate the modelling and simulation capability for a track model with this topology, allowing for simultaneous multiple points of contact of the wheel onto the track components. Point (3) is accounted for via separate track properties for the switch and crossing panel, but the continuous variation during simulation is not addressed in this Benchmark.
In addition to addressing these particular S&C features in simulation and comparing the obtained results between different software and modelling approaches, the S&C Benchmark has also (a) contributed to the creation of a set of reference simulation cases for dynamic vehicle-track interaction in S&C using multi body simulation tools and (b) knowledge sharing and advance of the state-of-the-art in S&C simulations within the railway dynamics community which should ultimately help improve the railway sector as a whole. It has to be noted, however, that this benchmark does not constitute an absolute reference for the validation of vehicle dynamic interaction in S&C as the participants results are compared amongst themselves for a set of nominal simulation cases and not against actual physical measurement of the same situation.

Overview of the Benchmark simulation cases
The simulation cases are listed in Table 1, where cases based on the 56E1 rail section comes from a UK design and the one based on the 60E1 from a Swedish design. The switch and crossing panels have been evaluated separately to prevent differences in modelling and simulation in the switch panel to propagate and give different initial conditions once the vehicle enters the crossing panel. To simplify modelling, the varying S&C rail profiles are always located on the right-hand side of the track, as shown in Figure 1. In this way, the only change that has been needed to change simulation set-up from the through to the diverging route has been a change in track layout and vehicle speed. Run 9 is introduced to allow for a baseline comparison of each participant's simulation set-up and vehicle model.
As the main purpose of the Benchmark has been to evaluate modelling and simulation of S&C rail profile geometry in different software, the requested results for the Benchmark submission are focused around the w/r interaction. A full description of the Benchmark simulation cases is available in the Benchmark statement [2]. Table 2 lists the participating institutions and software used. There is a total of 19 sets of results submitted by 18 independent participants (TTCI provided two sets of results using a different w/r contact coupling approach), using nine independent software, two of which being research programmes (SDITT and MUBODyn), while the others are commercially available and all but VOCO participated in the original Manchester Benchmark (VI-Rail was then ADAMS/Rail). For two of the software, there are multiple results sets submitted by different participants, seven for Simpack and four for VI-Rail. Amongst the Identical to Run #2 but with a constant 56E1 rail profile replacing the stock rail geometry and the switch rail geometry being removed (retaining the same track formulation)  participants, there are 10 university/research institutes, while the rest are either software developers and/or consultancy companies. While this benchmark was initially carried out blind by all participants, there was then a period of consultation during which interim results were made available, leading to revision and improvements, mainly because of errors of interpretation or implementation. The greatest differences in method and results between participants in the Benchmark are found between software. This is natural as some aspects of the method used to simulate the Benchmark cases are software rather than modeller specific. For the sake of presentation, two groups of participants have therefore been defined. The first one is software lead participants consisting of software developers and VDG as the sole user of Vampire. The second group consists of the users of Simpack and VI-Rail.

Compilation of method statements
This chapter is a compiled summary and comparison of the modelling and simulation methods employed by the Benchmark participants to complete the exercise. The information comes from participants' method statements that are formulated around the questions asked in the Benchmark statement [2]. The compilation is focused on the principal ideas used to model the Benchmark cases and is therefore by necessity succinct. For full details on each participants' contribution, see their full method statement [1].

Geometry implementation
The Benchmark rail geometry is provided as sets of transversal rail cross-sections with specified positions along the track. This section will give an overview of how participants built the 3D rail geometry from this input in their software. The discussion will be broken down into two parts: (a) the implementation of individual 2D cross-sections and (b) the construction of 3D rail sections and their representation in time domain simulations. Frequently used terms in this section are defined in Table 3.

Two-dimensional cross-sections
Among the Benchmark submissions, there are three methods for representing individual cross-sections in software: (a) spline fit, (b) discrete data points and (c) a fitted chain of circle segments. Regardless of the representation format in the software, modellers are faced with implementation choices. These can mainly be put into two categories: (a) smoothing or alteration of the profile shape or discretisation, particularly with a focus on the smoothing of the running surface for more regular contact conditions in the case of irregular curvature variations, and (b) cropping or splitting of the profiles in order to control the build-up and interpolation of 3D rail sections. The methods used by software to represent the crosssections are listed in Table 4. Comments in the implementation are found in Table 5 for software participants and in Table 6 for software users. The modifications include different levels of applied smoothing in the spline fit, different discretisation steps and one instance Table 3. Definitions of frequently used terms in the description of rail geometry implementation Term Definition (rail) cross-section or profile A 2D transversal cross-section of the rail head described by a curve defined by discrete data points. As supplied with the Benchmark input data set or interpolated from the provided data (rail cross-section) segment A partial length of a rail cross-section curve (rail) section A local length of rail with specified start and end points along the track (rail) section break An interruption (stop/start) in the otherwise continuous rail profile interpolation in the longitudinal direction (rail) running surface The rail surface area(s) where w/r contact is expected during vehicle running. Also applicable to the corresponding segment(s) on individual cross-sections The profiles are read as discrete data points.  of rounding the top of the 56E1 crossing profiles. The second category of modifications will be covered under 3D geometry representation. Minor comments are made with regard to the wheel profile implementations. These can be found in the method statements.

Three-dimensional geometry representation
Among the software in the Benchmark, there are four principal methods to represent rail geometry between given cross-sections in time domain simulations: (a) to perform geometrical interpolation between the provided cross-sections online during simulations to obtain the profile at a given longitudinal wheel position (the majority of software). (b) To  solve the w/r contact problem for cross-sections in advance for a range of contact conditions, tabulate the results, and then interpolate between the contact lookup tables online  during time domain simulations (NUCARS-WNT, Vampire). (c) Is analogous to the previous, but instead of contact tables the parameters needed for the contact calculations are tabulated instead and interpolated online to calculate the contact conditions (VOCO). (d) To implement the geometry as a series or rail sections where each section has a constant profile (NUCARS FIT).
For the software that interpolate on the geometry, the interpolation orders range from linear to cubic spline. In addition to the geometry interpolation method used by the software, the resulting geometry is highly dependent on how the modeller chooses to define the reference sections that specify the interpolation paths. The definition of the longitudinal geometry interpolation is most critical where there are large geometric changes from one cross-section to the next and where a direct interpolation would lead to a severely warped profile. In practice, these changes mainly concern the crossing geometries that feature several step changes in cross-section geometry along its length. Figure 2 shows the rail cross-sections corresponding to the 56E1 crossing. In this figure, the alternating colours for the rail sections visualise the two principal ways that modellers have defined distinctive rail sections for geometry interpolation. In the left plot, the crossing is divided into several longitudinal rail sections with interpolation breaks in between. In the right plot, the geometry is modified and split into two rails, i.e. the wing rail and the crossing nose, and then the geometry interpolation is performed separately on each continuous rail body which avoids the necessity of interpolation breaks at the beginning and end of the distance where the crossing nose and wing rail overlap in the longitudinal direction.
In order to create the rail sections in Figure 2(a) from the given rail geometry data, the first task is to identify step changes in cross-section geometry between adjacent crosssections. This can either be done via visual inspection or, as in MUBODyn, via a tolerance for when the profile arc length or another geometric measure of two adjacent cross-sections differ by a certain amount. To build the section break, the profile with the longer arc length is trimmed such that the remaining profile shape matches the profile with the shorter arc length. At the wing rail to crossing nose transition, for example, this means that the first cross-section including both the crossing nose and wing rail is trimmed such that only the wing rail remains and that its overall shape matches the preceding wing rail profile. Positioning the trimmed profile at the same longitudinal coordinate as the original profile, two cross-sections are defined at the same location. This allows for one geometry interpolation leading up to this location and another one leading from there while not introducing any irregularities in the running surface (i.e. continuous contact conditions) as the two cross-sections match one another exactly on the overlapping segment.
The software that interpolate between contact tables (Vampire, NUCARS-WNT) and contact parameters (VOCO) use linear interpolation between cross-section data. Examples of how the interpolation can be controlled with this approach are analogous to those found for the interpolation of geometry (a) via the selection of a subset of the provided rail cross-sections for interpolation and (b) the introduction of additional sections in certain locations to guide the interpolation and (c) the introduction of separate rail bodies for two adjacent rails (e.g. crossing nose and the wing rail) to perform the interpolation separately for these two bodies. This latter approach is analogous to the geometry implementation in Figure 2(b). For more information on the use of contact tables in MBS simulations, see e.g. [9,10]. For comments on the use of piecewise constant cross-sections, see the NUCARS FIT entry in Table 8. Table 7 presents how software account for the 3D rail geometry in time domain simulations. Comments on the implementation are given in Table 8 for software participants and in Table 9 for software users.

Track model
The track model topology specified in the Benchmark is presented in Figure 3. It is a co-running track model where the same system of masses and bushings is replicated independently under each wheelset. For simulations in the switch panel, masses 1, 3 and 4 are active to represent the opposite stock rail, switch rail and main stock rail, respectively. For simulations in the crossing panel, masses 1, 2 and 4 are active to represent the opposite stock rail, the check rail and the crossing rail, respectively. The ground reference is a track-following coordinate system running along with each wheelset and track model Table 7. Software methods to represent rail geometry between given cross-sections in time domain simulations Software Method GENSYS Linear longitudinal interpolation between adjacent cross-sections to obtain the profile at the contact location. The interpolation is performed between corresponding profile segments (not necessarily the full profile on both ends) to avoid profile warping at discontinuous changes in rail profile shape. Calculations are performed online. The method description concerns the w/r coupling creep_fasim_4 used in the Benchmark MEDYNA Longitudinal interpolation of cross-sections using Bezier splines to obtain the profile at the contact location. Calculations are performed online Simpack ditto VI-Rail ditto MUBODyn Longitudinal interpolation of the given cross-sections using cubic splines prior to analysis. These splines are then used to interpolate cross-sections (also described by cubic splines) online during simulations NUCARS Penetration model (FIT): Piecewise constant cross-section changing for each discretely given profile for online contact calculations Rigid contact model (WNT): Discrete cross-sections used to generate w/r contact tables. Linear interpolation between contact tables in time domain simulations SDITT Linear longitudinal interpolation between rail cross-sections to obtain the profile at the contact location. The interpolation is performed between corresponding profile spline segments that each have the same number of points in their discretisation. Calculations are performed online Vampire Discrete cross-sections are used to generate w/r contact tables.  Only a subset of the provided cross-sections that are deemed relevant for the rail geometry description are used to avoid interpolation and discretisation issues. Existing sections are duplicated and shifted laterally and longitudinally when necessary to create intermediate contours for contact table interpolations  VOCO The wing rail and crossing nose are discretised individually and give their own interpolation of contact parameters even though they are rigidly connected system. Most software uses a joint or constraint definition for each track mass and rail body that allow for displacement and rotation in a plane orthogonal to the track tangent direction while the longitudinal position along the track is prescribed by the vehicle speed (see details in Table 10).
For the participant entries where the software support the implementation of the full Benchmark track model (the majority of cases), the main difference in the implementation stem from whether the rail's rotational degrees of freedom (DOFs) are constrained to the track mass via constraints or stiff rotational springs, or whether they are constrained with Table 9. Comments on the implementation of 3D rail geometry for software users Simpack user Implementation details Chalmers Section breaks are introduced at rail profile discontinuities Polito ditto ViF ditto PROSE Section breaks are introduced at rail profile discontinuities Certain profiles are divided into multiple lateral sections where each section is associated with its own longitudinal interpolation. Intermediate profiles are created at some locations that possess characteristics of both adjacent profiles TUB Section breaks are introduced at rail profile discontinuities All stock rails are truncated below the gauge corner line to ensure a consistent interpolation UoB Section breaks are introduced at rail profile discontinuities 60E1 crossing profiles are trimmed at the ends for improved interpolation.
VI-Rail user

D2S
The provided rail cross-sections are imported in the software SpaceClaim to create rail surfaces that account for rail discontinuities. These surfaces are then discretised to generate input for VI-Rail GMUNIFI Section breaks are introduced at rail profile discontinuities UoHvi ditto respect to the track-following coordinate system. For entries where alterations are required due to software limitations, simplifications include, for example, the use of massless rails and one rail body per track side (Vampire (VDG), NUCARS-WNT) or the coupling of the check rail contact to ground (Vampire (VDG), GMUNIFI, NUCARS FIT). The track model implementation for software participants is presented in Table 8. The software user variability is presented in Table 11.

Comments on the calculation of track preloads
During the comparisons of initial results for this Benchmark exercise, a significant variability in results was observed for the vertical wheel position and relative w/r movement in the switch panel simulation cases due to the wheel transition from the stock rail to the switch rail. The origin of these discrepancies could be traced to different ways of initialising simulations and the calculation of track and w/r interface preloads.
In the set-up of MBS simulations, it is common that preloads are introduced in the model's force elements. In contrast to a static equilibrium simulation where the equilibrium of the system is sought in terms of system displacements, the preloads are calculated to achieve full or partial equilibrium of the system in a given (typically nominal) configuration. ditto NUCARS-WNT and FIT The rails and track structure are modelled as massless, with the rails coupled directly to the ground in the vertical direction. For this Benchmark study, an optional feature was used that inserts a massless track body between the rails and ground in the lateral direction. One massless body is used to represent the switch and stock rail, and one massless body is used for check rail and stock rail, respectively. The track model uses equivalent stiffness and damping between each rail and the ground and/or tie, and between tie and ground in the vertical and lateral direction, with the lateral stiffnesses corresponding to the Benchmark track model Vampire (VDG) The track model implementation includes the following simplifications: (a) there is only one rail body on each side of the track and these bodies are massless and (b) the check rail is mounted directly to ground and is modelled as a gap function between the wheelset lateral displacement and the check rail contour VI-Rail The  Due to the pairs of rail bodies on each side of the track model, care must be taken in the preload definitions to avoid the introduction of any unintended offset in the relative positioning of the rails. If a preload is introduced to the rail bushing that is loaded at the start of simulations (i.e. a stock rail), and none to the adjacent unloaded bushing (i.e. switch rail or check rail), an offset would have been introduced in the positioning of the two rail bodies in the nominal unloaded configuration as the stock rail would move upwards due to the preload force. Another possible reason that could introduce an offset is the functionality available in some software that shifts the rail profile with respect to the rail body to engage the w/r contact and reach a specified preload. Both effects can be cancelled by either duplicating the computed preload settings to the adjacent rail or to avoid preloads in the track and w/r contact altogether and instead solve for the static equilibrium before time domain simulations. This discussion mainly concerns vertical preloads. Preloads in the lateral direction would be much smaller than the vertical and practically irrelevant from the perspective of simulation initialisation. They are also not recommended as they would effectively correspond to a track gauge change.

Wheel/rail coupling
The wheel/rail (w/r) coupling is here defined as the methods for determining the w/r contact point quantities (such as forces, contact patch dimensions, etc.) and their location on wheel and rail. This compilation focuses on two main aspects of the w/r coupling: the modelling of individual contact patches and what aspects of relative wheel and rail movement are considered in determining the contact point locations.

Wheel/rail contact modelling
The theories used to model normal and tangential contact and the resulting forces for each contact point are listed in Table 12 for software participants and in Table 13 for software users. GENSYS, Vampire (VDG) and VOCO opted for a simplified modelling of the flangeback to check rail contact due to the flat and almost vertical nature of the check rail contact surface. They used an equivalent spring-damper for normal contact and a simple creepforce model for tangential force calculation with assumed full slip. Contact modelling parameters for all software are listed in Table 14. Further elaboration on wheel-rail contact modelling in S&C is presented in [11].Contact modellingparameters for all software are listed in Table 14 and for all participants in Table 15.

Wheel/rail contact detection
The purpose of this section is to give an overview of selected functionality in each software's w/r coupling algorithm that can have an impact on S&C simulation cases. As variability in the w/r coupling functionality is only found between software, there is no user variability to report. It is investigated whether the following aspects are accounted for in the w/r contact search: 1. Wheel(set) roll 2. Rail roll 3. Longitudinal contact point shifts with respect to the wheel centre due to wheelset yaw 4. Effective change in wheel profile contact geometry due to wheelset yaw A fifth aspect would be to account for a longitudinal shift in contact point location due to changing rail profile geometry. It is not listed as a dimension of comparison; however, as all software use a single rail cross-section for all contact points at a given longitudinal wheel position and thus do not account for this effect.
Wheel and rail roll, i.e. rotational movement in a plane perpendicular to the track direction, are illustrated in Figure 4(a) where α w is the roll rotation for the wheel and α r the roll rotation for the rail, respectively. While the wheel roll kinematics might be important to Table 12. Modelling of normal and tangential w/r contact for individual contact points among software participants Software participant Normal contact Tangential contact GENSYS One linear spring and damper per contact patch. Elliptical contact patch computed from normal force and wheel and rail geometry. The method description concerns the w/r coupling creep_fasim_4 used in the Benchmark [12] FASTSIM [13] MEDYNA Kik-Piotrowski [14] F A S T S I M [ 15,16]

MUBODyn
Hertzian approach that accounts for the restitution effect [17]. Equivalent elliptic contact patch computed from strip discretisation of the contact patch via an extension of the approaches in [9,18] Polach [19] NUCARS FIT Hertzian [20,21] Lookup creep-force table generated using CONTACT [22]. FASTSIM used for contact conditions outside of  Simpack's equivalent elastic contact * was used as the base choice for normal contact modelling. Chalmers and PROSE also used the discrete elastic contact # in some instances. Chalmers used it to model the wheel to checkrail contact and PROSE used it for all contacts in Runs 5,6 and 8. The discrete elastic contact was used in these instances for more stable contact modelling against the check rail (the near-vertical check rail face can make contact search in the z-direction difficult).
FASTSIM [13,27] * The equivalent elastic contact computes an equivalent interpenetration of the bodies in contact corresponding to an equivalent contact ellipse [28]. The normal contact force is then computed according to Hertzian theory [26]. #The discrete elastic contact uses a STRIPES-based approach [24,25] Chalmers consider as the wheelset moves laterally and the rolling radius difference increases, the rail roll is dynamically driven by the applied contact forces and the effective torques with respect to the rail fixing reaction point [29,30]. The longitudinal shift for a contact point that can occur, e.g. due to wheelset yaw, typically for a contact between the flange and the gauge corner, is illustrated in Figure 4(b) and can result in additional yaw torques on the axle. Here, γ w represents the wheelset yaw angle and x cp the longitudinal shift in contact point location with respect to the wheel centre. Related to the aspect of the longitudinal shift in contact point due to yaw is whether the change in the effective contact geometry   of the wheel profile as projected on to the contact plane is accounted for in the contact search. An illustration of how the projected wheel profile changes as a function of yaw angle is presented in Figure 5. It can be noted that the largest wheelset yaw angle observed in the Benchmark results is about 10 mrad where the influence of this effect is very small as observed in the figure. The conclusion that this amount of yaw has a small influence on the contact conditions is further supported by the findings in [31]. The functionality of the w/r coupling contact searches in the different software are listed in Table 16 for wheel and rail roll and for the influence of wheelset yaw in Table 17.

Track geometry implementation
A particular feature of S&C is a nominally instantaneous change from tangent track to a constant radius curve in the diverging route. Table 18 lists the different modelling options used by software participants. The variability between software users is reported in Table 19. All software are able to model a short transition distance.    Table 20 lists the numerical integration method used by software participants and settings where applicable. Simpack and VI-Rail users used the same integrators as the software participants.

Results
The results chapter is organised as follows. First, results are presented for the baseline comparison Run 9, followed by results from the switch panel (Runs 1-4) and the crossing panel (Runs 5-8), respectively. The selected results are chosen to illustrate the following aspects 1. Lead axle kinematic motion 2. Lead axle w/r coupling general forces 3. Lead axle w/r contact details The results are evaluated focusing on the aspects covered in the method statements presented previously, and key aspects of the dynamic w/r interaction in S&C, i.e. the transitions from stock rail to switch rail in the switch panel and from the wing rail to the crossing nose in the crossing panel.
The results are generally presented by discussing separately the behaviour in the through and the diverging directions, and also comparing the different behaviour observed between the shorter 56E1-245 m-1:9.25 and the longer 60E1-760m-1:15 switch and crossing cases presented in Table 1.
Results presented are comparing each of the software participants while the results from user variability, for Simpack and VI-Rail, are treated as a separate section at the end of this chapter. This separation is done to improve the visibility of the plots by keeping one result set per software as the user variability remains generally lower in comparison. The list of software compared in this section therefore includes: NUCARS-WNT, SDITT, GENSYS, MUBODyn, VI-Rail, Simpack, VOCO, MEDYNA and Vampire (VDG 1 ).
Note: A slash sign (/) is used throughout the results description to separate values of 56E1 cases (Runs 1,2,5,6 on the left of the sign) from those of 60E1 cases (Runs 3,4,7,8 on the right of the sign), unless specified otherwise. Plots showing the software comparison are presented with a unique line type and colour for each participant software. For the contact location and contact angles, where a drastic change of the position exist, lines are replaced with dots and unique markers are used to better identify each result set. Where contact results are concerned, there is no differentiation between contact points 1, 2 and 3 and where two or three co-exist, they are all represented in the same way.

Baseline model comparison (Run 9)
A baseline scenario is added to verify the correctness of the vehicle modelling and the output reference system and produced data from each participant, using the conditions of Run 2 but removing the variable rail, effectively simulating a non-transitioned, noncompensated tight curve. In this run, the expectation is to get an initial disturbance at s-distance equal zero, followed by a steadying of the output as the vehicle stabilises into the curve. The selected results from Run 9 are presented in Figures 6-8, focusing on the kinematic output, steering forces and vehicle acceleration response. From the results, it can be noted that: • The sudden curve entry with no (or negligible) transition (cf. Table 18) leads to a rapid lateral displacement offset of the leading axle ( Figure 6(a)). One participant (MEDYNA) shows a slightly early transition with respect to s = 0, indicating a slight error in the implementation of the curve definition and/or its position. After about 2 m, the quasistatic curving lateral offset position for the leading wheelset is obtained within a range of +7 to +7.7 mm (±5.5% difference) for all participants. • All but one result set (MEDYNA) show steady-state upward absolute vertical wheel displacements at around -0.8 mm (Figure 6(b)), which indicates that they all account for the w/r kinematic motion properly and that the track stiffness has been incorporated correctly. • Quasi-static curving forces (Q and Y-forces not plotted here) are generally similar indicating that the vehicle has been modelled correctly. However, looking at the longitudinal X-force (Figure 7(a)) additional variations are visible both in peak force at 2 m and in the steady-state value which are ranging between ±4.5 and ±9.5 kN. MEDYNA longitudinal creep forces build up earlier as already observed with the axle motion. MUBODyn, which is the only software using the Polach model, gives the lowest steering forces, meaning that lower accuracy is observed for this model at higher creepages [37]. • Car body lateral acceleration results show peak levels around -0.8 m/s 2 at about 10 m into the curve (peak slightly outside the range plotted here). MUBODyn, which models  the vehicle suspension slightly differently from what is specified, shows a slightly lower and delayed response to the curving. • A closer examination on the contact location (Figure 8(a), −35 mm on the y-axis being the rail gauge face and zero the top of the rail) reveals that all software predict a move of the point contact toward the gauge corner. There is about 2.5 mm differences in the range of predicted lateral contact coordinate which depend on contact detection and definition from each of the software. Simpack shows an intermittent 2nd point contact on the edge of the rail gauge face, while Vampire (VDG) and SDITT predict two-point contact throughout the curve. This could indicate the contact condition is at the boundary of two-point contact. One code predicts intermittent two-point contact, two codes predict continuous two-point contact, while the remaining code predict two-point contact conditions are not yet established. Looking at the contact patch size in Figure 8(b), the transfer of load from the rail crown to the gauge corner is identified at about x = 1.8 m with a simultaneous drop (crown) and increase (gauge corner) of contact patch size. In the curve, the secondary contact points are comparatively small because of the smaller transverse radius at the gauge corner and the resulting higher ellipticity. The predicted main contact area ranges from 63 to 130 mm 2 (at s = 4 m) supporting the previous comment on the rather large variability in contact prediction for this simple curving case with a constant rail profile. Yet, more variability can be expected from varying rail cases that follow.

Switch panel (Runs 1-4)
The simulation in the switch panels is characterised by a transfer of load between the stock rail and the switch rail with a momentary shared contact between the two, which is lasting longer in the diverging direction. This leads to deviation of the axle motion in the area of the switch toe near the start of the switch, resulting in dynamic amplification forces. In the diverging direction, these forces are higher, as earlier contact with the switch rail occurs, thus initiating the sudden change in the travel direction of the vehicle forced into the nontransitioned and non-compensated curve. This section describes the results comparison grouping them by through and then the diverging route.

Through route (Runs 1 and 3)
The selected results from Runs 1 and 3 are presented in Figures 9 and 10 focusing on the vertical motion, contact effect and vertical loads. The following are observed: • Despite a through route direction of movement, the entry into the switches leads to the fairly large lateral movement (Figure 9(a,b)) and oscillation of the lead axle toward the diverging direction initially, with the longer 60E1 switch leading to a larger offset of +7 vs +4 mm for the 56E1 one (peak not seen in Figure 10(a) due to selected x-axis limits). This corresponds to a longer length spent contacting the diverging stock rail (see contact coordinate plot in Figure 9(c,d)) and thus leading to a greater influence of the rolling radius difference affecting steering. The kinematic oscillation wavelength appears to match for all software with about 10/6.4 m (crest-to-crest semi-wavelength), and mainly depends on the w/r kinematics (equivalent conicity) and wheelset yaw inertia. • The contact coordinate (Figure 9(c,d)) clearly shows the diverging contact point on the stock rail and the sudden jump onto the top of the switch rail at about 3.5/8 m. All software show good agreement with a few millimetre difference in the contact position and a few more variations on the first contact with the switch in Run 3 between 8 and 10 m, as a split contact appears that then converges to the gauge corner position as the switch rail shape builds up to a full nominal rail section. • The absolute vertical wheel motion shows a lowering of the wheel of about 1.8/1 mm, with an ensuing rapid rising motion onto the switch rail. For Run 3, the wheel further rises by more than 1 mm above the original level in correspondence with a large lateral offset and two-point contact, before returning to zero vertical height. All results sets are in good agreement with the exception of Vampire (VDG) for which the wheel remains levelled on the stock rail and MEDYNA for which the wheel remains lower after contacting the switch rail. On further discussion with VDG, the values shown are based on a wheel lift output variable. This output is a relative value working properly for a single contact geometry (e.g. Run 9, Figure 6b) but incorrectly for cases where contact geometry is interpolated between multiple contact data files. • Vertical Q-forces (Figure 10(c,d) zoomed-in on 3-5 m) are characterised by a load transfer from the stock rail onto the switch rail at about 3.5/7.5 m. Following this, there is ensuing dynamic load amplification of higher frequency than otherwise seen from curving. The peak force generally occurs about 0.3/0.65 after the first contact with the switch rail and all software show a reasonable match with peak force value in the range 65-68.6/60-64.5 kN. This is excluding results from NUCARS predicting additional transient effect before load transfer and a peak force up to 83/75 kN, as well as Vampire (VDG) with 87 kN peak load for Run 1, out of phase reaction force and delayed peak response for Run 3. These higher forces are the expected result of the rigid connection between the two rails in the rigid profile definition adopted for these software (see Table 5).

Diverging route (Runs 2 and 4)
The selected results from Runs 2 and 4 are presented in Figures 11-13.   From the diverging results, it can be noted that: • Contact forces in the diverging route are characterised by a sustained two-point contact on both switch and stock rails, starting at about 1.9/4 m, until the load transfers fully onto the switch rail at about 3.5/7.5 m. The first pick up point can vary slightly depending on the results sets mainly due to curves set-up as already observed for Run 9. The vertical load sharing is roughly 2/3 on the stock rail and 1/3 on the switch rail (Q-forces not plotted). • The lateral Y-forces (Figure 11(a,b)) are characterised by a shared lateral load between the switch and the stock rail that lasts for about 2/5 m. During this two-point contact, the forces acting on the switch and stock rails are in the opposite direction, forcing the two rails together. The larger force is applied onto the switch rail (about 30-40 kN positive), which takes the contact nearer the flange of the wheel with a higher contact angle and larger lateral component to the load, while the force on the stock rail is about 10 kN negative. Quasi-static force levels are reasonably similar for all results. However, a number of additional transient effects occur especially after the switch rail returns to the standard rail profile for Run 2, which is indicative of rail profile interpolation effects at the back end of the switch. In the case of Run 4, the steering forces are reducing significantly after the load has fully transferred onto the switch rail due to the larger curve radius. • On contact patch location (Figure 11(c,d)), in the diverging route there is an earlier contact with the switch rail, which is shared with the stock rail for some distance ( ∼ 1.5/4 m). In both runs, the two-point contact continues onto the switch rail after the stock rail contact disappears. As previously indicated for the through route, the main differences observed between software are the very transient transfer of contact from the stock rail to the crown of the switch rail for Run 4 (s > 8 m), a few contact jumps are predicted by NUCARS on the stock rail for Run 2 and two-point contact are predicted by Vampire (VDG) on entry to the switch for Run 4. • On the contact patch size (Figure 12(a,b) showing contacts on both stock and switch rail), the trend is the same for all results and qualitatively explained: contact patch increase as the load transfers back to a one-point contact of a normal size (around 50-60 mm 2 ) after another 0.5/2 m. At this point, there remains a substantial variation between software (range 60-135/35-55 mm 2 ) that indicate different approaches in contact detection and estimation in curve. These discrepancies are inherently due to the various contact methods employed (e.g. Hertzian/elliptical vs non-Hertzian/non-elliptical, or conversion made to equivalent Hertzian while outputting value). • The contact angles (Figure 12(c,d), only showing switch rail contact) are fairly consistent for all software. However, larger variations existing in the transition area beyond the variable switch when returning to steady-state curving conditions, for the same reasons explained in the paragraph (4). • The normal contact forces (Figure 13(a,b), only showing switch rail contact) are generally consistent for the first initial contact on the switch rail (s = 1.7-3.6/3.6-8) with the exception of MEDYNA (higher force), due to reasons already highlighted with Run 9. Vampire (VDG) shows differing normal force with a sustained two-point contact passed 5 m on Run 2, while SDITT predicts a later return to single point contact at about 5.5 m. Normal forces in the area of two-point contact on the switch rail show higher variability, in correspondence with observation already made with contact patch size and contact angles. • The longitudinal creep force (Figure 13(c,d), only showing switch rail contact) applied by the wheel flange onto the switch rail is negative -opposite to the direction of travel, while the force generated by the wheel tread contact on the switch rail crown during two-point contact is positive, i.e. in the direction of travel. In the area of two-point contact and on the switch, the forces range between negative 6-13/8-19 kN, while on the tread, it ranges between positive 4.5-9.5/4-12 kN, which is a significant variation between software.

Crossing panel (Runs 5-8)
The simulation in the crossing panels are characterised by a very rapid transfer of load between the wing rail onto the crossing nose rail, associated with a large vertical dynamic impact-like load. This is due to the wheel dropping its height as the wing rail diverges and the contact with the wheel rapidly moves toward the field side leading to an effective drop of the wheel centre of mass due to its general conical shape. After the wheel transfers onto the crossing nose, it then rapidly rises up onto the nose topping (ramped shape of the nose), and this corresponds to a rapid change of vertical direction of motion of the wheel that generates high inertial vertical force. Other lateral and steering effects are also ensuing due to rapid changes in rolling radius difference.
In the diverging route, the dynamics is further complicated by the action of the check rail on the opposite side to the crossing, which maintains/forces the wheelset within a certain lateral offset position so that it does not make interference contact with the crossing nose. Similar vertical dynamic behaviour occurs, with some variations due to slightly offset contact path on both crossing rails and wheel, as well as a redistribution of forces between the opposite rail, check rail and crossing rail.
This section describes the key results comparison grouping them by through and then the diverging direction. Note that results from Vampire (VDG) participant are not supplied for Runs 5-8.

Through (Runs 5 and 7)
The selected results from Runs 5 and 7 are presented in Figures 14 and 15.
From the results, it can be noted that: • On wheel vertical kinematics (Figure 14(a,b)), two different patterns emerge from the 56E1 and the 60E1 crossings. For Run 7, the typical downward pointing triangular shape is clearly visible, with a drop of about 3 mm at the lowest point ( ∼ s = 0.5 m). For Run 5, the rapid change of shape of the wing rail from nominal rail, leads to an initial wheel upward motion of ∼ 0.75 mm which returns to its initial position at the other end of the crossing (vee to nominal). In the middle section, there is also the typical inverse triangular shape, but it is much shorter, due to the wider crossing angle, with an overall drop of about 1 mm. Overall, this means the 56E1 crossing generates additional transient dynamics at the entry and exit of the crossing element in addition to the load transfer from wing to the nose. In terms of the comparison between software, for Run 5, all software capture the same general trend, while a few show some degree of dynamic oscillation (wavy shape) after entry on the leg end and after load transfer. For Run 7, they all show good agreement and a smooth wheel vertical kinematic motion over the variable rail shape, while NUCARS predicts small step-like variations in places. Note the predicted dip angle from all software are about 16.5/5.6 mrad, which is typical for these two types of crossing designs considering a nominal wheel. • The contact patch location (Figure 14(c,d)) can be seen to move from a fixed position on the crown of the nominal rail to a sudden shifted position in the field side of the wing rail on entry into the crossing. For Run 5, this change is very sudden as the wing rail shape builds very quickly and this corresponds to the initial upward motion of the wheel in Figure 14(a) as well as the initial transient loads predicted in Figure 15(a). For Run 7, the movement is more gradual as the contact steadily follows the path of the diverging wing rail. Past the IP point (s = 0) and nose tip, the contact transfers over onto the crossing nose, initially very near its gauge corner and then gradually moves to a more centred position on the crossing vee. For Run 5, there is a further shift of contact point toward the field side of the crossing vee, which is due to the unworn crossing shape, and for  Run 7, there is a more gradual 'U' shaped movement of the contact toward the field side (s = 1.5) and then back onto the crown of the rail, to the initial nominal rail position. All software are predicting the contact location in a similar fashion, with only a minor offset in the contact shift/movement. For NUCARS, the contact location appears to be slightly shifted further away toward the field side in general and more particularly on the wing rail. • The vertical Q-forces on the crossings are characterised by a number of high amplification dynamic loading/unloading in three key areas: (a) the entry from nominal rail to the wing rail characterised by a rapid change in contact position, (b) the load transfer from the wing rail to the nose and (c) the exit from the crossing back to the nominal rail. The predicted unfiltered force amplification can reach up to 6/3 times the nominal wheel load in those areas. The highest amplification occurring for Run 5 wider crossing angle for which the wheel dip angle motion is three times higher than Run 7 (see wheel kinematics). Passed s = 0, the wheel contact with the wing rail shows a reduction of load for Run 7 and a complete loss of contact for Run 5. In the case of Run 5, the subsequent impact load on the nose is the highest and a rebound of the wheel against the rail mass (second loss of contact) and further transient effects follow. For Run 7, there is also a further loss of contact force at about s = 1.8 m attributed to the lateral shift of the point contact toward the top of the crossing vee (see Figure 14(d)) and the effective rapid reduction of wheel radius at the contact point. In terms of software comparison, all of them capture the overall behaviour; however, the magnitude of dynamics loading, the exact pick-up point with the nose and some resulting delays are observed. NUCARS result shows a lower frequency response to load transfer on both runs. Force magnitude in both cases is lower than the rest of the software prediction. The MEDYNA results  also show a slightly anticipated contact transition related to the vehicle effect discussed for Run 9. • The Dynamic Impact Factor (DIF) at the crossing nose (derived from Figure 15(c,d)) can be further investigated to highlight the dynamic response characteristic of vertical dip angle irregularity input as described by Jenkins for rail joint [38] and further elaborated by others in the case of crossings [39][40][41]. From most results, we can see that the load transfer impact force is made up of a first higher frequency peak (P1) and an overall longer wavelength dynamic response (P2). The P1 peak appears roughly in the range s ≈ 0.28-0.31/0.425-0.500 m and the overall P2 duration occurs in the range s ≈ 0.28-0.40/0.42-0.70 m. The estimated values from the plot are tabulated in Table 21 and an average calculated, to produce the bar chart in Figure 16. P1 force is taken at the maximum value following load transfer, and P2, the maximum value after low-pass filtering with a cut-off frequency of 200 Hz. This shows a general good agreement between software for Run 7, and more variations for the harsher crossing of Run 5.

Diverging (Runs 6 and 8)
The selected results from Runs 6 and 8 are presented in Figures 17 and 18. From the results, it can be noted that: • In the diverging direction, the lead axle is forced away from its offset curving position back toward the track centreline, under the influence of the check rail. This result in high lateral Y-forces throughout the crossing panels as plotted in Figure 17 (Figure 17(c,d)) shows a similar trend to the ones in the through routes, except that the contact on the crossing occurs much nearer the rail gauge corner because of the axle offset position nearer the crossing and the larger contact jump toward the field side of the rail on the diverging wing rail just before the load transfer occurs. This leads to the lateral transient loads explained earlier. Compared to most, there are slight offset with NUCARS and additional lateral shift, for example, Run 6 just before load transfer. • It is worth noting that the vertical Q-forces (not plotted) are partly taken up by the contact with the check rail (around 1/4 or 1/5) due to the 10/4°contact angle from vertical (see below). • The contact angle (Figure 18(a,b)) on the crossing can be seen to vary greatly depending on the contact proximity to the gauge corner or the top of the wing rail. As for the contact angle on the check rail, the values are all fairly steady at about 1400/1500 mrad, 80/86°. Note that to match the over results, the check rail contact from MEDYNA are divided by factor 2.5 and 10 for Runs 6 and 8, respectively, possibly explained by an extension of the inner wheel flange by the user having a significant effect on the projected wheel contour with larger axle yaw angle. The large variation in contact prediction from NUCARS also translate in a large variation of contact angles with the check rail which are visible across the figure as series of dots oscillating between zero and the check rail contact angle, in particular for Run 8. • The longitudinal creep forces on the check rail are following the same trend and overall shape as already described for the lateral Y force with 6-10 / 3-10 kN steady-state load at s = 0 m. MEDYNA shows higher forces particularly on the entry flare due to the previously observed differing vehicle curving behaviour. VOCO shows larger forces in the horizontal section of the check rail for Run 6, due to the simplified approach to modelling the check rail contact. VI-Rail forces are also slightly higher than the rest. MEDYNA and VI-Rail show the discontinuous force described previously for Run 8 due to poor capture of the check rail shape. NUCARS shows the highest transient peak forces and unstable forces already described focused around the area of change of path or shape of the check rail.

On user variation
User variability is presented in Figures 19 and 20 for Simpack and VI-rail users, respectively. Selected outputs are presented based on Run 8, given than the most relevant variabilities are observed on the most challenging cases of crossing (highly variable rail shape) and in the diverging route (interaction with the check rail). Generally, the variation based on user input is very small and all the initial obvious errors have been removed after some iterative discussion between users as a whole and within each software group.

Simpack user variation
The following observations are made: • PROSE indicate a slightly higher (+15%) vehicle lateral acceleration compared to the average of the rest of the users. This output is generated using an accelerometer element (167) which includes gravity effects and generates an offset of 0.13 m/s 2 on average, due to a car body roll angle of 13 mrad. • All wheels vertical kinematic motions are matching very well (0.19 mm range at the lowest point of ∼ 3.7 mm, therefore ±2.5%).
• The contact patch size prediction on the check rail and crossing wing rail show noticeable variations. In particular, Chalmers results using a different contact type for the check rail detects a hugely increased contact patch size. • Contact forces on the check rail show impact transient for most user, except Chalmers (smooth) using a different contact element with the check rail. UoB predicts a much earlier contact with the check rail. • Most users predict discontinuity in loads on the first and last profile of the crossing (−2 and 1.2 m) but with varying magnitude, indicating an influence of user handling of the input data for variable rails and use of section breaks. • PROSE shows high-frequency oscillation of small magnitude in the contact forces, due to the use of a different contact model (discrete elastic as opposed to equivalent elastic, combined with high damping settings).

VI-Rail user variation
The MEDYNA software has been added to this comparison because the contact algorithm used in VI-Rail is the one developed by ArgeCare for MEDYNA, although the rest of the multibody dynamics software implementation and solver in particular are different.
The following observations are made: • MEDYNA indicates a slightly higher (+12.5% average) vehicle lateral acceleration compared to the average of the rest of the users, indicating differences in the vehicle model, speed or curve definition. • All wheels vertical kinematic motions are matching very well (0.16 mm range at the lowest point of ∼ 3.7 mm, therefore ±2.1%).  • The contact patch size prediction on the check rail and crossing wing rail show noticeable variations. • All MEDYNA/VI-Rail users predict distorted contact with the check rail, due to the same algorithm used to capture the check rail surface and the errors introduced in the area of the flare ending because too few input cross-sections are used. • All users predict discontinuity in loads on the first and last profile of the crossing (−2 and 1.2 m) but with varying magnitude, indicating an influence of user handling of the input data for variable rails and use of section breaks. • GMUNIFI shows high transient contact forces in a few areas such as the exit from the crossing for Runs 7 and 8 at 1.25 m due to a lack of quality on the longitudinal interpolation at the end of the crossing (last section break definition).

Main observations from the method statement compilation
The following main observations are noted from the method statement compilation section.
• Among the software in the Benchmark, there are three methods presented for representing individual cross-sections in software; (a) spline fit, (b) discrete data points and (c) a fitted chain of circle segments. For the software using a spline fit to represent the rail cross-sections, participants commonly applied different levels of smoothing for more consistent profile curvatures.
• Among the software in the Benchmark, there are four principal methods to represent rail geometry between given cross-sections in time domain simulations: (a) perform geometry interpolation between the provided cross-sections online during simulations to obtain the profile at a given wheel position; (b) solve the w/r contact problem for cross-sections in advance for a range of contact conditions, tabulate the results, and then interpolate between the contact tables in time domain simulations; (c) tabulate the parameters needed for contact calculations for each cross-section in advance, and then interpolate between the tables in time domain simulations and use as input in solving the w/r contact problem; (d) implement the geometry as a series of piecewise constant sections. Different forms of profile modifications and sectioning to create interpolation breaks are created by participants to ensure a smooth and physical representation of the rail running surfaces, either in the form of geometry, contact table or parameter table interpolation, in time domain simulations. • Most software implemented the specified co-running track model with minor variations, but some simplifications are present, e.g. massless rails, merged rail bodies (e.g. switch rail and stock rail treated as one body) and rail bushing connections to ground instead of track mass as not all software readily support double track mass-connected rail bodies on the same side of the track. • A range of contact models are employed to model the normal contact. From linear spring stiffnesses and equivalent Hertzian contact patches to models utilising a discretisation of the contact patch. For the tangential contact, different variations of FASTSIM are the predominant choice. A few participants used a simplified linear contact stiffness and assumed full slip for the tangential contact with the check rail, due to the nearvertical face. All participants used standard steel properties in the contact modelling or a linear contact stiffness providing a similar contact stiffness. There is a greater variation in the contact damping settings used. • The w/r contact search functionality varies between software, but with the rail roll DOF locked to the track mass in the Benchmark track model and the relatively small wheelset yaw angles observed in the results, it should not matter significantly for the results in the Benchmark whether the contact search accounts for rail roll or a change in the wheel profile contact geometry is accounted for in the contact search. Whether wheelset roll and longitudinal shifts in contact point locations are accounted for can still be significant. • All participants used a short transition curve length to model the step change between tangent track and curve for the switch simulations in the diverging route. • A variety of numerical methods have been employed for time integration.

Main observation on the results comparison
The main conclusions on the ability of various software to represent the physical interaction of railway vehicle with switches and crossings are: • Most participants implemented the Benchmark exercise appropriately; however, small errors are still present from either the vehicle definition or the way in which the track curvature is implemented. This is made obvious with Run 9 results and has some impact for all the turnout cases with slightly offset or delayed attitude of the vehicle/axle to the curving situation, i.e. small deviation in lateral position and yaw give rise to variation in contact conditions and dynamic loading. • There is generally a better match in the switch cases than the crossing cases, which are more challenging to implement because of the more rapid and wider variation in rail shape, and therefore the highest variability seems to come from dealing with the effective 3D rail profiles representation in the various software. • Where the switch and stock rails are considered as one body, the load transfer occurs much more suddenly and the dynamic amplification load increases. • Another high variability comes from the contact detection and parameter determination where rapid changes occur and where high contact angles are present. Perhaps surprisingly, even Run 9 provides different contact detection and split of contact showing that tight curves analysis with a single rail profile can also be a subject of further investigations. • Large differences can be present from approaches using pre-calculated contact tables as opposed to the online calculation of the rail shapes to feed into the contact algorithm.
Where linear interpolation is used in between these tabular data, discrete 'step-like' results can be obtained particularly for the crossing cases. The precomputed tables of parameters (e.g. VOCO) with separate contact interpolation for each rail body provide steady and reliable results, however. • For the crossing cases, additional transient loads to those expected in the wheel transfer area (wing to nose) are apparent, either due to rapidly varying contact conditions (lateral shift or 'jumps') or due to the introduction of discontinuity from the 3D interpolation of the rail shape. • The definition of the check rail is very important for diverging cases in crossing panels.
Simplifying this interface with spring-dashpot representation to constrain the axle at the correct lateral position seems to work well in comparison to using a full w/r contact. In the latter case, using too few 2D profiles sections in areas of rapid change of shape (flare) leads to distortion of the 3D interpolated check rail shape with software MEDYNA and VI-Rail. NUCARS-WNT presents unstable contact and loading when dealing with varying shaped check rail. • The main conclusions from the user comparison are (note these conclusions are concerning only two of the commercial software participating in the Benchmark, as they are the only ones having multiple users): (1) Generally, the user variability is low, each commercial software providing a solid platform for the user to ensure the proximity of results within a few per cent, provided there are no errors in implementation (vehicle modelling or track modelling). (2) The two commercial software evaluated by multiple users provide a very similar approach to dealing with variable rail surface (i.e. 3D interpolation based on input 2D cross-sections and in varying number of joined sections) and contact algorithm. The only differences come from their usability and feedback in terms of handling the input data for the variable rail definition. Most discrepancies between users come from the details of where and how section break are defined to ensure suitable continuity of the rail surface for the contact algorithm. Also, the profile smoothing available to Simpack users can have an influence on the contact output, while this is not an available option to VI-Rail users.

Final statement
A benchmark of railway multibody dynamics software ability to model vehicle interaction with switches and crossings has been carried out for the first time, comparing all major commercially available software and a few independent codes developed within research institutes. Notwithstanding initial variation principally due to errors of interpretation and implementation due to the complexity of the task at hand, the final results show that all software offer a reliable and efficient way to understand the wheel/rail kinematic and the dynamics forces between the wheels and the track elements. The highest challenges are found when modelling a combination of multiple interacting rails (e.g. check-stock or switch-stock), high longitudinal variation in rail shapes (crossings) and high lateral steering forces (i.e. diverging cases in tight radius). In those cases, the codes able to account for the exact relative motion of each wheels with respect to each rails independently are the most apt. The most significant variations between software are found in the contact prediction (detection: position, size and angle) and this has some influence on the detailed contact tangential (creep) and normal forces, which are expected to lead to varying levels of damage if those are used for further damage estimation.
Finally, the user variability is found to be very small, with the most time-consuming and error prone being the task of handling the input data for the variable rails definition. For this task, it is believed that all software could benefit from some improvement to assist the user and ensure higher reliability and efficiency of the process.
It is hoped that this benchmark will serve the railway community to improve its practices and understanding of dynamics in switches and crossings, provide new generations with a reliable starting point for their study and eventually help improve S&C design and their robustness. Note 1. VDG is not a representative of Vampire but the sole user of the software in this Benchmark. Table 22 gives the full list of contributors that have actively supported the benchmark including those that cannot be named author of this paper due to space restriction. Source of funding where appropriate is acknowledged in italic.

Disclosure statement
No potential conflict of interest was reported by the author(s).