Aquib Qazi

Keywords: railway dynamics, wheel-rail contact, contact mechanics, rolling contact, switches and crossings, multibody

The wheel-rail contact constitutes the link between the rolling stock and the track dynamics. The variety of different applications within the railway network necessitates different levels of modelling. For simulations carried out over railway tracks that are several kilometres long, approximate wheelrail contact models are generally sufficient. Conversely, some applications require more sophisticated models, such as in the estimation of wear in urban rail networks, or during the crossing of singularities of the track such as in switches and crossings (S&C).

Three criteria may be used to broadly classify the different wheel-rail contact modelling approaches: the assumptions about the location of the contact, the calculation of the normal forces (from the Hertzian ellipse to non-Hertzian contact), and the calculation of tangent forces (from analytical expressions to Kalker's complete theory). Contact modelling is an active field of research, especially in railway dynamics, where certain case studies may require assessing the three-dimensional transient response of rolling stock over long distances, and simultaneously solving the contact problem for each wheel up to 100000 times per second. A compromise is, therefore, often required between the computational efficiency and the precision of the modelling techniques. Moreover, the implementation of wheel-rail contact methods in multibody systems (MBS) software requires efficient and reliable methods of coupling the wheel-rail contact solution and the transient dynamics of the bodies closest to the contact interface (wheelset and track). The load transfers between bodies must also be handled with care.

The aim of this thesis is to describe, test, and extend the capabilities in terms of wheel-rail contact modelling of the MBS software VOCO, the railway dynamics code developed at Université Gustave Eiffel. Currently, VOCO allows real-time simulations, following two approaches, one Hertzian, and a non-Hertzian method considering non-elliptical contact patches. The present work proposes a new semi-analytical boundary element (BE) method i for determining the wheel-rail contact zone and the associated normal stress distribution. The contact pressure is assumed beforehand to take the form of an elliptical distribution in the rolling direction and evaluated using the BE method. The number of system unknowns is reduced via the discretisation strategy, and using a semi-analytical methodology enables fast computation speeds, an indispensable requirement in railway dynamics. Moreover, comparisons with other methods commonly used in MBS codes show a significant improvement in the normal contact results using Kalker's CONTACT software as the reference.

To address the tangential contact problem, an extension of Kalker's FASTSIM algorithm is presented. The new algorithm permits the handling of non-Hertzian contact through a local approach. Several aspects of the tangential contact modelling are investigated, the proposed algorithm being validated through a design of experiments. The new normal contact method is combined with the tangential contact algorithm, which then represents a new rolling contact model for the whole contact problem and provides the possibility of using a more detailed approach in the context of multibody dynamics.

Parallel to the development of new contact methods, a comparison of the contact modelling approaches used in different commercial and academic MBS codes has been made within the scope of an international benchmark on S&C. Following the benchmark, a collaborative study has been initiated to assess the wheel-rail contact results in detail, emphasising some delicate configurations which may occur in S&C, such as the conformal contact, the sharp-edge contact and impact loads. Future works should consist of the implementation of the new methods developed in this study, in VOCO.

Résumé

Le contact roue-rail constitue la liaison entre le matériel roulant et la dynamique de voie. La variété des problématiques ferroviaires nécessite différents niveaux de modélisation. Pour les simulations effectuées sur des voies de plusieurs kilomètres de long, des modèles approximatifs de contact sont généralement suffisants. À l'inverse, certaines applications nécessitent des modèles sophistiqués, comme pour l'estimation de l'usure dans les réseaux ferroviaires urbains, ou lors du franchissement de singularités de la voie telles que les appareils de voie (ADV).

Trois critères peuvent être utilisés pour classer les différentes approches de modélisation du contact roue-rail : les hypothèses sur la localisation du contact, le calcul des forces normales (de l'ellipse Hertzienne au contact non-Hertzien), et le calcul des forces tangentes (des expressions analytiques à la théorie complète de Kalker). La modélisation du contact est un domaine de recherche actif, en particulier en dynamique ferroviaire, où certains cas d'études peuvent exiger le calcul de la réponse transitoire d'un matériel roulant sur de longues distances, et simultanément la résolution du contact pour chaque roue jusqu'à 100000 fois par seconde. Un compromis est donc nécessaire entre efficacité en temps de calcul et précision des résultats. De plus, l'implémentation des méthodes de contact dans les logiciels multicorps (MBS) nécessite des méthodes efficaces et fiables de couplage entre la résolution du contact roue-rail et la dynamique transitoire des corps les plus proches de l'interface de contact (essieu et voie). Les transferts de charge entre corps doivent également être traités avec soin.

L'objectif de cette thèse est de décrire, tester et d'étendre les capacités du logiciel MBS VOCO, le code de dynamique ferroviaire développé à l'Université Gustave Eiffel, en termes de modélisation du contact rouerail. Actuellement, VOCO permet des simulations en temps réel, suivant deux approches, une Hertzienne et une non-Hertzienne prenant en compte des empreintes de contact non-elliptiques. Le présent travail propose une nouvelle méthode semi-analytique par éléments de frontière (BE) pour iv déterminer la zone de contact roue-rail et la distribution de contraintes normales associée. La pression de contact est supposée à priori elliptique dans la direction du roulement et est évaluée à l'aide de la méthode BE. Le nombre d'inconnues du système est réduit et l'utilisation d'une méthodologie semi-analytique permet de réduire le temps de calcul. De plus, les comparaisons avec d'autres méthodes couramment utilisées dans les codes MBS montrent une amélioration significative des résultats en utilisant le logiciel CONTACT de Kalker comme référence.

Pour traiter le problème du contact tangentiel, une extension de l'algorithme FASTSIM de Kalker est présentée. Le nouvel algorithme permet de traiter le contact non-Hertzien par une approche locale. Plusieurs critères à satisfaire dans la modélisation du contact tangentiel sont considérés et l'algorithme proposé est validé par un plan d'expériences. La nouvelle méthode de contact normal est couplée à l'algorithme de contact tangentiel, ce qui constitue un nouveau modèle de contact roulant et offre la perspective d'utiliser une approche plus détaillée dans les codes de dynamique ferroviaire.

Parallèlement au développement de nouvelles méthodes de contact, une comparaison des approches utilisées dans différents codes MBS commerciaux et académiques est effectuée dans le cadre d'un benchmark international sur les ADV. À la suite de ce benchmark, une étude collaborative est initiée pour évaluer les résultats de contact roue-rail, en mettant l'accent sur certaines configurations délicates qui peuvent apparaitre dans les ADV, tels que le contact conforme, le contact à bord vif et les impacts. Les travaux futurs devraient comprendre l'implémentation dans VOCO des nouvelles méthodes développées dans cette étude.

Mot clés : dynamique ferroviaire; contact roue-rail; mécanique du contact; contact roulant; appareils de voie; multicorps. v Preface I would have never imagined when I arrived in France five or so years ago that I would be here, presenting my research in front of a jury. It has been quite a journey, even more so when we think about how our normal life came to an abrupt pause for several months due to a global pandemic (!!!), and yet somehow, this only seems to be the beginning.

It's funny how we add the acknowledgements and the thank-you notes at the beginning of any thesis, even though they are pretty much the conclusion of three years of hard work and pleasure shared with a large number of people, people without whom none of this of course would have been possible. The invaluable backing and encouragement I have received has been instrumental in this personal (and collective) endeavour to make a small difference in the world of railway engineering. It is thus only natural to acknowledge all those who have contributed directly or indirectly to the advancement of my work.

Honoré, you have supported me from day one of this project and have always pushed me to be better, to be thorough and clear in presenting my ideas. Thank you for having your office door open for in-depth discussions, for your advice (especially the one to pay closer attention to details!) and your guidance over the past three years. Hugues, I will treasure the chats we had about everything and more. Thank you for your insights, for your honest feedback, and all your invitations to go hiking which I wish I could have taken up more often. Michel, you have been the best support all this time whether it was professionally or personally (and even in administrative matters!). Supervising a thesis on a daily basis is never easy, thank you for the time you took to explain things that weren't always easy to explain. Whatever little I have learnt about railway dynamics, it is because of you, and for that I will forever be grateful. Thank you for opening up this wonderful field to me and I hope we will have opportunities to continue working together in the future.

To all my colleagues at COSYS, GRETTIA, and NAVIER: thank vii you for making coming into work such a fun and pleasant experience. Emmanuel, thank you especially for all your honest comments and feedback, for accepting the role of an "unofficial" supervisor. My thanks also goes out to all the wonderful people at ESI Group and the ANRT, who supported and financed this work with interest: Anne, Mustapha, Cédric, Paco... A special mention is also deserved for the industrial partners from SNCF Reseau, Vossloh, Alstom, and Vibratec who participated in the various follow-up committees and always provided useful suggestions. I want to thank Professor Yves Renard and Professor Frédéric Lebon for agreeing to be the referees for this work, for their detailed reports and their insights, and finally for the gainful discussions and questions during the thesis defense. I want to thank Ms. Vanessa Lleras for agreeing to be a part of the jury, for taking out the time to go through the manuscript, and for her interesting questions. I would finally like to give my thanks to Professor Yann Bezin for his perceptive comments and the rewarding exchanges, as well as for organising the benchmark on switches and crossings, which provided interesting modelling challenges and fruitful results for the entire railway dynamics community, including some of the work presented in this manuscript.

To my fellow Phd students, Mohamed, Maria, Younes, Demeng: you all have become close friends over this journey we shared. Mohamed, thanks for all the shopping deals and thieb/mafé recipes you shared, for being my office sparring partner for all of the past three years. Thank you guys for making it so enjoyable, time really does fly when you're having fun. To the priceless connections I've made since I landed on these once-foreign shores, Faizan, Mallory, Aniket, Sarah, Aakifah and so many more: thank you for being there all those times when I needed it the most. To my oldest and dearest friends, Faraz, Rushabh, Farheen, Kassim: we've grown up a lot but you all are still the best. A small note of love for someone who will definitely not be reading this (and to be fair, not reading anything), my cat Jeïna, who tolerated all the troublesome work-from-home hours imposed on her nap times! Finally, I reserve my biggest share of gratitude for my family. To all my uncles, aunts, cousins, nephews and nieces spread over the four corners of the world, my sincerest thanks and so much love. Mom, Dad, finally a doctor in the family! (Even if it isn't the kind you imagined all those years ago!) Thank you both for your constant and immense support, for your invaluable prayers, for believing in me even when I did not believe in myself. None of this could have been ever possible without your love and encouragement. Everything I am today, I owe it to the upbringing you gifted us with. Much viii love goes out to my sister who has always taken care of me, since I was a little kid on my way for the first day of school, and who, as fate would have it, was also the only one representing my family all the way in France on the "last" day. To my younger brother, with whom I share the pain of supporting a football team that hasn't been doing so well (cc Manchester United), the love of stand-up comedy, anime, and banter, and who took out the time to proofread the different chapters and articles even though he, in his own words, "did not understand a single thing". I love you guys!

Introduction

Rail is often seen as one of the most safe, innovative, and sustainable modes of transport, providing an attractive response to the continuing demand for a shift to greener alternatives. Many issues linked to the ever-increasing population density and urbanisation such as road traffic congestions or air and noise pollution can be addressed by a systematic change to rail for the transport of people and freight. With 2021 being declared as the European Year of Rail by the European Union, the onus is more than ever on the different players in the railway industry to deliver more competitive and resource-efficient solutions through research and innovation, replacing physical testing, which is often very expensive, with virtual prototyping and simulations to qualify and approve rolling stock and railway infrastructure. The challenge at hand in railway dynamics is to evaluate the dynamic response of the system in different running conditions, with a focus on the wheel-rail contact forces whose variations are mainly due to the changing geometry of the track.

or nonlinear springs and dampers. Bodies are usually rigid with 6 degrees of freedom (dof), but flexibility can be incorporated through modal stiffness and damping properties computed using FE tools as and when required. The specific design of rolling stocks leads to a distinction between the wheelset and the other bodies of the vehicle. The wheel may undergo high-frequency variations of wheel-rail forces, for instance, due to contact jumps between the tread and the flange or sudden changes of the rail geometry such as in turnouts. On the other hand, the other bodies reactions are low-pass filtered by the suspensions. An MBS modelling of the rolling stock is relevant in railway dynamics because most of its flexibility is ensured by the links between its bodies. The wheel-rail contact may be viewed as the link between the wheelset and the track infrastructure and is integral to all research related to vehicletrack interaction in railway dynamics. The contact forces govern the dynamic behaviour of the vehicles. The contact zone is also subjected to high levels of stress in small concentrated areas, leading to several damage phenomena. As a direct result of the numerous functions that must be ensured at the wheel-rail contact level (load transfer, guidance, brakingtraction), the modelling of wheel-rail contact is of utmost importance in MBS railway dynamics and must be handled with care.

The rolling contact model considered as the primary reference in this study is Kalker's variational theory, which is implemented in the CONTACT software. The contact problem in CONTACT is solved via the boundary element (BE) method. CONTACT has been developed extensively over the past years, broadening its range of applicability with some extensions such as falling friction or third body layer, which are out of the scope of the present work. Similar development efforts have also been made to enhance the program's speed versus accuracy trade-off. This has led to the implementation of CONTACT as an add-on in several railway MBS packages (SIMPACK, GENSYS, NUCARS and Universal Mechanism). The relatively high computational cost of CONTACT still inhibits MBS users from using it extensively in railway dynamics cases studies, and various approximate methods are instead implemented in most MBS software, including VOCO (Voitures en courbe), the explicit railway dynamics code developed at Université Gustave Eiffel over the past three decades.

An important hypothesis in Kalker's theory is the quasi-identity of the contacting bodies, which is linked to the identical material properties of the wheel and the rail. The contact is also assumed to be concentrated on a three-dimensional half-space, for which the surface response is known analytically through the theory of Boussinesq and Cerruti. The quasi-identity assumption results in the decoupling of the normal and tangent problems. Thus, wheel-rail contact modules implemented in MBS codes are a combination of a normal and a tangential contact method, that solve the whole rolling contact problem sequentially.

Chapter 1 starts with a brief overview of the different modelling strategies used in railway MBS, with a particular interest in the modelling of the wheel-rail contact. Some of the common contact methods implemented in railway dynamics codes are presented. As far as the normal contact is concerned, focus is made on methods which consider non-elliptical contact patches, as they provide a finer assessment of contact stresses when compared to the Hertzian methods where the patch is constrained to be elliptical. For the tangential contact, most methods, including the popularly used FASTSIM algorithm, are restricted to Hertzian contact patches, although a few non-elliptical variants exist in the literature.

Chapter 2 outlines the development of a new boundary element method for determining the wheel-rail contact zone and the normal stress distribution within it. A reduced BE formulation of the contact between two elastic half-spaces is proposed. Together with additional assumptions based on the Hertzian theory applied locally, the normal contact problem is solved following a semi-analytical methodology. The novel method is then validated using the results from CONTACT as a reference, with theoretical as well as wheel-rail profiles. A comparison of the obtained results and the required computational resources is also carried out versus some of the approaches presented in Chapter 1.

As the original FASTSIM algorithm is restricted to Hertzian contact patches, a strip-based approach is proposed to extend it to non-elliptical contact cases. Chapter 3 presents this local approach based on curvature properties in detail, which was first introduced by Ayasse & Chollet along with the semi-Hertzian normal contact method. A new variant is introduced using a weighing process, which is then shown to provide more consistent results through a design of experiments. The combination of this algorithm with the normal contact method presented in Chapter 2 represents a new solution for the whole rolling contact problem. This is shown to be promising, even if the contact location and the contact kinematics are supposed to be known beforehand in the cases tested here.

Chapter 4 deals with case studies in railway dynamics involving S&C within the context of an international benchmark of academic and commercial MBS codes. A collaborative study investigates the pros and cons of the different wheel-rail contact models when handling specific contact configurations that may occur in S&C simulations.

Chapter 1

Overview of wheel-rail contact models used with the multibody systems approach

Introduction

The last few decades have been characterised by exponential developments in the hardware capabilities of personal computers. Numerical modelling and simulation of complex mechanical systems have now become routine tasks that can be accomplished with relative ease. This is in stark contrast to only a few decades ago when computer systems were a novelty and highly specialised according to their intended applications. In general, the computational approach begins with the idealisation of a system of interest. This idealisation process introduces suitable hypotheses to ensure that the system of interest can be represented reliably to the desired level of accuracy, but at the same time remain simple enough to allow efficient and fast simulations. These two requirements are essentially contradictory, which further underlines the difficulty of modelling a real (non-idealised) physical system or phenomenon.

The model of the idealised system is expressed using mathematical equations that represent its physical behaviour, and the interactions between the different system parameters. Some parameter values, for instance, damping or friction coefficients, may remain unknown and must be introduced via statistical approaches. Solutions of the mathematical equations are sought using different numerical techniques and interpreted, with the modelling considered satisfactory if the numerical predictions agree with the experimental findings within an acceptable tolerance [START_REF] Popp | Ground vehicle dynamics[END_REF].

The various mechanical components of the system of interest may be modelled differently, either through continuous systems, finite element (FE) systems, or multibody systems (MBS) approaches. While there is some freedom of choice for the approach that is selected for a given application, the number of degrees of freedom (dof) characterising the behaviour of the physical model and the corresponding equations of motion do play an influencing role. For continuous systems, analytical solutions exist at any location but are restricted to simple geometries, for instance, simple beams, circular plates etc. Modelling with FE systems is very flexible as they do not present any restrictions concerning the geometry, which is usually represented through non-rigid elements like rods, beams, plates etc., with a finite number of dofs at the nodes. The behaviour of the system between the nodes is subsequently assessed through appropriate shape functions.

The MBS approach, on the other hand, consists of representing the kinematics of a complete system through a finite number of bodies or elements. These bodies or elements are interconnected with each other and with their surroundings through linear or nonlinear springs, viscous dampers, force-controlled actuators or a combination thereof. MBS bodies are usually rigid, but flexibility can be incorporated through modal stiffness and damping properties computed using FE tools as and when required. The equations of motion of an MBS are written as the following initial value problem:

Mÿ + C ẏ + Ky = F , (1.1)
where M is the lumped mass matrix (in the case of rigid bodies), K is the stiffness matrix, C is the damping matrix, y is the vector of the vehicle dofs, and F is the external force vector. For railway vehicle dynamics, if one neglects the aerodynamic forces and the traction/brake forces, the external force F consists of the wheel-rail contact forces. The problem at hand, therefore, is to evaluate the dynamic response of the system when a train circulates in the railway network, with a focus on the wheel-rail contact forces. Their variations are mainly due to the changing geometry of the track when the train runs at a given speed, as well as the time-history of the y vector. The equations of motion may be solved fully at each of a series of very small time steps to update the initial value problem for the succeeding time step, all the while ensuring that the size of each step is small enough for the solution to remain numerically stable [START_REF] Polach | Simulation of railway vehicle dynamics[END_REF]. The number of dofs of MBS can be reduced by introducing inflexible, kinematical connecting elements like joints and bearings, which enforce constraints on the movements of the bodies. The masses representing different components of the system are interlinked through the different springs, dampers and actuators, each of which itself is assumed to be massless. The specific design of rolling stocks leads to a distinction between the wheelset and the other bodies of the vehicle. The wheel may undergo high-frequency variations of wheel-rail forces, for instance, due to contact jumps between the tread and the flange or sudden changes of the rail geometry such as in turnouts. On the other hand, the other bodies reactions are low-pass filtered by the suspensions. An MBS modelling of the rolling stock is relevant in railway dynamics because most of its flexibility is ensured by the links between its bodies.

The multibody systems approach for railway applications

Several MBS software, specifically adapted for railway vehicle dynamics applications, have seen substantial developments in recent years. These include, but are not limited to, GENSYS, MEDYNA, NUCARS, SIMPACK, Universal Mechanism, VAMPIRE, VI-Rail, and VOCO. These simulation packages allow the automatic formulation of the equations of motion for many types of mechanical systems, via a user interface for the input of the vehicle and track parameters, which aids the modelling procedure. In the following sections, some key aspects of the railway multibody systems approach are highlighted, using the MBS code VOCO as an example when required. This is by no means a complete coverage of all the formalisms used in different MBS codes, but rather an overview of how railway systems are commonly modelled and analysed in multibody system dynamics, with a particular focus on the modelling of the wheel-rail contact. Some modelling and simulation techniques used in MBS software are also highlighted later on in Chapter 4 when handling specific case studies in railway dynamics with VOCO. For a more generalised in-depth analysis, the interested reader may refer to the books by Shabana et al. [START_REF] Shabana | Railroad vehicle dynamics: a computational approach[END_REF], Shabana [START_REF] Shabana | Dynamics of multibody systems[END_REF], or Popp and Schiehlen [START_REF] Popp | Ground vehicle dynamics[END_REF].

Coordinate systems and associated degrees of freedom

The equations of motion describing the behaviour of each component are of central importance in any MBS model. Equally important is the definition of a coordinate system that is used to formulate these equations. As opposed to a standard multibody approach where bodies may undergo large displacements and rotations [START_REF] Shabana | Dynamics of multibody systems[END_REF], guided transportation allows the following additional assumptions:

• The relative displacements and rotation with respect to the track geometry are limited as the vehicle motion is prescribed to follow the • Although the absolute displacements and rotations may be large, the relative ones between the bodies constituting the vehicle are generally small, except for the wheel rotation.

Choosing an appropriate coordinate system that is adapted for railroad vehicle formulations ensures the sparse matrix structure of the dynamic equations and the ease of implementation in MBS algorithms. In VOCO, the general curvilinear track coordinate system is implemented using the theory from [START_REF] Ayasse | Dynamique ferroviaire en coordonnées curvilignes: tracé de voie[END_REF]. This formulation is advantageous as it enables vehicle dynamics analysis on straight, curved as well as transition curves without making the kinematic and dynamic equations of motion too complex.

For the analysis, several coordinate systems are needed. The inertial coordinate system R 0 (X 0 , Y 0 , Z 0 ) is an absolute Cartesian reference where the origin is fixed at the beginning of the track and located in the track plane at the track centre-line, as shown in Figure 1.1. Following the righthand system, X 0 denotes the initial rolling direction which is tangent to the track centre-line, Y 0 is the lateral axis pointing towards the left rail with respect to the rolling direction, and Z 0 is pointing upwards. The mobile track coordinate system R f (X f , Y f , Z f ) is located at the track length of a given body at a given time. Its origin is at the track level on the track centre-line, with the longitudinal direction X f , the lateral direction Y f and the vertical direction Z f similarly following the right-hand system. The origin of the track system moves with the speed of the vehicle V along the track centre-line, following the curves and the track gradients. The mobile body coordinate system R v (X v , Y v , Z v ) is a local system belonging to each body, as shown in Figure 1.1, where the origin of the coordinate system is the same as in the previous case, and the axes are aligned with the principal axes of the body [START_REF] Hoffmann | Dynamics of european two-axle freight wagons[END_REF]. The orientation of the local mobile body system is the result of three successive rotations as shown in Figure 1.2 with the yaw angle β around Z 0 , the pitch angle θ around Y f , and the roll angle ψ around Figure 1.3: Origin of the mobile track coordinate system and the cant angle ψ X v . The value of the yaw angle may be higher than 2π radians, for instance, in a closed test track. However, the pitch and roll angles are kept small in most railway applications. If β is large, the value of the relative yaw α of a body with respect to the track is small. θ is related to the slope, which does not usually exceed three per cent in most of the railway tracks. The roll angle ψ, also called the cant angle, is prescribed to help negotiate curves via an elevation of the outer rails, as shown in Figure 1.3, with 2e being the nominal track gauge.

Railway vehicle dynamics and motion scenarios

The dynamics of railway vehicles can be categorised into vertical, longitudinal, and lateral dynamics. The study of vertical dynamics deals with the railway vehicle's response to track irregularities, transition zones (turnouts, bridges etc.) or rail corrugation. Longitudinal dynamics is related to the effects related to the acceleration and braking actions of a train set as it travels along the track, which may lead to complex interactions between the different vehicles [START_REF] Wu | International benchmarking of longitudinal train dynamics simulators: results[END_REF]. The lateral dynamics on the other hand is concerned with the vehicle behaviour in straight or curved tracks and is generally related to motion stability, safety and track loading, as required by the European Standard EN 14363 [START_REF]Railway applications-Testing and simulation for the acceptance of running characteristics of railway vehicles-Running behaviour and stationary tests[END_REF].

Hunting is one of the important characteristics of the motion of railroad vehicles, which is linked to the conicity of the wheel profiles as well as the rigidity of the wheelset assembly. The hunting phenomenon can be defined as the lateral oscillatory motion of the rigid wheelset with respect to its initial position on the track, as shown in Figure 1.4. In the absence of friction forces, the wheelset undergoes a periodic sinusoidal motion called

λ = 2π Re Γ , (1.2)
where R is the nominal rolling radius, and Γ is the equivalent conicity calculated according to the European Standard EN 15302 [START_REF]Railway applicationswheel-rail contact geometry parameters -definitions and methods for evaluation[END_REF] or the UIC 519 leaflet [START_REF]Method for determining the equivalent conicity[END_REF].

If frictional forces are taken into account, they tend to mitigate the hunting motion. These forces depend on the vehicle speed. The critical speed of the vehicle is defined as the speed at which this hunting motion becomes unstable, which may also lead to derailment. Derailment refers to the complete loss of contact between the wheel and the rail following a wheel climb due to high lateral forces. A commonly used criterion for derailment is the ratio between the resultant wheel-rail lateral Y and vertical Q forces acting on the wheel. The limit of the Y /Q ratio can be determined using Nadal's formula [START_REF] Shabana | Railroad vehicle dynamics: a computational approach[END_REF]:

Y Q = tan γ f -µ 1 + µ tan γ f , (1.3)
where γ f denotes the flange angle, as shown in Figure 1.5, and µ is the coefficient of friction.

In steady curving, the cant plays an important role to ensure that the centrifugal forces are balanced out by the lateral forces. If the vehicle travels with a speed that is above the speed where the two forces balance out, the vehicle is said to have a cant deficiency. The cant deficiency is defined as:

cd = 2e V 2 gR c -ψ , (1.4)
where V is the train speed, R c is the radius of the curve, and g is the gravitational acceleration. The cant deficiency is at most 160 mm in France. When cd is negative, the train is said to be in cant excess. Cant excess represents the amount by which the vertical elevation between the two parallel rail needs to be reduced so that the vehicle speed becomes equal to the balance speed.

Vehicle and track modelling

In the context of railway multibody systems, a standard vehicle model consists of the carbody, the two bogie frames, and the four wheelsets, as shown in Figure 1.6. Each pair of wheelset is connected to its corresponding bogie frame through the primary suspension system, while the secondary suspension accounts for the connection between the bogies and the carbody.

For the modelling of low-frequency motion, the bodies are generally assumed to be rigid, with six dofs assigned to the centre of mass. This is considered to be sufficiently adequate to investigate problems such as ride comfort, vehicle stability, and curving behaviour. VOCO normally considers each body of the vehicle as a rigid structure, although the flexibility of suspended bodies (carbody, bogie frame, etc.) can also be incorporated through FE results for studies involving ride comfort, stress analysis of lightweight vehicles, etc. As far as the wheelset is concerned, only its torsional flexibility is taken into account. For the vehicle shown in Figure 1.6, each component has six degrees of freedom1 , resulting in a total of 42 dofs for the entire model. Wheel axle structural flexibility may have significant effects on the vehicle-track interaction, especially in applications where the interest lies in evaluating the damage, wear, and fatigue effects, or in the vibrations of the wheelsets and rails [START_REF] Bruni | State-of-the-art and challenges of railway and road vehicle dynamics with multibody dynamics approaches[END_REF]. The tasks of modelling the primary and secondary suspension systems, as well as the different connecting elements, is an important sub-aspect of the vehicle modelling procedure. A traditional ballasted track consists of rails, sleepers, rail pads and fasteners connecting the rails and the sleepers, the ballast and sub-ballast layers and the foundation. In order to model the different components shown in Figure 1.7, the rail is modelled using beams that are discretely supported by a bed of springs and dampers. For low-frequency studies, a co-running track model may be considered, which is a lumped parameter representation of different components and linked to each wheelset, moving with the same speed as the vehicle. In VOCO, it is possible to take into account the dynamic flexibility of the track via a model where the rail is modelled using finite elements beams, the sleepers by rigid solids, and other components (ballast, subgrade etc.) through linear or nonlinear support elements. A different possibility can be to adopt a co-simulation approach together with a FE software modelling the track.

Wheel-rail contact modelling in the context of multibody systems

The modelling of the dynamic interaction between the wheel and the rail is a complex task, and one of the fundamental problems that must be addressed while developing any MBS formulation for railroad vehicle dynamics. The wheel-rail contact problem is characterised by large forces that are transmitted in small zones, which may change abruptly depending on the contact geometry, the kinematics, and the frictional properties. In general, the constraint approach and the elastic approach are the two main approaches that can be used in MBS formulations to implement the wheelrail contact module [START_REF] Shabana | Railroad vehicle dynamics: a computational approach[END_REF].

The constraint approach describes the wheel-rail contact using nonlinear kinematic constraint equations, and the normal contact forces are calculated as constraint forces. Virtual interference or local deformation between the wheel and the rail is not permitted, and the contact is considered to be rigid. This leads to a model in which the wheelset has four degrees of freedom, with the vertical and roll motions being dependent on the lateral motion of the wheelset. The constraint approach also only allows one contact point per wheel-rail pair, and a non-conformal contact is assumed. In the formulations based on the elastic approach, no additional constraint is introduced. The normal contact force is expressed as a compliant force with assumed stiffness and damping between the two bodies influenced by Hertz' contact theory. The wheelset has six degrees of freedom, and separation between the two is allowed. This in turn permits more precise modelling of many simulation scenarios, such as derailment, which is not always possible with the constraint approach. Regardless of the approach that is used, the wheel-rail contact module must broadly handle three main functions [START_REF] Bruni | State-of-the-art and challenges of railway and road vehicle dynamics with multibody dynamics approaches[END_REF]:

• Determining the contact points on the wheel and rail surface profile through a contact point search algorithm;

• The contact kinematics, which involves the calculation of the creepages at the contact points;

• The contact mechanics, which is the problem of determining the normal and tangential components of the contact forces.

The modelling of the impact force also constitutes an important aspect, especially in simulations containing switches and crossings, due to the impact-contact seen at the crossing nose. The most basic contact force models are those based on Hertz' theory [START_REF] Hertz | Über die Berührung fester elastischer Körper[END_REF], where no energy dissipation is considered during the impact process. However, even a collision between two perfectly elastic bodies dissipates energy due to the internal damping, and some initial kinetic energy is lost through vibrations, heat, sound and other forms [START_REF] Goldsmith | Impact: The theory and physical behavior of colliding solids[END_REF]. Thus, elastic models must be further enhanced to incorporate this energy loss. One of the most widely used concepts to take energy dissipation and the nonlinear viscous-elastic behaviour into account is through a coefficient of restitution, and several models have been developed using this approach in recent years [START_REF] Machado | Compliant contact force models in multibody dynamics: Evolution of the Hertz contact theory[END_REF][START_REF] Skrinjar | A review of continuous contact-force models in multibody dynamics[END_REF].

Contact search

The contact may be treated as a two or three-dimensional problem, depending on the parameters that are used to describe the wheel and rail profiles. Determining the location of the contact points online is not straightforward, and most elastic force models reduce the three-dimensional contact problem to a two-dimensional one. In the two-dimensional approach, the contact points are searched solely in the vertical plane. The location of the points of contact can be determined by using lookup tables or by using a discrete nodal search or else by solving a set of algebraic equations [START_REF] Magalhães | Implementation of a non-Hertzian contact model for railway dynamic application[END_REF]:

• Lookup tables may be used when the computational efficiency is of primary importance, and the interpolation of a pre-computed table replaces the contact search, given the relative kinematics of the wheel and rail;

• In the nodal search approach, the profile of the wheel and rail is described using discrete nodal points, and the distance between these nodes is used to determine the points on the wheel and the rail that may come into contact;

• In the algebraic approach, the geometric constraints are described using normal and tangential vectors at the wheel and rail that determine the potential points of contact by imposing a minimum distance condition.

A nodal search with two curves representing the transversal profiles of the wheels and rails, or lookup tables are usually the strategies considered for MBS formulations when solving algebraic equations online may not be very efficient.

Contact kinematics

Following the contact detection, the longitudinal, lateral and spin creepages at the contact point are determined through kinematic analysis of the two contacting surfaces. First introduced by Carter in 1926 [START_REF] Carter | On the action of a locomotive driving wheel[END_REF], "creepage" refers to the small slip in the contact area due to the difference between the tangential strains of two elastic bodies in the presence of friction. If there is a relative translation motion between the two bodies, the rolling motion is accompanied by sliding, which splits the contact area into two zones. In the adhesion or stick region, the relative velocity of the particles is compensated by the elastic deformation of the bodies and the two bodies "stick" to each other. In the remaining region, also within the contact area but outside the adhesion zone, the relative velocity of each particle exceeds the elastic contribution and the particles start to slide over each other, a phenomenon referred to as micro-slip [18, p. 242]. This micro-slip effect was also observed by Reynolds, who was studying the case of a rubber cylinder rolling on a metal plane [START_REF] Reynolds | On rolling-friction[END_REF].

For railway applications, creepages are normalised relative velocities between both contacting bodies considered to be rigid. The expressions for the creepages described in the tangent plane used in the MBS context are given as:

ν x L,R = V -r L,R ω ∓ e α V , (1.5)
ν y L,R = ẏ V -α cos γ , (1.6)
φ L,R = - sin γ R , (1.7)
where ν x1 and ν y are called longitudinal and lateral creepages, and φ is called the spin. The subscripts L and R denote the left and right wheel, respectively. R is the nominal rolling radius of the wheel, r is the local rolling radius of the wheel, V is the translational velocity, ω is the angular velocity, α is the yaw angle of the wheelset, and γ is the contact angle. The description of the different creepages, commonly considered for railway applications, is shown in Figure 1.8. It should be noted that some terms of the relative velocities have been neglected in the expressions (1.5)-(1.7): some authors prefer to compute the exact relative velocities to evaluate the creepages [START_REF] Vollebregt | Detailed wheel/rail geometry processing using the planar contact approach[END_REF].

Contact mechanics

Contact mechanics is defined as the study of unilateral interaction of two or more bodies over a point or a surface. As all real bodies are essentially soft or flexible, the surface deformation implies that the resulting forces are transmitted between them through a "contact area", which is now under the influence of a stress distribution. This contact area can, in turn, be concentrated in a small region when compared to the dimensions of the contacting bodies, e.g., the contact between a wheel and a rail, or spread over a larger or diversified area, e.g., the contact between a cylinder and an elastic plane. Considering the contact to be locally concentrated permits the approximation of the contacting bodies as elastic half-spaces bounded by a plane. Solving the problem of normal contact between two elastic bodies consists of determining the zone of contact, i.e., the contact area, as well as the associated normal stress distribution in the contact patch. The solution of the tangential problem, on the other hand, involves using the normal contact results and contact kinematics to evaluate the tangential components of the wheel-rail contact forces, as well as the division of the slip and adhesion regions within the contact patch.

The first elastic theory of contact mechanics is credited to the seminal analytical work of Hertz [START_REF] Hertz | Über die Berührung fester elastischer Körper[END_REF], although with strong assumptions concerning the curvatures of the contacting surfaces. These assumptions may lead to approximate results in railway applications where the curvatures of the profiles change along the width of the rail section. The Hertzian solution is characterised by an elliptical contact patch, as shown in Figure 1.9, with a semi-ellipsoidal normal pressure distribution acting within the contact area. The solution of nonlinear contact problems using detailed finite element (FE) or boundary element (BE) methods has grown together with the innovations in modern computers. In decreasing order of complexity, it is possible to classify the various approaches of computational contact mechanics for continuous media as follows:

• FE methods, for cases with concentrated as well as diversified contact problems;

• BE methods, with or without the half-space assumption for locally concentrated contact problems;

• Analytical or semi-analytical approaches based on Hertz' theory.

The use of discrete element methods offers a more realistic approximation of the contact conditions in real-life applications. Discrete methods are based on the principle of expressing the elastic field in terms of the loading in the potential contact, which then leads to an influence function method [21, p. 49]. Variational and BE methods allow approximating the contacting bodies as half-spaces, for which the influence functions may be determined analytically. The elastic response of the contacting bodies to surface loads deviates from the response of the half-space, and there exists the possibility of introducing the true behaviour of the bodies using numerically computed influence coefficients via FE methods. Several authors have assessed the contact problem in this manner, including Ahmadi, Keer and Mura [START_REF] Ahmadi | Non-Hertzian contact stress analysis for an elastic half space-normal and sliding contact[END_REF] who presented a quadratic programming solution for non-Hertzian normal and sliding elastic contact. This method was subsequently further developed and rigorously proved by Kalker [START_REF] Kalker | Contact mechanical algorithms[END_REF] and is now implemented in the commercial program CONTACT [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF]. One can additionally cite the matrix inversion method (MIM) by Johnson [START_REF] Johnson | Contact mechanics[END_REF][START_REF] Toumi | Modélisation numérique du contact roue-rail pour l'étude des paramètres influençant les coefficients de Kalker : Application à la dynamique ferroviaire[END_REF], the methods by Kalker-van Randen [START_REF] Kalker | A minimum principle for frictionless elastic contact with application to non-Hertzian half-space contact problems[END_REF], de Mul-Kalker-Fredriksson [START_REF] De Mul | The contact between arbitrarily curved bodies of finite dimensions[END_REF], and Paul and Hashemi [START_REF] Paul | Contact pressures on closely conforming elastic bodies[END_REF] where various types of elements are chosen to discretise the contact area. The ad hoc methods presented by Reusner [START_REF] Reusner | Druckflächenbelastung und oberflächen-verschiebung im wälzkontakt von rotationskörpern[END_REF] and Nayak-Johnson [START_REF] Nayak | Pressure between elastic bodies having a slender area of contact and arbitrary profiles[END_REF] for roller bearings, and by Knothe & Le The for wheel-rail contact [START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with nonelliptical contact area[END_REF], also enable a good estimation of the contact properties. The full normal problem has also received much attention for the contact of rough surfaces [START_REF] Polonsky | A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques[END_REF][START_REF] Vollebregt | A new solver for the elastic normal contact problem using conjugate gradients, deflation, and an FFT-based preconditioner[END_REF].

FE methods can be used for a wide range of applications and are not restricted to the half-space assumption generally used with BE methods. Thus, they allow the modelling of the real 3D geometry and can be used to handle both small and large deformation contact problems. The inequality constraint on the deformation field means that the contact problem is highly nonlinear even when we consider the case of linear elasticity with small strain. Frictional contact presents more difficulties, as not only are the inequality constraints in the normal direction present, but there is also a special constitutive behaviour in the tangent direction at the contact interface [34, p. 112]. Most FE codes nowadays handle contact problems using either the penalty or the Lagrange multiplier method. Each of the methods has its advantages and disadvantages, the Lagrangian multiplier method enforces the contact constraints precisely, but the additional unknowns may lead to issues with the numerical stability. The penalty method on the other hand is easier to implement from a numerical viewpoint, but the contact constraints are satisfied exactly only in the limit of infinite penalty values. In comparison, the augmented Lagrangian method offers more stability and robustness without any additional unknowns but is slow to converge. Certain recent works have also used Nitsche's method for frictional contact problems, which offers good convergence and robustness as compared to the classical penalty methods or mixed methods [START_REF] Mlika | An unbiased Nitsche's formulation of large deformation frictional contact and self-contact[END_REF][START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF]. The book by Kikuchi and Oden [START_REF] Kikuchi | Contact problems in elasticity: A study of variational inequalities and finite element methods[END_REF], the monograph by Wriggers [START_REF] Wriggers | Computational contact mechanics[END_REF] or the book by Laursen [START_REF] Laursen | Computational contact and impact mechanics[END_REF], as well as the references mentioned within them offer an in-depth study for readers who wish to understand the theoretical and numerical background of modern finite element techniques for contact mechanics.

Normal contact: Several researchers have treated the normal wheel-rail contact problem using FE methods. Yan and Fischer [START_REF] Yan | Applicability of the Hertz contact theory to rail-wheel contact problems[END_REF] have investigated the applicability of the Hertz contact theory to wheel-rail contact problems by applying three-dimensional finite element models. The contact pressure distributions for different types of rails concerning different initial contact positions were found to agree with the Hertz theory only if the surface curvature did not vary, or if plastification did not occur. Telliskivi and Olofsson [START_REF] Telliskivi | Contact mechanics analysis of measured wheel-rail profiles using the finite element method[END_REF] have studied the influence of the half-space and elasticity assumptions in the Hertzian as well as Kalker's complete theory when compared to elastoplastic FE solution, as shown in Figure 1.10. The results using Hertz and CONTACT methods provided a maximum contact pressure that was larger than the FE solution. A similar study has been carried out by Wiest et al. [START_REF] Wiest | Assessment of methods for calculating contact pressure in wheel-rail/switch contact[END_REF] for the contact pressures obtained in the case of railway switches. The results, in this case, were found to be in good agreement with the Hertz theory and CONTACT within the elastic limit, whereas plastification was shown to decrease the maximum contact pressure considerably.

The effect of residual stress distributions, due to plastic deformations From [START_REF] Telliskivi | Contact mechanics analysis of measured wheel-rail profiles using the finite element method[END_REF] in the subsurface area, on the fatigue limit and cracking tendencies of the material has been investigated by Bijak-Żochowski and Marek [START_REF] Bijak-Żochowski | Residual stress in some elasto-plastic problems of rolling contact with friction[END_REF]. The influence of wheel and rail profiles on the distribution of contact zones and stresses has been studied by Sladkowski and Sitarz [START_REF] Sladkowski | Analysis of wheel-rail interaction using FE software[END_REF] and used for the design improvement of operating wheels. Zhao and Li [START_REF] Zhao | The solution of frictional wheel-rail rolling contact with a 3D transient finite element model: Validation and error analysis[END_REF] have used a threedimensional transient FE model to solve the rolling contact problem for arbitrary geometry. This model has been extended to include the influence of plastification on the contact patch in [START_REF] Zhao | A three-dimensional finite element solution of frictional wheel-rail rolling contact in elasto-plasticity[END_REF]. The coupling between the normal and tangential contact problems was found to become stronger when plastic flow occurs. Blanco-Lorenzo et al. [START_REF] Blanco-Lorenzo | A contact mechanics study of 3D frictional conformal contact[END_REF] have presented the study of three-dimensional contact with friction when taking into account the effects of conformity. Toumi et al. [START_REF] Toumi | Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods[END_REF] have developed a three-dimensional rolling wheel-rail contact model to study the normal and the tangential contact problems as well as the creep force characteristics in elasticity and elastoplasticity using explicit and implicit integration schemes. Approaches focused on the local tribological analysis of the wheel-rail contact and the rail corrugation process have also been presented in recent years by Saulot [START_REF] Saulot | Analyse tribologique du contact roue-rail : modélisation et expérimentations : cas de l'usure ondulatoire[END_REF] and Duan [START_REF] Duan | Numerical tribology of the wheel-rail contact : Application to corrugation defect[END_REF].

Even though the use of FE and BE methods for rolling wheel-rail contact problems has been around for around half a century now, from the point of view of practical computational times, discrete methods continue to struggle in terms of their capability to be used online in dynamic MBS simulations.

A very fine discretisation is usually required with FE approaches to satisfactorily represent the contact patch boundaries, while classical BE methods normally require iterative techniques over large potential contact grids to verify the stress constraints. In MBS codes, a compromise is sought between the level of accuracy that is deemed acceptable, and the time it takes to solve the contact problem. For the normal contact problem, most MBS software use either Hertz' theory, or methods based on the theory of virtual penetration (VP) of the contacting surfaces [START_REF] Piotrowski | Wheel-rail contact models for vehicle system dynamics including multi-point contact[END_REF]. The contact area in VP approaches is approximated using the area in which the surfaces would interpenetrate geometrically (overlap) if there was no deformation. This overlapping area is bigger than the actual contact patch, and a scaling factor is used to prescribe the virtual penetration such that the interpenetration zone approaches the real contact area. The technique used for determining the scaling parameter is what broadly differentiates these VP approaches [START_REF] Sichani | Comparison of non-elliptic contact models: Towards fast and accurate modelling of wheel-rail contact[END_REF]. There exists in parallel another family of models wherein the nonelliptical part of the contact area is approximated using a series of individual ellipses. This so-called multi-Hertzian method, proposed by Pascal and Sauvage [START_REF] Pascal | The available methods to calculate the wheel/rail forces in non Hertzian contact patches and rail damaging[END_REF], is further developed by Ayasse et al. [START_REF] Ayasse | Paramètres caractéristiques du contact roue-rail[END_REF] to an analytical approach. Some questions raised regarding these models, when compared against VP methods in [START_REF] Piotrowski | Wheel-rail contact models for vehicle system dynamics including multi-point contact[END_REF], have been addressed recently by Pascal [START_REF] Pascal | Solving conformal contacts using multi-Hertzian techniques[END_REF].

Tangential contact: One of the first theories to evaluate the creep forces was presented by Carter [START_REF] Carter | On the action of a locomotive driving wheel[END_REF]. Nowadays, the tangential creep forces in MBS codes are commonly computed with the FASTSIM algorithm [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF], assuming a steady state. The simplified theory behind FASTSIM is based on the computation of the creepage coefficients c ij with the help of Kalker's exact theory [START_REF] Kalker | On the rolling contact of two elastic bodies in the presence of dry friction[END_REF]. These coefficients are derived from the initial slope of the tangent forces applied on an elliptical contact patch for small creepage values. A newer version of the FASTSIM algorithm offers second order accuracy by taking precise care of the contact conditions in the numerical integration scheme employed [START_REF] Vollebregt | FASTSIM2: A second-order accurate frictional rolling contact algorithm[END_REF]. The FaStrip algorithm [START_REF] Sichani | An alternative to FASTSIM for tangential solution of the wheel-rail contact[END_REF] is proposed as an alternative to FASTSIM, and is based on the strip theory, which extends the two-dimensional rolling contact solution to three-dimensional contacts [START_REF] Kalker | A strip theory for rolling with slip and spin[END_REF]. Corrections are introduced in the original strip theory to obtain accurate estimations for any contact ellipse size, and it is combined with a FASTSIM type algorithm to handle spin. Both FASTSIM and FaStrip are iterative algorithms that require a discretisation of the contact ellipse into rectangular elements, and express the tangent forces by integrating the tangential shears over the discretised contact area. As opposed to this, fast analytical methods such as the Vermeulen-Johnson model [START_REF] Vermeulen | Contact of non-spherical elastic bodies transmitting tangential forces[END_REF], the Shen-Hedrick-Elkins model [START_REF] Shen | A comparison of alternative creep force models for rail vehicle dynamic analysis[END_REF] and Polach's model [START_REF] Polach | A fast wheel-rail forces calculation computer code[END_REF] use formulas giving the tangent forces directly as a function of the unsaturated shears. Approaches such as the Book of Tables by Kalker (USETAB) [START_REF] Kalker | Book of tables for the Hertzian creep-force law[END_REF] use detailed methods like Kalker's complete theory [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF] to develop lookup tables that are used in vehicle dynamics simulations.

The coefficient of friction in vehicle dynamics is generally considered as a known parameter and is a constant for the entire simulation. This assumption is clearly not always accurate as the friction coefficient depends on the operating conditions, and varies a lot between different circumstances like the weather, or the presence of sand, foliage, lubricants or contaminants, water, ice and snow etc. A detailed investigation on the various tribological aspects of creep force modelling, arising from third body layers, can be found in the state-of-the-art paper by Vollebregt et al. [START_REF] Vollebregt | Challenges and progress in the understanding and modelling of the wheel-rail creep forces[END_REF]. Numerical algorithms based on Kalker's simplified theory, which enable the modelling of variable friction, have been presented by different authors [START_REF] Giménez | Introduction of a friction coefficient dependent on the slip in the Fastsim algorithm[END_REF][START_REF] Piotrowski | Kalker's algorithm Fastsim solves tangential contact problems with slip-dependent friction and friction anisotropy[END_REF][START_REF] Rovira | Application of Fastsim with variable coefficient of friction using twin disc experimental measurements[END_REF][START_REF] Vollebregt | Fastsim with falling friction and friction memory[END_REF]. The Modified FASTSIM algorithm presented in [START_REF] Spiryagin | Creep force modelling for rail traction vehicles based on the Fastsim algorithm[END_REF] incorporates a variable contact flexibility which depends on the ratio of the slip area to the area of adhesion, together with a slip-velocity-dependent friction coefficient. A stepping method based on FASTSIM to solve non-steady rolling contact problems has been presented in [START_REF] Eddhahak-Ouni | Modélisation de l'usure superficielle dans les contacts roulants : mise en place des méthodes adaptées au cas non stationnaire des contacts came galet[END_REF]. A falling friction variant is also available with Polach's creep force model [START_REF] Polach | Creep forces in simulations of traction vehicles running on adhesion limit[END_REF].

The main drawback for all of the above cited fast tangential methods is that they remain restricted to elliptical contact patches, although FaStrip has been extended to non-Hertzian patches [START_REF] Sichani | A fast wheel-rail contact model for application to damage analysis in vehicle dynamics simulation[END_REF]. Piotrowsky et al. [START_REF] Piotrowski | The Kalker book of tables for non-Hertzian contact of wheel and rail[END_REF] have extended USETAB to be used in non-Hertzian cases, approximating the contact patch by a single double-elliptical contact. Broadly, two approaches may be used to extend FASTSIM type algorithms to non-elliptical contact cases [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF]. The first approach consists of regularising the non-Hertzian contact patch to a single equivalent ellipse and using the global creepages to determine the creep forces using one of the various Hertzian approaches [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF][START_REF] Vermeulen | Contact of non-spherical elastic bodies transmitting tangential forces[END_REF][START_REF] Shen | A comparison of alternative creep force models for rail vehicle dynamic analysis[END_REF][START_REF] Polach | A fast wheel-rail forces calculation computer code[END_REF]. While this approach works well in dynamic studies, it cannot be used to study the contact stresses. The second approach consists in extrapolating the original algorithm to a non-elliptic patch and using creepages defined locally in cells or strips [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF][START_REF] Knothe | Determination of the tangential stresses and the wear for the wheel-rail rolling contact problem[END_REF][START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF].

For the wheel-rail contact problems addressed in this thesis, the results of the contact detection and contact kinematics processes are supposed to be known beforehand, assessed through the MBS code VOCO or taken from the literature. To investigate the normal and tangential contact mechanics, Kalker's variational theory for elastic half-spaces [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF] is used as the reference. This approach, also referred to as Kalker's complete theory for the rolling contact problem, is outlined in the following section.

The rolling contact problem with Kalker's complete theory

The boundary element method is an appealing and computationally viable technique that allows solving various problems, including the contact between two bodies. It is especially relevant when the interest of the user lies in the evaluation of only surface stresses and surface displacements, rather than determining these quantities at every point in the interior of the bodies. BE methods are particularly advantageous when the domain of study is large, as only the surface boundary needs to be discretised now. If we consider for example the problem of contact between two elastic bodies, the mesh only applies to the potential contact surface unlike in the finite element method which requires a mesh of the 3D geometry of both bodies. Different formulations, including direct, indirect, and variational approaches, of BE methods are available in the literature [START_REF] Johnson | Contact mechanics[END_REF][START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF][START_REF] Andersson | The boundary element method applied to twodimensional contact problems with friction[END_REF][START_REF] Björklund | A numerical method for real elastic contacts subjected to normal and tangential loading[END_REF]. The following section presents Kalker's variational contact theory, which is implemented in the commercially available software CONTACT [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF]. CONTACT has been developed extensively over the past years as compared to Kalker's original version, in particular by Dr E. A. H Vollebregt, broadening its range of applicability by including velocity-dependent friction laws [START_REF] Vollebregt | Quasi-static analysis of twodimensional rolling contact with slip-velocity dependent friction[END_REF], the effects of roughness and contamination (the so-called "third body layer") [START_REF] Vollebregt | Numerical modeling of measured railway creep versus creep-force curves with contact[END_REF], extensions for solving conformal contact problems [START_REF] Vollebregt | Solving conformal wheel-rail rolling contact problems[END_REF], and via new enhanced numerical solvers [START_REF] Vollebregt | A new solver for the elastic normal contact problem using conjugate gradients, deflation, and an FFT-based preconditioner[END_REF][START_REF] Zhao | A fast nonlinear conjugate gradient based method for 3D concentrated frictional contact problems[END_REF]. Some added enhancements are out of the scope of the present work, and we are mainly interested in the contact problem of three-dimensional homogeneous elastic bodies, where the contact occurs in a concentrated zone under the influence of dry friction.

In 1979, Kalker published his exact three-dimensional rolling contact theory, which applies the principles of virtual work and complementary virtual work to solve the contact problem. This was based on the previous works by Fichera [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: il problema di signorini con ambigue condizioni al contorno[END_REF] or Duvaut and Lions [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] who established the variational theory of frictional contact. The exact theory has been presented in detail by Kalker in his monograph [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF], and the following assumptions are made initially:

• The bodies are elastic, homogeneous and isotropic;

• The contacting surfaces have a continuous profile;

• The contact area is essentially flat and concentrated with respect to the overall dimensions of the contacting bodies;

• The inertial effects are small compared to the contact stresses and may be neglected which implies a quasi-static state.

The assumptions made here permit using the half-space approach, where the two bodies are assumed to be semi-infinite elastic half-spaces. The analytical solutions for the response of the bodies to surface loadings have already been presented in detail by Boussinesq and Cerruti [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF][START_REF] Cerruti | Richerche intorno all equilibrio de corpi elastici isotropi[END_REF] and are summarised in Appendix A.

The elastic half-space and quasi-identity

We consider the case of two bodies of arbitrary surfaces that come into contact. If the contact is non-conformal, the contact stresses are highly concentrated in a small region of interest which lies close to the point of contact. In many practical applications such as the wheel-rail contact, the dimensions of this region of interest, i.e., the contact area, are small as compared to the overall dimensions of the two bodies. One may thus assume that these localised stresses are largely independent of the overall shape of the bodies, and rather depend on the geometry of the contact interface. The contact stresses can be consequently evaluated by considering each elastic body as a semi-infinite space bounded by a plane surface. This assumption of each body as an elastic half-space is highly useful, as it opens up the possibility of using a wide variety of theories developed for the loading of elastic half-spaces. The half-space approximation can be visualised in Figure 1.11. The half-space consideration for frictional contact problems is generally accompanied by the quasi-identity assumption. For two homogeneous and isotropic bodies with modulus of rigidity and Poisson's ratio (G 1 , ν 1 ) and (G 2 , ν 2 ), respectively, the combined modulus of rigidity G is given as:

1 G = 1 2 1 G 1 + 1 G 2 , (1.8)
the combined Poisson's ratio ν as:

ν G = 1 2 ν 1 G 1 + ν 2 G 2 , (1.9)
Figure 1.11: The half-space approximation. From [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF] and the difference parameter K as:

K G = 1 4 1 -2ν 1 G 1 - 1 -2ν 2 G 2 . (1.10)
Two bodies are said to be quasi-identical if the difference parameter K = 0. This is true in the case of contact between a steel wheel and steel rail, and a small value of K is often neglected. Note that two geometrically symmetrical bodies are also quasi-identical, irrespective of whether they are considered as half-spaces or not. The addition of the quasi-identity assumption enables the separation of the normal and tangent contact problems in the half-space theory, which can then be solved sequentially. The existence and uniqueness of quasi-identical frictional contact problems can also be proved, as done in the works of Fichera [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: il problema di signorini con ambigue condizioni al contorno[END_REF] and Duvaut and Lions [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF], as opposed to existence only for non quasi-identical frictional contact problems [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF].

The general contact problem

According to Hooke's law, the linear elasticity relations between the stresses σ ij and the strains e hk for the interiors of the two contacting bodies may be written as [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF]:

σ ij = E ijhk e hk , (1.11)
where E ijhk are the elastic constants.

Considering the half-space assumption, the contact problem can be brought into surface-mechanical form, where no attention is paid to the stresses, strains or displacements in the interior of the contacting bodies. The relationship between the surface Figure 1.12: The general contact problem stresses p and the surface displacements u in a steady state can be written as:

u(x) = C A(x, x ′ )p(x ′ )dC , (1.12)
where p = [p n , p x , p y ], and A(x, x ′ ) is the influence function matrix describing the displacement difference at x arising due to a unit load acting on x ′ . The influence coefficients depend on the material properties, as well as the geometries of the contacting bodies, and are obtained through the theory of Boussinesq and Cerruti.

Several quantities need to be defined to express the general contact problem. Consider two bodies of arbitrary shape defined by z 1 (x, y) and z 2 (x, y), respectively, that come into contact as shown in Figure 1.12. The origin of the Cartesian coordinate system is at the geometrical point of contact. When the two contacting bodies are subjected to compression, the surfaces deform locally, and the deformed distance e z (x, y) can be given as:

e z (x, y) = h(x, y) + u z (x, y) -δ , (1.13)
where h(x, y) = z 1 (x, y) -z 2 (x, y) is the separation or undeformed distance between the two bodies, u z (x, y) = u z1 (x, y) + u z2 (x, y) is the normal elastic displacement, and δ is the rigid body approach at the point of geometric contact. For the tangential contact problem, the important quantity to be defined is the relative slip velocity s t of two particles on the surfaces of the two contacting bodies:

s t = w t + ut /V , (1.14)
where the relative rigid slip w t describes the velocity by which the undeformed surfaces move with respect to each other, u is the time derivative of the surface displacements and V is the rolling velocity. For concentrated contact problems, w t is given by:

w t = [ν x -yφ ν y + xφ] . (1.15)
The complete contact problem consists of determining the contact area C, its subdivision into adhesion and slip areas A and S, respectively, and the surface tractions p such that the following contact conditions are satisfied for the considered domain [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF]:

normal problem: ∀x ∈ C : e z = 0, p n ≥ 0 , (1.16) ∀x / ∈ C : e z > 0, p n = 0 , (1.17)
tangential problem:

∀x ∈ A : ||s t || = 0, ||p t || ≤ g , (1.18)
∀x ∈ S :

||s t || > 0, p t = -g s t ||s t || , (1.19) ∀x / ∈ C : p t = 0 , (1.20)
p t = [p x , p y ] , (1.21)
g(x) = µp n (x) . (1.22)
From the above condition, the contacting bodies cannot interpenetrate in the deformed state, and the normal pressure within the contact area is assumed to be always compressive. The surface tractions also vanish at the boundary of the contact area. The tangential tractions are limited by Coulomb's friction law, and no slip occurs where the tangential traction falls below the traction bound. Finally, the tangential tractions in the slip area are opposite to the slip direction. The quasi-static frictional contact problem can now be considered as a variational inequality and written as a minimisation problem [21, p. 175]:

min p ϕ = P 1 2 p(x)A(x,x')p(x ′ ) + h(x)p n (x) + [W t (x) -u t ′ (x)]p t (x) dP (1.23) sub ∀x ∈ P : p n (x) ≥ 0, ||p t (x)|| ≤ g(x) , W t = w t V (t -t ′ ) .
Here, u t ′ (x) denotes the tangential displacements [u x , u y ] at the previous time step and P denotes the potential contact area. With the quasiidentity assumption and the accompanying decoupling of the normal and tangential problems, the contact problem can be solved sequentially. Solving the contact problem between two elastic solids sequentially first consists in determining the contact area, as well as the normal stress distribution acting on it. The normal contact results are subsequently used to determine the tangential surface tractions. This one-step outer procedure is called "Johnson's process". When the two bodies have different elastic properties, the tangential tractions affect normal displacement differences and vice versa, and the two problems are strongly coupled through the difference parameter K (see Appendix A). Solving this coupled problem consists of evaluating normal and tangential problems alternatively until the iteration process converges to a tolerance value with the so-called "Panagiotopoulos" process [START_REF] Kalker | Rolling contact phenomena: Linear elasticity[END_REF].

As the contact area is not known in advance, the minimisation problem is solved for an arbitrary potential contact area P, wherein the contact conditions are evaluated. Several discretisation strategies for the potential contact area have been used by different authors: the triangular elements used by Kalker and van Randen [START_REF] Kalker | A minimum principle for frictionless elastic contact with application to non-Hertzian half-space contact problems[END_REF], equal rectangular elements by de Mul-Kalker-Fredriksson [START_REF] De Mul | The contact between arbitrarily curved bodies of finite dimensions[END_REF], and adapted rectangular mesh by Paul and Hashemi [START_REF] Paul | Contact pressures on closely conforming elastic bodies[END_REF]. In CONTACT, the potential contact area is discretised using rectangular elements. Kalker developed the original active set algorithms NORM, TANG, and KOMBI that were used to solve the frictional contact problem [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF]. Using rectangular elements results in influence coefficient matrices that are of the block-Toeplitz-Toeplitz-block (BTTB) form. In the recent versions of CONTACT, a purpose-built fast numerical solver for BTTB matrices called NormCG is used for the normal contact problem, which is based on a bound-constrained conjugate gradient method with a fast Fourier transform (FFT) pre-conditioner [START_REF] Vollebregt | A new solver for the elastic normal contact problem using conjugate gradients, deflation, and an FFT-based preconditioner[END_REF]. This approach replaces the active set NORM algorithm [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF]. Similarly, the nonlinear tangential problem with prescribed creepages is solved using case-specific iterative solvers relying on the active set strategy, all of which enforce the traction bounds in all the elements [START_REF] Zhao | A fast nonlinear conjugate gradient based method for 3D concentrated frictional contact problems[END_REF][START_REF] Vollebregt | Improving speed and accuracy of the frictional rolling contact model "contact[END_REF][START_REF] Vollebregt | A Gauss-Seidel type solver for special convex programs, with application to frictional contact mechanics[END_REF].

Normal contact problem using fast approaches

The case of normal contact between two continuous and non-conforming bodies was first treated analytically by Hertz in 1882 [START_REF] Hertz | Über die Berührung fester elastischer Körper[END_REF], by assuming friction-less contact and each body as an elastic half-space. The Hertzian solution is characterised by an elliptical contact patch, with a semi-ellipsoidal normal pressure distribution acting within the contact area. Analytical solutions of the normal contact problem remain restricted to few cases with simple geometries of the contact bodies. One may refer to the books by Johnson [START_REF] Johnson | Contact mechanics[END_REF] or Popov [START_REF] Popov | Contact mechanics and friction: Physical principles and applications[END_REF] for an overview of the different approaches present in literature. For non-elliptic or non-Hertzian contact conditions, several "virtual penetration" or VP methods have been developed in recent years to avoid the computation of the surface deformation, including the Linder model [START_REF] Linder | Verschleiss von eisenbahnrädern mit unrundheiten[END_REF], the Kik-Piotrowski model [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF], the semi-Hertzian model (STRIPES) [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF], the extended Kik-Piotrowski (EKP) model [START_REF] Liu | A non-Hertzian method for solving wheel-rail normal contact problem taking into account the effect of yaw[END_REF], the modified Kik-Piotrowski (MKP) model [START_REF] Sun | A robust non-Hertzian contact method for wheel-rail normal contact analysis[END_REF], and the modified semi-Hertzian method [START_REF] Zhu | Modification of the semi-Hertz wheel-rail contact method based on recalculating the virtual penetration value[END_REF] amongst others. The method ANALYN [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF] adds an approximation of the surface deformation to the virtual penetration methods. In this section, we look at the Hertzian approach in detail, along with some commonly used virtual penetration approaches.

Hertz' theory

For someone so closely associated with the field of contact mechanics, Heinrich Hertz' major research developments came in electromagnetism theory, achievements which have of course been immortalised by the naming of the standard unit of frequency in his honour. It is also within the context of this research that Hertz was first attracted to the contact of elastic solids. While working as a research assistant for Hermann von Helmholtz at the Physical Society of Berlin between 1880-1883, discussions on the optical phenomena of Newton's rings piqued Hertz' interest. While a great deal was already known about the interference patterns between the contacting glass lenses, Hertz was interested in the influence of the local elastic deformation of the contacting lenses when subjected to a normal contact force. The questions he poses in the introduction to his now-famous paper [START_REF] Hertz | Über die Berührung fester elastischer Körper[END_REF][START_REF] Hertz | Miscellaneous papers[END_REF]:

"What surface is it, of which the surface of pressure forms an infinitesimal part? What is the form, and what is the absolute magnitude of the curve of pressure? How is the normal pressure distributed over the surface of pressure?", give us an insight into his mindset while working on this problem. Hertz subsequently proceeds to address these questions by drawing on his background in electromagnetism and optical theory and thorough mathematical proofs. The main results from Hertz' classical theory are presented here. The Hertzian contact solution begins with the following assumptions:

• The surfaces in contact are perfectly smooth and non-conforming, with no friction and only a normal pressure acting between them;

• The surfaces are homogeneous and isotropic, subjected to linear kinematic equations and following the linear elastic material law;

• The two surfaces can be represented using second-order polynomials locally;

• The dimensions of the contact area are very small as compared to the contacting surfaces and thus may be considered locally to be the case of contact between two semi-infinite half-spaces.

Two different mathematical concepts are then combined [START_REF] Johnson | One hundred years of Hertz contact[END_REF]:

1. The geometry of two curved surfaces which touch without deformation;

2. The theory of potential applied to an elastic half-space bounded by a plane surface.

As the two surfaces near the vicinity of the contact area can be represented using second-order polynomials, the geometric separation between the two surfaces in the x-y plane can be expressed as:

h(x, y) ≈ Ax 2 + By 2 , (1.24)
where A and B represent the relative curvature of the two surfaces in the longitudinal and lateral directions. If the principal axes of both surfaces are aligned, the longitudinal A and lateral B curvatures are given by:

A = 1 2R x1 + 1 2R x2 , B = 1 2R y1 + 1 2R y2 , (1.25)
with R x1,2 and R y1,2 the radii of the two contacting surfaces in the x and y directions, respectively. The contact and separation conditions are given by Equation (1.16) and Equation (1.17), respectively, which have been defined in previous sections. Hertz hypothesised that the case of non-conformal contact allows restricting ourselves to the area which is very close to the point of contact, as the local stresses are extremely great compared with those occurring elsewhere in the contacting bodies. The strains in this zone are also assumed to be sufficiently small to lie within the scope of the linear theory of elasticity. The stresses and surface deformations can subsequently be found, considering each body as an elastic half-space bounded by the plane surface z = 0.

At the time when Hertz was working on his research, the theory of potential had already been used to assess the deformation of semi-infinite half-spaces, eminently by Boussinesq [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF] for the case of point loading. Hertz conjectured that the contact area is elliptical with semi-axes a and b, perhaps influenced by his observations of interference fringes at the contact of cylindrical lenses. He then proceeded to draw on the theory of potential to show that an elliptical pressure distribution of the form:

p n (x, y) = p 0 1 - x a 2 - y b 2 , (1.26)
acting on the aforementioned elliptical contact area induces a normal elliptical deformation given as:

u z = (L -M x 2 -N y 2 ) πE ⋆ , (1.27)
where L, M , and N are functions of the eccentricity and the maximum normal pressure p 0 . E ⋆ denotes the combined Young's modulus of elasticity and is given as:

1 E ⋆ = 1 -ν 2 1 E 1 + 1 -ν 2 2 E 2 .
(1.28)

For bodies with the same material properties such that E 1 = E 2 = E and ν 1 = ν 2 = ν, the combined Young's modulus is defined as:

E ⋆ = E 2(1 -ν 2 )
.

(1.29)

An appropriate choice of values for a, b and δ in Equation (1.27) will then satisfy the contact condition given by Equation (1.16). The expressions for these parameters for different profiles of contacting bodies have been detailed in various literature. The interested reader may refer to [START_REF] Johnson | Contact mechanics[END_REF]Chap. 4] for in-depth derivations. Here, we proceed to give the final results for a selected few contact cases.

Line contact of cylindrical bodies

For the case of line contact between two cylindrical bodies with their axes parallel to the y-axis pressed together by a vertical force N ′ per unit length, the contact area will be a strip of width 2a such that:

a = 4N ′ R eq πE ⋆ , (1.30)
p 0 = N ′ E ⋆ πR eq , (1.31)
where 1/R eq is the relative curvature given as:

1 R eq = 1 R 1 + 1 R 2 .
(1.32)

1.4.1.

Point contact of spheres

For point contact between two solids of revolution, with

R x1 = R y1 = R 1 and R x2 = R y2 = R 2 as shown in Figure 1.
13, the contact area will be circular with radius a such that:

a = 3 3N R eq 4E ⋆ , (1.33)
δ = 3 9N 2 16R eq E ⋆2 , (1.34)
p 0 = 3 6N E ⋆2 π 3 R 2 eq , (1.35)
where N is the applied vertical load.

Contact between general profiles

In the case of contact between general profiles, where the separation is given by Equation (1.24), the expressions for the different contact parameters involve elliptic integrals whose values must be found from tables. These expressions are generally pre-tabulated for engineering applications by replacing the elliptic integrals with equivalent algebraic equations in terms of the curvature ratio (A/B): where m, n, and r are Hertz' coefficients expressed as functions of the angle θ, as shown in Table 1.1:

a = m 3 3 2 N 1 2E ⋆ 1 A + B , (1.36) b = n 3 3 2 N 1 2E ⋆ 1 A + B , (1.37)
δ = r 3 3 2 N 1 2E ⋆ 2 (A + B) , (1.38)
cos θ = |A/B -1| A/B + 1 . (1.39)
The maximum normal pressure p 0 is calculated using:

p 0 = 3N 2πab . (1.40)
For wheel-rail contact, the longitudinal A and lateral B curvatures from Equation (1.25) can be rewritten as:

A = 1 2R wx = cos γ 2R , (1.41)
B = 1 2R wy + 1 2R ry , (1.42)
where the subscripts w and r denote the wheel and the rail, respectively, and R is the nominal rolling radius of the wheel. The term 1/2R rx in the longitudinal A curvature vanishes as the radius of the rail in the longitudinal direction is assumed to be infinite. For Equation (1.39) to be valid, only positive values of the curvature ratio A/B can be considered. A negative value for the curvature ratio would imply a situation which is not geometrically possible. In wheel-rail contact, the A curvature is always positive, as evident from Equation (1.41), which then implies that the B curvature should also remain positive for the Hertzian theory to remain valid.

Virtual penetration methods for non-Hertzian contact

Although Hertzian and multi-Hertzian methods are commonly used in MBS packages, a non-Hertzian contact model is required when the contact patch is not elliptical [START_REF] Liu | Comparison of wheel-rail contact models in the context of multibody system simulation: Hertzian versus non-hertzian[END_REF]. In this section, the non-Hertzian contact is addressed using virtual penetration (VP) methods. VP methods are based on the idea that the actual contact area can be approximated using the area in which two surfaces touching each other at the geometrical point of contact would intersect, if they are shifted towards each other by a distance δ. The geometrical point of contact is defined as the point where the two surfaces touch without transmitting any load. As the bodies cannot penetrate each other in reality, local surface deformations arise and the contact area is formed. It has been shown by Hashemi and Paul [START_REF] Hashemi | Contact stresses on bodies with arbitrary geometry: Applications to wheels and rails[END_REF] that for the same δ, the zone of intersection or the interpenetration region encircles this actual contact area. This is true when the influence function describing normal deflection of bodies is unidirectional, as in the case of Boussinesq's influence function [START_REF] Piotrowski | Wheel-rail contact models for vehicle system dynamics including multi-point contact[END_REF]. Not only the area of the interpenetration region is larger, but in a Hertzian case the aspect ratio of the contact ellipse is not the same as the one given by Hertz' theory.

In Equation (1.13), the surface deformation is totally neglected, and is instead implicitly considered by reducing the rigid body approach δ to a fraction of it. Thus, the deformed distance is now approximated as:

e z (x, y) ≈ h(x, y) -ϵδ , (1.43)
Scaling factors and shape correction strategies are subsequently used to prescribe a virtual penetration ϵδ such that the interpenetration zone approaches the real contact area. The interpenetration g(y) is defined as:

g(y) = ϵδ -h(y) if h(y) ≤ ϵδ , 0 if h(y) > ϵδ .
(1.44)

A point at a given y coordinate is said to be in contact if g(y) is greater than zero. The technique used for determining ϵ is what broadly differentiates the different virtual penetration approaches available in the literature [START_REF] Sichani | Comparison of non-elliptic contact models: Towards fast and accurate modelling of wheel-rail contact[END_REF].

The normal stress distribution in the direction of rolling is also assumed to be elliptical in such methods in order to be consistent with Hertz' theory. In the following sections, some available VP methods are outlined briefly, with a focus on the Kik and Piotrowski and the semi-Hertzian method, which are used later in the case studies of Chapter 4.

Kik and Piotrowski's method

In Kik and Piotrowski's (KP) method [START_REF] Piotrowski | Wheel-rail contact models for vehicle system dynamics including multi-point contact[END_REF][START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF][START_REF] Kik | A fast, approximate method to calculate normal load at contact between wheel and rail and creep forces during rolling[END_REF], the scaling factor is taken as a constant value with ϵ = 0.55. This value is obtained heuristically by carrying out several realistic numerical simulations using CONTACT. A similar value (ϵ = 0.5) was also found to give a good approximation of the contact area by Vohla [START_REF] Vohla | Werkzeuge zur realitätsnahen simulation der laufdynamik von schienenfahrzeugen[END_REF].

The origin in KP method is located at the geometrical contact point and contact conditions are satisfied only at the origin. The normal pressure distribution over the contact patch is assumed to be semi-ellipsoidal, such that:

p n (x, y) = p 0 a(0) a 2 (y) -x 2 , (1.45)
where p 0 is the maximum normal pressure. The limits of the interpenetration region in the longitudinal direction, used to define the contact area are calculated as: a(y) ≈ 2Rg(y) .

(1.46)

Equation (1.45) and Equation (1.46) imply that the normal pressure is always maximum at the point of geometrical contact, which is not necessarily true in the case of non-Hertzian contact [START_REF] Sichani | Comparison of non-elliptic contact models: Towards fast and accurate modelling of wheel-rail contact[END_REF]. Using Boussinesq's theory from Equation (1.12), the normal displacement at the origin is described by the integral:

u z (0, 0) = 1 -ν 2 πE p 0 a(0) C a 2 (y) -x 2 x 2 + y 2 dxdy . (1.47)
with u z (0, 0) = δ/2, as evident from Figure 1.12 and Equation (1.13) at the geometrical point of contact. The maximum pressure can now be numerically evaluated using the expression:

p 0 a(0) = δ πE 2(1 -ν 2 ) C a 2 (y) -x 2 x 2 + y 2 dxdy -1
.

(1.48)

The normal contact force is evaluated by integrating Equation (1.45) over the contact area. In a Hertzian case, the curvatures remain constant, and the curvature ratio can be calculated as:

A B = W L 2 , (1.49)
where L and W denote the length and the width of the uncorrected contact patch, as shown in Figure 1.14. These length and width values are in turn used to introduce the shape correction. The corrected length and width of the contact area are defined as: where β 0 is the aspect ratio a/b of the ellipse with curvature ratio A/B as given in Table 1.1.

W c = LW β 0 , (1.50)
L c = β 0 LW , (1.51)
The width of the resulting contact patch is corrected by introducing a "stretching/contracting" of the separation function across lateral y coordinate. This correction is performed such that the geometrical point of contact stays at its place. The length of the contact patch on the other hand is corrected by the modification of the Hertzian radius of curvature in the longitudinal direction [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF]. Correcting the contact patch using such an approach implies that each contact patch must be treated separately. The KP method has been improved recently by using a different strategy for the shape correction [START_REF] Sun | A robust non-Hertzian contact method for wheel-rail normal contact analysis[END_REF]. The method has also been improved to take into account the effects of the wheelset yaw, which leads to asymmetric contact patches [START_REF] Liu | A non-Hertzian method for solving wheel-rail normal contact problem taking into account the effect of yaw[END_REF].

Linder's method

Linder's method [START_REF] Piotrowski | Wheel-rail contact models for vehicle system dynamics including multi-point contact[END_REF][START_REF] Linder | Verschleiss von eisenbahnrädern mit unrundheiten[END_REF] uses the same scaling factor as the KP method, with ϵ = 0.55. The leading and the trailing coordinates of the interpenetration region are similarly computed as: a(y) ≈ 2Rg(y) .

(1.52)

Each pair of points [±a(y), y] is assumed to lie on an ellipse, and the lateral semi-axis of each ellipse is assumed to be the same and equal to:

b e = W 2 . (1.53)
The longitudinal semi-axis at each y co-ordinate is then calculated as:

a e (y) = g(y) b e b 2 e -y 2
.

(1.54)

The penetration δ e (y) is calculated geometrically for each ellipse using the interpenetration function [START_REF] Piotrowski | Wheel-rail contact models for vehicle system dynamics including multi-point contact[END_REF]. With a e (y), b e , and δ e (y), it is possible to calculate the normal load N (y), and finally the maximum pressure at each lateral co-ordinate is given as:

p 0 (y) = 3N (y) 2πa e (y)b e 1 - y b e 2 .
(1.55)

The semi-Hertzian method

The semi-Hertzian or STRIPES method proposed by Ayasse and Chollet [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF] can trace its origins back to the same time as when the KP method was being developed in the mid-90s, although it was only officially published in 2005.

Rather than having a fixed value, the scaling parameter ϵ in the STRIPES method may be defined in terms of the curvature ratio. In a Hertzian case, g(b) = 0 as seen from Figure 1.15. Furthermore, using h(b) = Bb 2 from Equation (1.24), we now have:

ϵ = Bb 2 δ = n 2 r 1 1 + A/B , (1.56)
with the semi-axis b defined using Equation (1.37). In MBS simulations, when the resultant normal force N is the known parameter, the term ϵδ rather than δ is the unknown quantity which is evaluated iteratively.

The STRIPES method differs from other VP methods in the fact that it evaluates the normal contact force N by introducing the idea of a local contact stiffness k l . The contact stiffness is independent of the contact patch dimensions for a Hertzian ellipse with semi-axes a and b, as shown in Figure 1.15. The normal force is assumed to be proportional to the product of this contact stiffness with the virtual interpenetration. The virtual interpenetration in the Hertzian case can be written as:

g(y) = B(b 2 -y 2 ) , (1.57)
with the normal force N :

N = k l b -b
g(y)dy .

(1.58)

Using Equation (1.37), Equation (1.57) and Equation (1.58), the contact stiffness can be evaluated as:

k l = E 2(1 -ν 2 ) 1 + A/B n 3 . (1.59)
For non-Hertzian contact, each strip is assigned a local stiffness which depends on the local curvature values:

k(y) = k l δy = E 2(1 -ν 2 ) 1 + A(y)/B(y) n 3 (y) δy , (1.60)
where δy is the width of each strip. The normal contact force is finally evaluated as the sum of the normal force per strip for all strips where g(y) is greater than zero:

N = N (y) = k(y)g(y) .
(1.61) Figure 1.16: Determination of the virtual interpenetration g(y) for wheel-rail profiles in STRIPES

The half-lengths of the contact area a(y) can be subsequently evaluated to compute the contact pressure distributions. In order to be consistent with Hertz' theory, the shape of contact area is corrected using the longitudinal and lateral curvatures A(y) and B(y), respectively. The chosen procedure is to correct the A(y) curvature only, which corresponds to the rolling direction. The corrected curvature A c (y) is given as:

A c (y) A(y) = [n(y)/m(y)] 2 A(y)/B(y) . (1.62)
The half-length of the contact patch in the rolling direction a(y) is calculated using the expression:

a(y) = g(y) A c (y) . (1.63)
In a Hertzian case, replacing A c (y) with A in Equation (1.63) would give the same semi-axis a as Equation (1.36). A second possibility involving the correction of both the longitudinal and lateral curvatures has also been presented in [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF]. For small values of the aspect ratio a/b, this correction leads to inconsistent results, and is thus not taken into consideration. As discussed in Section 1.4.1.3, the Hertzian parameters are not defined when the lateral B(y) curvature is negative. Consequently, B(y) is enforced to be always greater than zero by replacing any negative values by infinitesimal positive values. Moreover, in the case where the curvature B(y) is discontinuous, the normal forces will also be discontinuous. The KP method avoids this issue, as only the geometry of the profile is used to assess the normal contact force. This is more robust in dynamic simulations as it avoids the computation of the curvatures. In STRIPES, a further step involves the smoothing of the B(y) curvature before the calculation. The smoothing filter is based on the Boussinesq approach, where the characteristic width depends on the curvature value [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF]. The smoothed curvature B s (y) is approximated as:

dB s (y) dy = B(y) -B s (y) 2 3 L o (y) , (1.64)
where L o is the characteristic length:

L o (y) = 3 3 2 N 1 -ν 2 E 1 A(y) + B(y) . (1.65)
Assuming that the pressure distribution is elliptic in the rolling x direction, the normal pressure distribution is given as:

p n (x, y) = p 0 (y) 1 - x a(y) 2 ,
(1.66)

with the maximum pressure at the centre of each strip:

p 0 (y) = 2 π N (y) a(y)δy . (1.67)
As opposed to the KP method, the geometrical point of contact does not play an important role as far as the normal contact results are concerned, and the maximum pressure value is independent of its location. Additionally, this also implies that multiple contact points can be handled with relative ease using Equation (1.44) and need not be separated as required in the KP method. The whole wheel-rail normal contact problem, including multiple contact, can be handled simultaneously by introducing the vertical interpenetration δtz1 , as shown Figure 1.16. For a given position of the wheel and a given geometry of the rail, z(y) is defined as the vertical separation between the wheel and the rail. z(y) is linked to h(y) through a projection using the contact angle γ(y): h(y) = z(y) cos γ(y) .

(1.68)

Equation (1.44) defining the virtual penetration g(y) is similarly amended by using a projection:

g(y) = [δtz -z(y)] cos γ(y) if z(y) ≤ δtz , 0 if z(y) > δtz .
(1.69)

ANALYN

The method ANALYN proposed by Sichani et al. [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF] attempts to enhance the virtual penetration approaches by introducing an approximation for the surface deformation. The principle idea of this approximation is based on the observation that there is a similarity between the elastic deformation on the surface and the separation in a Hertzian contact. The expression for the interpenetration function, in this case, is given as:

g(y) = δ -h(y) -β(y)h(y) , (1.70)
with the additional term β(y)h(y) taking into account the influence of the surface deformation. Taking only positive values of g(y), the half-length of the contact area a(y) is:

a(y) = g(y) (1 + α(y))A(y)
.

(1.71)

The coefficients α(y) and β(y) are obtained using the results from Hertz' theory (see Appendix in [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF]):

α(y) = r(y) m 2 (y) 1 + B(y) A(y) -1 , (1.72)
β(y) = r(y) n 2 (y) 1 + A(y) B(y) -1 . (1.73)
The maximum pressure value p 0 (y) at the centre of the contact patch is:

p 0 (y) = E π(1 -ν 2 ) 1 n(y)r(y) g(y) a(y) , (1.74)
and the pressure distribution p n (x, y) in elliptic in the rolling direction. As in STRIPES, the Hertzian parameters are limited to positive values of the B(y) curvature. The non-positive part of the curvatures is therefore replaced by a polynomial function heuristically [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF]. This inhibits straightforward implementation in MBS codes, as seen in the recent switches and crossing benchmark where only Hertz, KP and STRIPES methods were used [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF].

Tangential contact problem in Hertzian conditions

1.5.1 Carter's theory

Frederick William Carter's work on the contact of locomotive wheels and rails in 1926 [START_REF] Carter | On the action of a locomotive driving wheel[END_REF] was the first to introduce the idea that the driving wheel subjected to a tangential force in any direction, besides rolling, also creeps slowly in the general direction of the force. He defined this "creepage" as the ratio of the distance gained by one surface over the other in the longitudinal direction, to the distance traversed. The railway track is assumed to be an elastic half-space, and the wheel an elastic cylinder. The two bodies are also considered quasi-identical. With the further consideration that the wheel is large compared to the contact area, the problem is one of an infinite elastic medium bounded by a plane, under the influence of a certain local distribution of pressure and tangential traction. Using Professor Love's mathematical analyses [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF] on this subject, Carter then proceeded to establish a relation between the longitudinal creepage and the tangential contact force. Carter's theory is significant from the point of view that analytical solutions to the rolling contact problem can be obtained easily. A similar problem was solved almost at the same time by Hans Fromm in his doctoral thesis [START_REF] Fromm | Berechnung des schlupfes beim rollen deformierbarer scheiben[END_REF]. However, Fromm assumed that the two bodies were elastic cylinders and not half-spaces. The contact area in the rolling direction is split into an adhesion (or stick) zone near the leading edge, where the creepage is constant; and a slip zone near the trailing edge, where the tangential stress is bound by Coulomb's law of friction. Denoting the coefficient of friction as µ, the tangential traction over the stick zone can be obtained by subtracting an elliptic traction distribution from the traction bound:

p x = p ′ x + p ′′ x , (1.75)
where p ′ x is the elliptic traction bound over the contact strip of length 2a:

p ′ x = µp n = µp 0 1 - x a 2 , -a ≤ x ≤ a , (1.76)
and p ′′ x is the elliptic traction bound over the adhesion zone when the slip length is 2d: (1.78)

p ′′ x = - c a µp 0 1 - x -d c
The longitudinal stress direction opposes the direction of the slip. The tangential stress distribution over the entire contact strip in Carter's theory is illustrated in Figure 1.17. The width of the slip zone is determined by the magnitude of the creep force F x such that [18, p. 253]:

d a = 1 - c a = 1 -1 - F x µN , (1.79)
where N is the normal force. The resultant tangential strain in the strip and Hertz' relations for the maximum pressure p 0 are subsequently used to establish the relationship between the creep ratio ν x , and the tangential force [18, p. 253]: where R eq is defined using Equation (1.32). The relationship between the creep force F x and the creep ratio ν x , as shown in Figure 1.18, is known as a "creep curve". Finding the tangential traction by assuming full slip solution over the entire contact area and subtracting a traction from the adhesion area, as done in Equation (1.75), was also subsequently used by other researchers, including in the theory of Johnson [START_REF] Johnson | The effect of a tangential contact force upon the rolling motion of an elastic sphere on a plane[END_REF], the Vermeulen-Johnson [START_REF] Vermeulen | Contact of non-spherical elastic bodies transmitting tangential forces[END_REF] solution, and the strip theories of Haines-Ollerton [START_REF] Haines | Contact stress distributions on elliptical contact surfaces subjected to radial and tangential forces[END_REF], Halling [START_REF] Halling | Microslip between a rolling element and its track arising from geometric conformity and applied surface tractions[END_REF], and Kalker [START_REF] Kalker | A strip theory for rolling with slip and spin[END_REF]. However, neglecting the effect of the lateral and spin creepages, as well the restriction to a one-dimensional analysis, implies that the Carter theory is largely limited in terms of its application to the analysis of contact problems involving complex railroad vehicles.

ν x = - µa R eq 1 -1 - F x µN , (1.80)

Johnson's theories

In 1958, K. L. Johnson published two articles dealing with the problem of three-dimensional elastic rolling contact of quasi-identical bodies. The first paper addresses the case of rolling with pure creepage where the spin is assumed to be absent [START_REF] Johnson | The effect of a tangential contact force upon the rolling motion of an elastic sphere on a plane[END_REF]. The contact area is circular, which was later extended to ellipses by Vermeulen and Johnson [START_REF] Vermeulen | Contact of non-spherical elastic bodies transmitting tangential forces[END_REF]. The contact area is divided into two regions as done previously, the slip zone where the tangential traction reaches the traction bound, and the adhesion or stick zone where the slip vanishes. The adhesion zone is assumed to be an ellipse Figure 1.19: Stick-slip zone division according to Johnson [START_REF] Johnson | The effect of a tangential contact force upon the rolling motion of an elastic sphere on a plane[END_REF], and Vermeulen and Johnson [START_REF] Vermeulen | Contact of non-spherical elastic bodies transmitting tangential forces[END_REF] theories with the same aspect ratio as the contact area, meeting the contact ellipse at the leading edge as shown in Figure 1.19. The adhesion region does not include the leading edge of the contact ellipse, and the slip and tangential traction are in the same direction in the small area highlighted in red near the leading edge. Thus, Equation (1.19) is not satisfied, and the solution is only approximate. The input of the method are Ft , the non-saturated creep forces expressed as linear functions of the longitudinal and lateral creepages. An exact derivation of Ft has been developed by Kalker, as will be seen in the next section. Similar to Carter's theory, the resultant tangential traction is obtained by considering the sliding traction over the entire contact area and subtracting a second traction acting only over the adhesion area. Consequently, the no-slip condition is satisfied in the adhesion area, and the tangential traction in the slip area equals the traction bound. The relations for the creep force acting on the elliptical contact area shown in Figure 1.19 from the Vermeulen-Johnson theory are given as [21, p. 75]:

p t = µp 0 Ft Ft   1 - x a 2 - y b 2 -1 - x ′ a ′ 2 - y ′ b ′ 2   if x ∈ A , (1.81)
p t = µp 0 Ft Ft 1 - x a 2 - y b 2 if x ∈ S , (1.82)
where

p 0 = 3N 2πab , (1.83)
Ft = || Ft || , FT t = [ Fx , Fy ] . (1.84)
The semi-axis of the adhesion ellipse a ′ is defined such that:

x ′ = x + a ′ -a . (1.85)
From the Vermeulen-Johnson theory, we have:

a ′ a = 3 1 - F t µN , (1.86)
where the relationship between F t and Ft is given by a cubic saturation law:

F t =      µN Ft µN -1 3 Ft µN 2 + 1 27 Ft µN 3 if Ft ≤ 3µN , µN if Ft > 3µN .
(1.87)

Johnson's second paper [START_REF] Johnson | The effect of spin upon the rolling motion of an elastic sphere on a plane[END_REF] is significant because it introduced the idea of a spin creepage that tends to "twist" the contact interface for the first time. Johnson showed that in the case of pure spin, with the longitudinal and lateral creepage equal to zero, the spin gives rise to a lateral creep force and a moment about the z-axis. This is a no-slip theory, unlike partial slip in the precedent case, where the coefficient of friction is assumed to be infinitely large such that the slip vanishes in the entire contact area. The first to apply the no-slip theory in three-dimensional frictional contact was Mindlin [START_REF] Mindlin | Compliance of elastic bodies in contact[END_REF], who solved the no-slip problem of two quasi-identical Hertzian bodies pressed together. The two bodies are subsequently subjected to a displacement, and a rotation about the z-axis relative to each other, without rolling.

The no-slip theory is also called the linear theory of rolling as the longitudinal and lateral creep forces, as well as the creep moment, are linear functions of the creepages and the spin. In [START_REF] Johnson | The effect of spin upon the rolling motion of an elastic sphere on a plane[END_REF], Johnson considers the case where F y = 0. For a circular contact area, the relationship between the lateral creepage and the spin is given as:

ν y = - 2 3 (2 -ν) (3 -2ν) aφ , (1.88)
The moment about the z-axis is shown to be:

M z = 32 9 (2 -ν) (3 -2ν) Ga 4 φ . (1.89)
Kalker has compared Johnson's linear spin theory and Vermeulen and Johnson's no-spin theory with the exact values in [21, p. 66], where he presents the relative error between the exact and approximated coefficients of proportionality for the creepages and spin.

Kalker's theories

Along with Johnson, Professor Joost Kalker is perhaps one of the most significant contributors to the theories of rolling contact mechanics in the last fifty years or so. Many of his theories developed in the latter half of the twentieth century are still in use today as references for new rolling contact models. A good summary of all of Kalker's theories can be found in his monograph [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF] or in several review publications [START_REF] Kalker | Survey of wheel-rail rolling contact theory[END_REF][START_REF] Vollebregt | Survey of programs on contact mechanics developed by[END_REF][START_REF] Zaazaa | Review of Joost Kalker's wheel-rail contact theories and their implementation in multibody codes[END_REF].

Strip theory

The original strip theory for steady-state rolling by Haines and Ollerton in 1963 [START_REF] Haines | Contact stress distributions on elliptical contact surfaces subjected to radial and tangential forces[END_REF] was extended to include the effect of lateral and small spin creepage by Kalker in 1967 [START_REF] Kalker | A strip theory for rolling with slip and spin[END_REF]. The shortcomings of the Vermeulen-Johnson theory, as highlighted in Figure 1.19, are avoided as the contact area is divided into thin strips parallel to the rolling x-axis. The twodimensional Carter's solution is then applied to each strip, considering it to be independent of the influence of its neighbour, as shown in Figure 1.20.

The strip theory is confined to slender ellipses, i.e., when the contact area is considerably wider (twice or three times [START_REF] Kalker | Survey of wheel-rail rolling contact theory[END_REF]) in the lateral y direction than in the rolling x direction, and to small spin values. The failure to approximate the three-dimensional creep-force law meant that the strip theory was superseded by other relevant theories shortly after it had been proposed. The strip theory has recently been revisited by Sichani et al. in the method FaStrip [START_REF] Sichani | An alternative to FASTSIM for tangential solution of the wheel-rail contact[END_REF] which will be addressed in Section 1.5.4. The linear theory was based on an idea by de Pater in 1956 which was then developed by Kalker during his doctoral thesis [START_REF] Kalker | On the rolling contact of two elastic bodies in the presence of dry friction[END_REF]. Here, Kalker assumes that the true slip vanishes at the leading edge and that the traction must be zero when the particles enter the contact area [START_REF] Kalker | Survey of wheel-rail rolling contact theory[END_REF]. The linear theory gets its name from the fact that the creep-force law curve here is linear, assuming no slip occurs over the contact area. Kalker's non-dimensional coefficients c 11 , c 22 , c 23 , and c 33 are also introduced, which represent the constants of proportionality between the creepages and the creep forces and moment.

The creep forces and the spin moment can be expressed in matrix form as:

   F x F y M z    = -Gab    c 11 0 0 0 c 22 √ abc 23 0 - √ abc 23 abc 33       ν x ν y φ    . (1.90)
The longitudinal creep force F x is independent of the lateral creepage ν y and the spin φ, while the lateral creep force F y and φ are coupled, which also corresponds to the definition of spin by Johnson [START_REF] Johnson | The effect of spin upon the rolling motion of an elastic sphere on a plane[END_REF]. Kalker's coefficients are functions of the Poisson's ratio ν and the aspect ratio a/b of the contact ellipse, as shown in Table 1.2. The values of these coefficients for different ν and a/b have been computed by Kalker using the complete theory in [START_REF] Kalker | On the rolling contact of two elastic bodies in the presence of dry friction[END_REF] and compared to the approximate values found by Vermeulen and Johnson [START_REF] Vermeulen | Contact of non-spherical elastic bodies transmitting tangential forces[END_REF]. Kalker also subsequently provided empirical expressions for the creepage coefficients in [START_REF] Kalker | The tangential force transmitted by two elastic bodies rolling over each other with pure creepage[END_REF].

The linear theory is applicable in the case of quasi-identical bodies. The expressions for the combined modulus of rigidity G, and the combined Poisson's ratio ν are defined using Equation (1.8) and Equation (1.9). When it was put forth in the 1970s, the interest of this simplified approach was to provide a stiffness coefficient for evaluating the critical speed of a linearised bogie model, with the contact assumed to be of constant form under a constant load. The most important contribution by Kalker was the rigorous calculation of the c ij coefficients, which were also verified experimentally for the wheel-rail contact cases on test-benches [START_REF] Chollet | Etude en similitude mecanique des efforts tangents au contact roue-rail[END_REF], as this was an important aspect in determining the risk of instability and derailment of railway vehicles. With respect to more recent theories, the main drawback of the linear theory is that it does not take into account saturation, and hence is only applicable in the case of infinitesimal creepages.

Simplified theory: The FASTSIM algorithm

Kalker's simplified theory [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF] is perhaps the one that is most commonly used in the railway industry codes for evaluating the tangential contact results. The idea of the simplified theory is to replace the elastic body with a set of springs. The tangential surface tractions and the tangential surface displacements at a given point are assumed to be linearly proportional through three flexibility parameters in the adhesion region, as shown in Figure 1.21, and limited by a parabolic traction bound according to Coulomb's theory. This parabolic traction bound is shown to provide a more accurate division of the stick-slip zone when compared to an elliptic traction bound [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF][START_REF] Vollebregt | FASTSIM2: A second-order accurate frictional rolling contact algorithm[END_REF]. The flexibility parameters depend on the creepage and spin coefficients c 11 , c 22 , and c 23 of the linear theory. The tangent problem of the simplified theory is solved iteratively using the FASTSIM algorithm developed by Kalker [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF], which will be addressed in Chapter 3.

The FaStrip algorithm

The FaStrip [START_REF] Sichani | An alternative to FASTSIM for tangential solution of the wheel-rail contact[END_REF] algorithm provides an alternative to FASTSIM, and is based on the strip theory for three-dimensional rolling contact [START_REF] Kalker | A strip theory for rolling with slip and spin[END_REF]. This is combined with a FASTSIM type algorithm to determine the tangential stress directions. The nonlinear stress distribution in the stick area, and

                                                 0.0 π 2 /4(1-ν) π 2 /4 π √ g{1+ν( 1 2 Λ+ln 4-5)} 3(1-ν) π 2 /
                                                     1.0 3.
(Λ-2ν)g (1+ 3-ln 4 Λ-2ν ) 2π g [1+ (1-ν)(3-ln 4) (1-ν)Λ+2ν ] (1-ν)Λ+2ν 2π 3g √ g[(1-ν)Λ-2+4ν] π 4 [1-ν(Λ-2) (1-ν)Λ-2+4ν ] Λ = ln(16/g 2 ); g = min(a/b; b/a); ln 4 = 1.386
Table 1.2: Creepage coefficients c ij from Kalker's linear theory of rolling contact for elliptical contact areas. From [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF] Figure 1.21: Stick-slip zone division and the traction distribution over each strip according to the simplified theory, with φ = 0 the use of an elliptic traction bound in the slip area in the strip theory, results in better agreement with Kalker's complete theory [START_REF] Kalker | Rolling contact phenomena: Linear elasticity[END_REF] as compared to the simplified theory [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF]. Nevertheless, the accuracy of the strip theory is limited to slender contact ellipses, as well as to small values of spin. These drawbacks are addressed in the method FaStrip to extend the strip theory to more general contact cases. The strip theory incorporates Carter's exact solution for two-dimensional contact [START_REF] Carter | On the action of a locomotive driving wheel[END_REF], and the tangential tractions p x and p y in the adhesion region are given as:

p x (x, y) = µp 0 a κ a(y) 2 -x 2 -κ ′ [a(y) -d(y)] 2 -[x -d(y)] 2 , (1.91) p y (x, y) = µp 0 a λ a(y) 2 -x 2 -λ ′ [a(y) -d(y)] 2 -[x -d(y)] 2 , (1.92)
where µ is the coefficient of friction, a(y) is the half-length of the contact patch, d(y) is the half-length of the slip area, and κ, κ', λ, and λ' are the parameters dependent on the creepage and spin values:

a(y) = a 1 - y b 2 , (1.93)
d(y) = a η 2 + (1 -ψ 2 )(ξ -ψy a ) 2 + ηψ (1 -ν)(1 -ψ 2 ) , (1.94)
κ = ξ -ψy a ξ -ψy a 2 + η + ψd(y) a 2 = κ ′ , (1.95)
λ = η + ψd(y) a 2 ξ -ψy a 2 + η + ψd(y) a 2 = λ ′ + ψ , (1.96)
where ν is the Poisson's ratio, and ξ, η, and ψ are non-dimensional parameters for the creepages and spin, respectively:

ξ = - G 2µp 0 ν x , (1.97)
η = - G 2µp 0 (1 -ν)ν y , (1.98)
ψ = - G 2µp 0 aφ . (1.99)
These expressions are found to give good results in Hertzian cases when the semi-axes ratio is small. To improve the accuracy for all contact patches, FaStrip imposes the assumptions of Kalker's linear theory [START_REF] Kalker | Rolling contact phenomena: Linear elasticity[END_REF] such that the integration of the tangential tractions results in the same creep forces [START_REF] Sichani | A fast wheel-rail contact model for application to damage analysis in vehicle dynamics simulation[END_REF]. This is achieved by multiplying the creepage and spin parameters ξ, η, and ψ with correctional terms such that:

ξ ′ = - 2 π 2 G(1 -ν) µp 0 c 11 ν x , (1.100)
η ′ = - 2 π 2 G(1 -ν) µp 0 c 22 ν y , (1.101
)

ψ ′ = - 3 2π G µp 0 √ abc 23 φ . (1.102)
The terms ξ, η, and ψ in Equations (1.94)-(1.96) are subsequently replaced by ξ ′ , η ′ , and ψ ′ , respectively. The stress directions in the slip area are evaluated using the FASTSIM algorithm [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF] detailed in Chapter 3, as it takes into account the effect of the local spin value. However, unlike FASTSIM, FaStrip uses an elliptic traction bound in the slip area. For strips covered entirely by the slip area, the three flexibility parameters from the simplified theory are replaced by a single equivalent flexibility parameter as suggested by Kalker [START_REF] Kalker | Rolling contact phenomena: Linear elasticity[END_REF] to improve the accuracy for higher values of creepages and spin [START_REF] Sichani | An alternative to FASTSIM for tangential solution of the wheel-rail contact[END_REF].

The original version of FaStrip expresses the tangential tractions in terms of a and b, the longitudinal and lateral semi-axes of the Hertzian contact patch, and the maximum normal pressure p 0 . In a later publication, these parameters are replaced with their local values a(y), b(y) and p 0 (y) to extend the approach to non-elliptical contact patches [START_REF] Sichani | A fast wheel-rail contact model for application to damage analysis in vehicle dynamics simulation[END_REF]. For the contact patch half-length a(y) and the maximum pressure distribution p 0 (y), the procedure is more or less straightforward, as these are the output quantities of any normal contact method. The assessment of b(y) is not as simple because a clear definition in non-Hertzian cases is not readily available. A possible solution is to assess b(y) using the relations from Hertz' theory locally [START_REF] Ayasse | Wheel-rail contact[END_REF]:

b(y) a(y) = n(y) m(y) . (1.103)

Heuristic solutions considering spin

Several heuristic solutions for the creep forces have been presented in the literature. In the following sections, two of these heuristic approaches, namely the one by Shen et al. [START_REF] Shen | A comparison of alternative creep force models for rail vehicle dynamic analysis[END_REF] and Polach's model [START_REF] Polach | A fast wheel-rail forces calculation computer code[END_REF], are outlined briefly.

Shen-Hedrick-Elkins model

The method of Shen-Hedrick-Elkins (SHE) [START_REF] Shen | A comparison of alternative creep force models for rail vehicle dynamic analysis[END_REF], is a combination of the linear theory of Kalker and the cubic saturation law of Vermeulen and Johnson. However, unlike the Vermeulen and Johnson approach [START_REF] Vermeulen | Contact of non-spherical elastic bodies transmitting tangential forces[END_REF], spin creepage is taken into account in the SHE method. The creep forces Fx and Fy are first evaluated using Equation (1.90) from the linear theory. Then using the cubic saturation law, the resultant tangential force Ft is limited by the nonlinear value F t as follows [START_REF] Shen | A comparison of alternative creep force models for rail vehicle dynamic analysis[END_REF]:

F t =      µN Ft µN -1 3 Ft µN 2 + 1 27 Ft µN 3 if Ft ≤ 3µN , µN if Ft > 3µN .
(1.104)

The nonlinear creep-force model is given by:

F x F y = ε Fx Fy , (1.105)
where ε = F t Ft .

(1.106)

The results using this heuristic model have been shown to be in good agreement with those from Kalker's simplified theory (FASTSIM), although it does lead to unsatisfactory results when the spin is significant.

Polach's model

In 1999, Polach presented a fast algorithm (as well as the computer code) for the computation of wheel-rail creep forces for a given Hertzian contact patch and the global creepages acting on it. Similar to the simplified theory, the tangential stress grows linearly over the adhesion area, starting from the leading edge, and is limited by an elliptic traction bound using Coulomb's law, as shown in Figure 1.22. According to these assumptions, the analytical expression of the tangential force in the absence of spin is given by [START_REF] Shabana | Railroad vehicle dynamics: a computational approach[END_REF][START_REF] Polach | A fast wheel-rail forces calculation computer code[END_REF]:

F t = 2µN π ε 1 + ε 2 + arctan ε , (1.107)
where ε is the gradient of the tangential stress in the adhesion area:

ε = π 4 GabC h µN ν c . (1.108)
The first term in Equation (1.107) roughly represents the contribution of the adhesion zone and tends to 0 if ε tends to infinity. The second term is associated with the slip zone and tends towards π/2 when ε tends towards infinity. If ε is small, we obtain Kalker's linear theory. C h is a constant that depends on Kalker's coefficients c ij such that:

C h = c 11 ν x ν t 2 + c 22 ν y ν t 2 .
(1.109) ν c represents the corrected resultant creepage that takes the influence of spin into account heuristically through a corrected lateral creepage ν yc computed using:

ν yc = ν y if |ν y + φa| ≤ |ν y | , ν y + φa if |ν y + φa| > |ν y | , (1.110)
with

ν c = ν 2 x + ν 2 yc , (1.111)
ν t = ν 2 x + ν 2 y .
(1.112) The contribution of the spin creepage to the lateral tangential force is given by:

F yc = 9 16 aµN K 1 + 6.3(1 -e -a/b ) , (1.113)
where K, δ and ε c are defined using the following expressions:

K = |ε c | - δ 3 3 + δ 2 2 - 1 6 + 1 3 (1 -δ 2 ) 1.5 , (1.114)
δ = ε 2 c -1 ε 2 c + 1 , (1.115
)

ε c = 8 3 Gb √ ab µN c 23 ν yc 1 + 6.3(1 -e -a/b ) . (1.116)
The longitudinal and lateral creep forces F x and F y , respectively, are finally given as:

F x = -F t ν x ν c , (1.117)
F y = -F t ν y ν c -F yc φ ν c .
(1.118)

Lookup tables

Lookup tables (LUT) offer the fastest and most straightforward manner of evaluating the creep forces. This approach consists of pre-calculating a table with a detailed program, and then interpolating the table values for the desired input conditions. Limited computing abilities in the early 1980s was partly the motivating factor for British Rail to develop one of the first LUT for evaluating wheel-rail contact forces, using the solutions provided by Kalker's DUVOROL program [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF]. Here, the aspect ratio of the ellipses, the longitudinal creepage, the lateral creepage and the spin were the varying input parameters used to construct the lookup table. Later, Kalker enlarged this table using the results from CONTACT [START_REF] Kalker | Book of tables for the Hertzian creep-force law[END_REF] with a total of 115,000 table entries. This table, called USETAB, was itself superseded by an enhanced version where a finer discretisation was used for the contact area [START_REF] Vollebregt | Survey of programs on contact mechanics developed by[END_REF]. USETAB has been shown to provide more precise results when compared to other simplified approaches for Hertzian cases, using Kalker's complete rolling contact theory as the reference [START_REF] Vollebregt | Assessing the accuracy of different simplified frictional rolling contact algorithms[END_REF]. The LUT presented so far are only applicable to Hertzian contact patches. For non-elliptical contact patches, an equivalent ellipse must first be evaluated. Piotrowski et al. [START_REF] Piotrowski | The Kalker book of tables for non-Hertzian contact of wheel and rail[END_REF] have developed the Kalker book of tables for non-Hertzian contact by the regularisation of non-elliptical contact patches to single double-elliptical contact (SDEC) patches. The number of entries to store, and the technique used to interpolate them, play an important role in deciding the accuracy of the LUT. Marques et al. [START_REF] Marques | On the generation of enhanced lookup tables for wheel-rail contact models[END_REF] have presented two procedures for the generation of enhanced LUT for wheel-rail contact models, where either the interpolation accuracy is improved for a similarly sized LUT, or the size of the LUT is condensed for a similar degree of accuracy. While LUT first made their appearance when the hardware capabilities of the available computers were largely restricted, they continue to be used in several MBS codes for their ease of implementation and accurate1 results for wheel-rail contact cases.

Conclusions

Safety, track fatigue analyses, and maintenance of railway vehicles are only some of the applications that highlight the importance of the choice of wheelrail contact model used in MBS simulation software. The complexity of different operations necessitates very different levels of modelling, and a simplified and coarse model may often be sufficient for some applications.

For online simulations where the calculations must be carried out in realtime over railway tracks which are several kilometres long, very approximate wheel-rail contact models are generally used. Other applications require more sophisticated models, as in the estimation of wear in urban rail networks, or the assessment of rolling contact fatigue. The shape of the contact area and the stress distribution within it categorise the different normal contact models available in the literature [START_REF] Piotrowski | Wheel-rail contact models for vehicle system dynamics including multi-point contact[END_REF][START_REF] Kalker | Survey of wheel-rail rolling contact theory[END_REF][START_REF] Knothe | History of wheel/rail contact mechanics: from Redtenbacher to Kalker[END_REF][START_REF] Meymand | A survey of wheel-rail contact models for rail vehicles[END_REF].

Most MBS codes use one of the different VP methods described in this chapter to evaluate the normal contact results. A proposed improvement to these models is through the introduction of an analytical approximation for the surface deformation in the ANALYN method [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF], although it does not seem to have been used in an MBS code yet. While existing methods reliably predict contact stresses in Hertzian conditions, the research of more accurate solutions is still an open field in the context of railway dynamics, where a compromise is required between the computing time and accuracy. Chapter 2 is dedicated to the development of a new semi-analytical approach for the normal contact problem, which aims to improve the prediction of normal contact stresses in non-Hertzian conditions.

Introduction

The following chapter outlines the development of a new semi-analytical boundary element method for determining the wheel-rail contact zone and the normal stress distribution within it. The chapter largely follows the same outline as in the original publication [START_REF] Qazi | A semi-analytical numerical method for modelling the normal wheel-rail contact[END_REF]. A reduced formulation of the contact between two elastic half-spaces is used in the first part together with a Hertzian model to solve the normal contact, following a semi-analytical methodology. The novel method is then validated using the results from CONTACT [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF] as a reference, with theoretical as well as wheel-rail profiles. A comparison of the obtained results and the required computational resources is also carried out using the STRIPES [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF] and ANALYN [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF] methods described in Chapter 1. Two alternative modelling strategies are proposed as a part of the perspectives, before the conclusions.

Methodology

Theoretical background

We consider the problem describing the stress field in a semi-infinite halfspace, subjected to a concentrated normal surface force, which was studied in detail by Boussinesq [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF]. As seen previously, for the half-space consideration to be valid locally, the contact area dimensions must be significantly smaller than the principal radii of the contacting surfaces. This is generally true when considering several common contact scenarios, including the tread contact between wheels and rails. The following assumptions are also made:

• The bodies are elastic, homogeneous and isotropic;

• The contacting surfaces have a continuous profile;

• The two bodies are quasi-identical, which is the case of steel-steel contact between the wheel and the rail.

These assumptions subsequently allow the normal and tangential problems to be solved sequentially. One possible solution is using the potential theory of Boussinesq [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF], described exhaustively in the book by Johnson [START_REF] Johnson | Contact mechanics[END_REF]. The normal elastic displacement u z (x, y) is related to the normal pressure distribution p n (x, y) by the integral equation [START_REF] Cesbron | Numerical and experimental study of multi-contact on an elastic half-space[END_REF] ∀(x, y) ∈ P, u z (x, y) E ⋆ is the combined Young's modulus of elasticity for the two contacting bodies described using Equation (1.29). With a defined rigid approach δ, the only geometrical input required to solve the contact problem is the separation between the two undeformed surfaces h(x, y), with the deformed distance e z (x, y):

= C p n (ξ, η)A zz (x, y; ξ, η)dξdη , (2.1)
e z (x, y) = h(x, y) + u z (x, y) -δ . (2.
3)

The conditions for contact and separation are as defined in Equation (1.16) and Equation (1.17):

∀(x, y) ∈ C : e z = 0, p n ≥ 0 , (2.4) ∀(x, y) / ∈ C : e z > 0, p n = 0 , (2.5)
which ensure that the normal pressure can only act within the contact zone, vanishing at the boundary of the contact area. The two bodies are also prevented from interpenetrating. The discretised form of the variational inequality eventually gives rise to a linear complementarity problem (LCP).

If the rigid approach δ is given, the problem to be solved can be completely described by Equation (2.1) and Equation (2.4).

If the total resultant force N is known in advance, the following equilibrium condition is added:

∀(x, y) ∈ C, N = C p n (ξ, η)dξdη .
(2.6)

Towards a reduced approach

The direct method to solve the contact problem is the matrix inversion method (MIM) found in [18, p. 144]. The potential contact area in the boundary element problem is overestimated in both the lateral and the rolling directions, and divided into a rectangular grid, wherein the contact constraints of Equation (2.4) and Equation (2.5) are evaluated.

The elements in which the pressure distribution has negative values are discarded, and the procedure is repeated until all the remaining elements satisfy the contact conditions. The elastic deformation is thus calculated at each iteration for each element within the potential contact area, resulting in a sizeable expenditure of computational resources. The drawback of using such an approach is clear, as the direct method does not scale to large-sized problems. Certain studies also use iterative techniques such as Gauss-Seidel, which however remain restricted to relatively small contact grids. Most recent works use either a conjugate gradient (CG) algorithm or a multi-grid (MG) strategy, combined with fast algorithms such as multi-level summation technique (MLMS) or fast Fourier transform (FFT) for equation solving [START_REF] Polonsky | A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques[END_REF][START_REF] Vollebregt | A new solver for the elastic normal contact problem using conjugate gradients, deflation, and an FFT-based preconditioner[END_REF]. The principle essence of these approaches remains the same: using different numerical techniques for a faster solution of the complete LCP.

In the proposed new approach, the pressure distribution is additionally assumed to be symmetric and elliptic about the x = 0 plane (i.e. in the rolling direction) and the potential area of contact is discretised only in the lateral y direction. A similar strategy for discretisation is also used by Reusner for the treatment of roller bearings in [START_REF] Reusner | Druckflächenbelastung und oberflächen-verschiebung im wälzkontakt von rotationskörpern[END_REF], and by Knothe and Le The for arbitrary elastic bodies in [START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with nonelliptical contact area[END_REF]. Along with the normal stress distribution acting on each strip, Reusner considers the lengths of the contact patch strips as additional unknowns, while Knothe and Le The attempt to further reduce the computational complexity by showing that the strip lengths in neighbouring elements should vary almost proportionally with the variations in the load or the deflection.

The novel approach presented here is to instead consider the contact Adapted from [START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with nonelliptical contact area[END_REF] patch boundaries as a quasi-known quantity dependent on the form of the normal stress distribution. The potential contact area is first divided into thin strips with the larger dimension in the rolling direction, as shown in Figure 2.2, from where the name MIM-1D is chosen for the new method. The unknowns, in this case, reduce to the maximum pressure values p 0 i at the centre of each strip i. At each iteration, the half-length of the contact strip a i in the rolling direction is then computed as a function of the maximum pressure distribution p 0 i , using an approximate analytical formulation based on Hertz' theory.

Discrete problem

The half-length and the half-width of each strip element are given as a i and b i ,1 respectively. The pressure distribution over each strip is assumed to be elliptical in the x direction, and constant in the y direction. Thus, the expression for the pressure distribution over each strip may be written as:

p n (x, y i ) = p 0 i 1 - x a i 2 , (2.7)
where p 0 i is the maximum pressure at the centre of the strip i. The deformation and the separation at the centre are denoted as u i and h i , respectively. Equation (2.1), Equation (2.3), and Equation (2.4) in discrete form can be written as:

∀i ∈ [1, n], u i = n j=1 A ij p 0 j , (2.8) ∀i ∈ [1, n] ∋ (x i , y i ) ∈ C, δ -u i -h i = 0 and p 0 i ≥ 0 , (2.9)
where n denotes the total number of elements in the potential contact area.

The terms A ij from Equation (2.8) are called the influence coefficients1 , and are defined as:

A ij = 2 πE ⋆ y j +b j y j -b j a j 0 1 -( ξ a j ) 2 (x i -ξ) 2 + (y i -η) 2
dξdη .

(2.10)

The expression for A ij describes the influence of normal stress distribution in the jth element to induce elastic displacement in the ith element. In the previous work by Reusner [START_REF] Reusner | Druckflächenbelastung und oberflächen-verschiebung im wälzkontakt von rotationskörpern[END_REF], the influence factors are expressed in the form of complete elliptic integrals. Knothe and Le The [START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with nonelliptical contact area[END_REF] determine these functions using numerical integration, however, no further details have been provided. In MIM-1D, the A ij terms are also evaluated numerically. The x i terms are ignored, as they always remain zero. Equation (2.10) is written as:

A ij = 2 πE ⋆   y j +b j y j -b j a j 0 1 -( ξ a j ) 2 -1 ξ 2 + (y i -η) 2 dξdη + y j +b j y j -b j a j 0 dξdη ξ 2 + (y i -η) 2 . (2.11)
The first integral is regular and can be evaluated numerically using Gaussian quadrature. The second integral is singular when the denominator approaches zero. This expression represents the case of uniform normal pressure acting on a rectangular area of 2a i × 2b i , and an analytical solution of this problem has been presented in detail by Love [START_REF] Johnson | Contact mechanics[END_REF][START_REF] Love | The stress produced in a semi-infinite solid by pressure on part of the boundary[END_REF]. We have: 2

y j +b j y j -b j a j 0 dξdη ξ 2 + (y i -η) 2 = (y + b j ) ln   a j + (y + b j ) 2 + a 2 j -a j + (y + b j ) 2 + a 2 j   + (y -b j ) ln   -a j + (y -b j ) 2 + a 2 j a j + (y -b j ) 2 + a 2 j   + 2a j ln   (y + b j ) + (y + b j ) 2 + a 2 j (y -b j ) + (y -b j ) 2 + a 2 j   , (2.12)
where y = y i -y j .

An additional advantage of decomposing the original integral expression in such a manner is the introduction of an analytical solution into the numerical results. This consequently reduces the integration error linked to the quadrature method. It also permits the use of lesser number of integration points, and consequently faster calculation times. Prior to this decomposition of the integral, an alternate strategy was to use a higher number of integration points closer to the diagonal terms A ii . Equation (2.8) in matrix form is written as:

u z = A p n , (2.13)
where the vectors p n and u z are given as {p 0 1 , ..., p 0n } T and {δ -h 1 , ..., δh n } T , respectively. The matrix of influence coefficients A is:

A =    A 11 . . . A 1n . . . . . . . . . A n1 . . . A nn    . (2.14)
This is the method of resolution generally followed when the rigid approach δ is known in advance. Unlike in the conventional approach, the matrix A needs to be evaluated at the beginning of each iteration, as the size of the elements in the x direction does not remain the same.

In the case where the normal force is prescribed instead of rigid body approach δ, the normal contact problem is solved with an additional iteration for δ. The resultant force at the end of each iteration is calculated using Equation (2.6) and Equation (2.7):

N = π n i=1 (a i b i )p 0 i .
(2.15)

The initial value of δ can be taken as the Hertzian rigid body approach. The subsequent values of δ may be evaluated using a dichotomy or an iterative scheme based on Hertz' relations [START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with nonelliptical contact area[END_REF]:

δ (m+1) = δ (m) N Ñ (m) 2 3 , (2.16)
where Ñ (m) represents the resultant normal force at the end of mth iteration.

Estimation of the contact patch half-length

The half-length of the contact patch a i is approximated using the equations that apply in Hertzian cases. a i is updated at each iteration to define the new potential contact zone, which is then used to construct the matrix A.

The equation for the Hertzian contact ellipse is:

x a 2 + y b 2 = 1 , (2.17)
where a and b are the semi-axes of the ellipse as defined in Section 1.4.1.3:

a = m 3 3 2 N 1 2E ⋆ 1 A + B , (2.18)
b = n 3 3 2 N 1 2E ⋆ 1 A + B . (2.19)
The normal pressure distribution p n (x, y) over the contact area is elliptical, and given by:

p n (x, y) = p 0 H 1 - x a 2 - y b 2 , (2.20)
with the maximum normal pressure at the centre of ellipse:

p 0 H = 3 2 N πab .
(2.21) Equation (2.17) over each strip can be written as:

a i a 2 + y b 2 = 1 . (2.22)
Combining Equation (2.20) and Equation (2.22), the normal pressure distribution in the longitudinal x direction for the ith strip is:

p n (x, y i ) = p 0 H a i a 1 - x a i 2 .
(2.23)

Comparing Equation (2.7) and Equation (2.23), it is possible to deduce:

p 0 i = p 0 H a i a = 3 2 N πab a i a .
(2.24)

Using the expressions for the semi-axes a and b defined previously, the contact patch half-length can be found using:

a i = π 2E ⋆ m 2 n A + B p 0 i . (2.25)
This Hertzian expression for the contact patch boundary depends on the curvatures of the profiles in contact, and the normal pressure at the centre of the strip under consideration. In non-elliptic cases, the curvatures and the Hertzian coefficients are replaced by their respective local values:

ãi = π 2E ⋆ m 2 i n i A i + B i p 0 i . (2.26)
The same expression for the half-length can also be obtained using the theory of ANALYN (see Appendix B). Hertz' solution remains valid only for positive values of B i . If the lateral curvature is negative at a given point, a correction must be carried out. Moreover, if the curvature is discontinuous, smoothing is also applied. This correction and smoothing are done using the procedure described in [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF]. The procedure involves two distinct steps described in Section 1.4.2.3:

1. First, the negative values of the lateral B curvature are limited to a minimal positive value.

2. Secondly, a smoothing filter based on the Boussinesq approach is used.

The assessment of ãi is a delicate matter in MIM-1D. Equation (2.26) is found to give good results in most cases when the curvature does not exhibit extreme variations along the profile. This sensitivity to the curvatures is a characteristic of all methods that use local values to estimate the half-length of the contact patch in the rolling direction, as shown in [START_REF] Sun | A robust non-Hertzian contact method for wheel-rail normal contact analysis[END_REF]. Several other solutions have also been tested, one of which will be presented later on in Chapter 4 when treating the specific case of sharp-edge contact encountered in switches and crossings.

To improve numerical stability, it is possible to modify Equation (2.13) using the normal force over each strip N i . From Equation (2.15), we have:

N i = πã i b i p 0 i .
(2.27)

The system of equations in matrix form in this case is:

u z = Ã N , (2.28)
where N = {N 1 , ..., N n } T . From Equation (2.13), we have:

A p n =    A 11 . . . A 1n . . . . . . . . . A n1 . . . A nn       p 0 1 . . . p 0n    =    A 11 . . . A 1n . . . . . . . . . A n1 . . . A nn       N 1 πã 1 b 1 . . . Nn πãnbn    = Ã N .
(2.29)

The modified matrix of influence coefficients à is now given by:

à =    A 11 πã 1 b 1 . . . A 1n πãnbn . . . . . . . . . A n1 πã 1 b 1 . . . Ann πãnbn    . (2.30)
From Equation (2.26) and Equation (2.27), the half-length of the contact patch at the end of each iteration can be computed using the normal force per strip:

ãi = 1 2E ⋆ n i m 2 i A i + B i 1 b i N i 1 2
.

(2.31)

Iterative resolution

Equation (2.26) is incorporated into the iterative algorithm presented in Figure 2.3 to solve the contact problem described by Equation (2.13).

In the presented test cases, the potential contact area is taken as the interpenetration zone using the Hertzian rigid body approach, as this is sufficiently large to enclose the actual contact area [START_REF] Paul | Contact pressures on closely conforming elastic bodies[END_REF]. The choice of this potential contact zone is not found to have a significant impact on the convergence of the proposed algorithm. When the problem is specified with only a given penetration, the interpenetration area is the optimal choice. To speed up the computation, the matrix of influence coefficients A is constructed only using the elements i where the separation h i is less than a predefined maximum value h max . When the resultant normal force N is known, the same algorithm is repeated for each value of δ (m) , evaluated using Equation (2.16). The total normal force at the mth iteration Ñ (m) is computed using Equation (2.15). The algorithm, in this case, converges when a user-defined tolerance value ϵ = | Ñ (m) -N | is attained.

Results and discussion

The proposed new approach, MIM-1D, is implemented as a Matlab function.

The approximate surface deformation method ANALYN [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF], and the VP method STRIPES [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF] are also programmed using Matlab. The results from the commercial version of the program CONTACT (v20.2) [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF] are used as the reference. It should be noted that the ANALYN results are sensitive to the method used for the negative-curvature correction. In the original publication, the negative lateral curvature values are replaced using a fifthdegree polynomial. However, this correction strategy has not been explained in further detail, and here this is done heuristically to obtain results as close as possible to the ones presented in [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF].

Theoretical profiles

Hertzian

To validate the new approach, a Hertzian case is considered with the contact between a sphere (R = 40 mm) and a flat surface. The separation curve between the two surfaces is shown in Figure 2.4a. Both bodies are made of steel, with ν = 0.3 and E = 208 GPa. The rigid approach between the two bodies is taken as 1 mm. The contact patches obtained using different approaches and the maximum pressure distribution p 0 (y) are shown in Figure 2.5. All the methods can be observed to be in good agreement with each other. The relative error in the contact area for MIM-1D is found to be within 1% of Hertz's analytical solution, which can be attributed to the accuracy of the numerical procedure used.

Non-Hertzian

The contact between a flat surface and a body of revolution (rolling radius R = 400 mm) with a non-Hertzian theoretical profile developed using two different radii R 1 = 40 mm and R 2 = 500 mm on either side of the point of the first contact is considered next. The separation curve between the two surfaces is shown in Figure 2.4b. The material properties are the same as in the Hertzian case, with the bodies pushed 1 mm towards each other.

The results for the contact area and the maximum pressure distribution are presented in Figure 2.6. From Figure 2.6a, MIM-1D and ANALYN correspond reasonably well with the reference results from CONTACT. The relative error in the contact area for MIM-1D is within 1% of the reference results. STRIPES notably underestimates the width of the contact patch: this is expected, as neglecting the surface deformation should lead to a smaller contact zone. Although the pressure distribution curve for ANALYN in Figure 2.6b follows the same trend as that of the reference, the peak of the pressure curve remains significantly higher. The results using MIM-1D can be observed to be in a better agreement with the reference method. A characteristic "bottleneck" region in the contact shape may be observed at around y = 0 mm. This can be attributed to Equation (2.26), which is used to estimate the contact patch length. Even though the stress distribution is relatively continuous over the entire contact zone, the ãi expression is essentially Hertzian. The terms m i and n i depend on the procedure used to compute the longitudinal and lateral curvatures A and B, respectively, the negative curvature correction as well as the smoothing strategy.

When the profile changes abruptly (e.g. R 1 = 40 mm to R 2 = 500 mm at y = 0 mm in Figure 2.4b, the discontinuity in curvatures might introduce a visible discontinuity in the contact area as well. These discontinuities may ultimately be treated in the pre-processing by using a more suitably adapted curvature smoothing process. However, these do not seem to be too significant usually in the case of wheel-rail contact, where the change in profiles is more gradual and not restricted to a sole point as in the considered theoretical cases.

Simulations are carried out using different potential contact zones to check the influence of this input parameter on the convergence of the proposed algorithm. The resultant normal forces are taken as 963 kN and 2645 kN (corresponding to δ = 1 mm) for the Hertzian and the non-Hertzian profiles, respectively. These results are presented in Figure 2.7. The outer loop iterations denote the number of iterations on the value of the applied rigid body approach δ (m) using Equation (2.16) to obtain the required resultant normal force. On the other hand, the inner loop iterations represent the average number of iterations required for the convergence of the algorithm presented in Figure 2.3 for each value of δ (m) . The choice of the potential contact zone visibly does not have a significant impact on the convergence.

Wheel-rail contact

The case of wheel-rail contact is presented using the standard wheel profile S1002 over the rail profile UIC60, with an inclination of 1:40. The profiles and their relative curvatures are presented in Figure 2.8. The material for both the wheel and the rail is steel. The resultant normal contact force is taken as 78500 N and the nominal rolling radius of the wheel R = 460 mm. The results are presented for various positions of the wheel, displaced from its centre position over the rail, denoted by ∆y. The sign convention is taken the same as in [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF], where a positive ∆y implies an outward movement of the wheel. The discretisation size is taken as approximately 0.2 mm for all the tested cases, which can be considered as a fine discretisation for wheel-rail contact applications.

The test cases are chosen to remain in the tread region with low contact angles, where the half-space assumption is not violated. From Figure 2.9, the MIM-1D results match closely with the reference CONTACT results in all the presented cases. This is highlighted more prominently in the maximum pressure distribution over the contact patch length. It can also be remarked that MIM-1D manages to accurately capture the characteristic slight variations in the pressure distribution, such as those presented in the tail end of the case ∆y = -1, a trait missing in the other simplified methods. ANALYN remains more precise compared to the VP methods, as neglecting the surface deformation in STRIPES again leads to an underestimated contact zone. The negative curvature compensation procedure used in STRIPES may have an effect as well [START_REF] Quost | Assessment of a semi-Hertzian method for determination of wheel-rail contact patch[END_REF]. This sensibility of STRIPES related to the processing of the curvature is also found to be true for ANALYN, with some fine-tuning required to obtain the desired results. A correction strategy dependent on the applied contact force and the separation may improve the results [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF]. comparison of the contact area for a range of ∆y values to the CONTACT results, emphasising the improvement using the new method as compared to existing fast approaches.

Computational cost

To make a representative comparison of the computational cost, MIM-1D is tested for different mesh sizes using a Hertzian profile. This comparison is done against the other presented methods i.e. STRIPES [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF], ANALYN [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF], and CONTACT [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF]. The computations are carried out using a 64 bits 2.70 GHz Intel processor. Only the time for the normal contact problem is considered. The CPU time for CONTACT is taken directly from the generated output file, while the other methods are measured using the elapsed CPU time-averaged over a finite number of runs. These results are shown in Table 2.1. STRIPES and ANALYN report similar and the fastest CPU times, which is expected as both methods are innately analytical. The latest version of CONTACT, implemented in Fortran, uses a bound-constrained conjugate gradient (BCCG) method with an FFT pre-conditioner which permits reasonably quick solutions even with a very fine discretisation [START_REF] Vollebregt | A new solver for the elastic normal contact problem using conjugate gradients, deflation, and an FFT-based preconditioner[END_REF]. The MIM-1D implementation here simply employs Matlab's inbuilt solver yet permits a gain in runtime compared to CONTACT, in no small part due to the reduced semi-analytical formulation of the method. Re-evaluating the matrix of influence coefficients at each iteration means that about two-thirds of the CPU time can be spent on the Gaussian quadrature. This may be an interesting guide for future developments, wherein accuracy can be traded for faster computation speeds.

Alternative pressure distributions in the rolling contact direction with MIM-1D

The main advantage of using strip elements to discretise the contact area lies in the reduction of the number of unknowns from n 2 to n (when considering the case of an equally spaced grid). Nevertheless, the requirement to numerically evaluate the influence functions of the elements at each iteration, as opposed to only one pre-computation in the case of conventional approaches, still leaves a little to be desired in terms of potential improvement of the semi-analytical method proposed here. In this section, we proceed to look at two possible alternatives to the elliptic pressure distribution assumed in the rolling x direction, for which closedform expressions of the displacement may be obtained.

Uniform pressure distribution

Considering the pressure to be distributed uniformly over the entire strip is perhaps the simplest alternative to the elliptic pressure distribution assumed in Equation (2.7). In this case, we have

p u n (x, y i ) = p u 0 i . (2.32)
The new influence coefficients A u,ij are defined as:

A u,ij = 2 πE ⋆ y j +b j y j -b j a j 0 1 (x i -ξ) 2 + (y i -η) 2 dξdη , (2.33)
for which the analytical solution is given by Equation (2.12). Using a uniform pressure distribution over the length of the strip, of course, implies that the maximum pressure at the centre of the strip is underestimated as compared to the expected value when considering an elliptical distribution.

In order to continue using the same iterative approach with the half-length of the contact area evaluated using Equation (2.31), we proceed to correct the normal force per strip, obtained through the new influence coefficients, by the following relation:

N i = N u i 0.9284 . (2.34)
This relation is approximated heuristically by comparing the results obtained using the standard procedure for several Hertzian cases with different aspect ratios of the contact ellipse.

Parabolic pressure distribution

Another interesting possibility is to consider that the normal pressure in the rolling x direction follows a parabolic distribution. Kalker in his monograph [21, p. 56] cites unpublished work where he implemented the same approach in an ad hoc manner for industrial applications. In [START_REF] De Mul | An orthogonal uniform/parabolic surface pressure half-space element for the analysis of non-Hertzian contact[END_REF], de Mul et al. use the parabolic pressure distribution for the analysis of non-Hertzian contact cases. This parabolic pressure distribution is integrated analytically over the surface element to obtain influence coefficients for the normal surface displacement. Their approach consists in developing two contact equations which are used to determine the contact pressures and the half-lengths of the contact zone simultaneously. We use the MIM-1D approach as presented previously, where the contact patch half-lengths are estimated using the local expression for a i given in Equation (2.26). The parabolic expression for the pressure distribution over each strip may be written as:

p p n (x, y i ) = p p 0 i 1 - x a i 2 , (2.35)
which replaces Equation (2.7). The new influence coefficients A p,ij are defined as:

A p,ij = 2 πE ⋆ y j +b j y j -b j a j 0 1 -( ξ a j ) 2 (x i -ξ) 2 + (y i -η) 2 dξdη , (2.36)
which can then be determined analytically [START_REF] De Mul | An orthogonal uniform/parabolic surface pressure half-space element for the analysis of non-Hertzian contact[END_REF] as:

2 y j +b j y j -b j a j 0 1 -( ξ a j ) 2 (x i -ξ) 2 + (y i -η) 2 dξdη = 1 a 2 j 4 3 a 3 j sinh -1 y + b j a j + 2a 2 j + (y + b) 2 3 (y + b) sinh -1 a j |y + b j | - a j 3 (y + b) a 2 j + (y + b) 2 - 1 a 2 j 4 3 a 3 j sinh -1 y -b j a j + 2a 2 j + (y -b) 2 3 (y -b) sinh -1 a j |y -b j | - a j 3 (y -b) a 2 j + (y -b) 2 .
(2.37)

Using the parabolic expression means that the maximum pressure p p 0 i at the centre of each strip exceeds the maximum pressure obtained using the elliptic expressions. Similar to the case with the uniform pressure distribution, the following relation is used to correct the normal force values at each iteration:

N i = N p i 1.054 , (2.38)
and hence the overall algorithm remains the same as the one presented previously.

To present the two alternative approaches, as previously, we consider the same wheel-rail contact cases for the different wheel positions presented in Figure 2.9. The comparison between the different MIM-1D algorithms are presented in Figure 2.11 and can be observed to be in relatively good agreement with each other. Although the real pressure distribution in the rolling x direction is close to a semi-elliptic function, the three chosen distribution functions give close results for the contact areas and the maximum pressures in the x = 0 plane. This is highlighted in Figure 2.12, where the contact area and maximum normal pressure values are compared with the reference results from CONTACT for a wide range of wheel positions ∆y. Nevertheless, it is important to point out that these alternatives should be looked upon purely as possible research avenues to be explored in the future, as they have not yet been tested and validated extensively. The zeroth-order approximations made heuristically in Equation (2.34) and Equation (2.38) means that the algorithm is not always as robust as in the case with semi-elliptic pressure distribution. These alternatives however still open up the possibility of even faster solutions .12: Comparison of the contact patch area (left), and the maximum normal pressure (right) for a range of wheel positions using different alternatives for the pressure in the rolling direction in MIM-1D when using a semi-analytical approach, with the view of implementation in MBS codes.

Conclusions

A simplified boundary element formulation is presented in this chapter and tested against two existing approximate methods, and a complete numerical method is used as the reference. In essence, this approach further develops the strip discretisation strategy [START_REF] Reusner | Druckflächenbelastung und oberflächen-verschiebung im wälzkontakt von rotationskörpern[END_REF][START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with nonelliptical contact area[END_REF] by using a semianalytical methodology to determine the contact patch dimension in the rolling direction. The novel method, implemented in the algorithm MIM-1D, provides a precise approach comparable to more rigorous complete methods such as CONTACT, with lesser computing effort. The results from the theoretical and wheel-rail test cases presented here allow the following conclusions to be drawn:

• The proposed approach MIM-1D enables an improved approximation of the pressure distribution and the contact area as compared to the other existing VP [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF] and approximate surface deformation [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF] methods;

• Comparison with a complete numerical method for wheel-rail contact shows close agreement, with the relative error in the contact area as compared to the reference results being less than 3% in all the tested cases using theoretical and wheel-rail profiles;

• From Figure 2.9, the normal stress distribution obtained using MIM-1D for the different test cases can be observed to be in better agreement with the fully detailed method. The peaks of the maximum pressure curve are also noted to conform well with the reference results. This may be seen as a significant advantage in the calculation of the tangential tractions and ultimately wear studies, which require the normal contact stresses instead of the total normal forces;

• The reduction in the number of system unknowns as compared to CONTACT due to a semi-analytical approach provides on an average, a 10-fold speed up with the current Matlab implementation. This should improve further by using a programming language closer to machine language such as Fortran;

• Two alternative possibilities are presented for the semi-elliptic pressure distribution assumed in the rolling x direction, considering it to instead be either uniform or parabolic. The availability of closed-form expressions of the displacement for these pressure distributions opens up the possibility of further optimisations when using a semi-analytical approach, with the view of implementation in MBS codes.

It is important to keep in context that the latest versions of CONTACT incorporate advanced numerical techniques to enable faster solving of the normal contact problem [START_REF] Vollebregt | A new solver for the elastic normal contact problem using conjugate gradients, deflation, and an FFT-based preconditioner[END_REF]. Using similar numerical optimisation strategies should permit further improvement in the performance of MIM-1D, and these will be tested in future developments.

Any simplified method of course brings its own set of drawbacks. Limitations concerning the generality of contact problems, especially when using rough profiles, are the most obvious. Considering the pressure distribution to be symmetric about the x = 0 plane also brings restrictions on taking into account the effect of the yaw angle. The bottleneck regions observed in certain cases, related to the abruptly changing curvatures, may demand supplemental inspections. Further developments must also be implemented to accurately treat the flange contact, where the contact angle varies significantly in a small zone, thus violating the half-space assumption.

In terms of the implementation of MIM-1D in the MBS code VOCO, the algorithm still needs refinement and robustness to compete with an analytical method such as STRIPES. The eventual implementation should account for the coupling with the tangential contact problem while taking the effect of friction into account. This requires specific tangential contact methods which enable the handling of non-Hertzian contact patches, an aspect which will be investigated further in Chapter 3. The new approach does not replace any of the previous ones but adds a new method in the spectrum of fast versus detailed methods. With proper optimisation, solving the normal contact problem using only strip elements should enable MIM-1D to be used as a good reference for other coarse models commonly employed in dynamic vehicle simulations.

Introduction

The quasi-identity assumption enables the separation of the normal and the tangential wheel-contact problems. The rolling contact problem in the context of MBS is then solved sequentially: the normal contact followed by the tangent one [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF]. The FASTSIM algorithm [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF] is widely used in MBS software packages for the evaluation of the tangential wheel-rail contact forces in a steady state. As the algorithm is restricted to Hertzian contact patches, a strip-based approach is proposed to extend FASTSIM to nonelliptical contact cases.

The following chapter presents this local approach based on curvature properties in detail [START_REF] Qazi | An extension of FASTSIM for steady state non-Hertzian contact[END_REF], which was first introduced by Ayasse & Chollet along with the semi-Hertzian method [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF]. Different settings for the traction bound are explored to determine their influence on the contact stresses, creep forces, and the limits of the saturation zone in the case of wheelrail contact. To validate this so-called FASTSIM SH approach, a design of experiments is constructed for different non-Hertzian contact cases, with different combinations of the longitudinal, lateral, and spin creepages. The absolute error in the normalised creep forces is used as the quantity of interest and compared with the reference results from CONTACT. The semi-analytical method presented in the previous chapter is combined with both FASTSIM SH and FaStrip algorithms to investigate their performance compared to the fully detailed method with a comparison of the tangential tractions, their directions and the stick-slip zone divisions. Finally, the influence of the normal contact modelling on the tangential results is also investigated by using two damage prediction models with the results from MIM-1D [START_REF] Qazi | A semi-analytical numerical method for modelling the normal wheel-rail contact[END_REF] and STRIPES [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF].

Extension of FASTSIM for steady state non-

Hertzian contact

The simplified theory

The FASTSIM algorithm based on Kalker's simplified theory [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF][START_REF] Vollebregt | FASTSIM2: A second-order accurate frictional rolling contact algorithm[END_REF] is perhaps the most widely used method in railway industry codes for evaluating the tangential contact parameters. The idea of the simplified theory is to replace the elastic body by a set of independent springs, as shown in Figure 3.1. Similar to a Winkler foundation, the tangential surface tractions p t and the tangential surface displacements u t at a given point are assumed to be linearly proportional through a flexibility parameter L. Thus:

u t = Lp t . (3.1)
In a steady state, the relative slip s t is defined as:

s t = w t - ∂u t ∂x , (3.2)
with x being the rolling direction in the local reference system. For concentrated contact problems, the creepages w t at a given point with coordinates [x, y] are defined using Equation (1.15):

w t = [ν x -yφ ν y + xφ] , (3.3)
where ν x , ν y , and φ are the longitudinal, lateral, and spin creepages respectively.

The slip s t is first assumed to be zero in Equation (3.2) and the adhesion region is supposed to cover the entire contact area. From Equation (3.2) and Equation (3.3), the longitudinal and lateral displacements u x and u y are then given as:

u x = (ν x -φy)(x -a i ) , (3.4
)

u y = ν y (x -a i ) + φ 2 (x 2 -a 2 i ) , (3.5)
where a i denotes the half-length of the contact patch at the ith ycoordinate. The contact patch considered in the FASTSIM theory is Hertzian. Integrating the expressions for the displacements over this elliptical contact area C with semi-axes a and b gives,

C u x dS = - 8 3 a 2 bν x , (3.6)
C u y dS = - 8 3 a 2 bν y - π 4 a 3 bφ . (3.7)
Kalker's linear theory establishes the following expressions for the creep forces in terms of the three creepages:

F x = C p x,ns dS = -Gabc 11 ν x , (3.8
) The expressions of the non-saturated shears p tns are given as:

F y = C p y,ns dS = -Gabc 22 ν y -G(ab)
p x,ns = - 3 8 Gc 11 ν x - 4 π Gc 23 b a yφ a i -x a , (3.11)
p y,ns = - 3 8 Gc 22 ν y a i -x a + 2 π Gc 23 √ abφ a 2 i -x 2 a 2 . (3.12)
These expressions correspond to the formulation with three flexibilities given in Equation (3.10), although an alternate formulation is also possible using only one flexibility parameter [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF]. For the total forces, the use of three flexibilities has been shown to offer better accuracy compared to using a single flexibility for a range of parameter values (creepages and aspect ratio) that occur for realistic vehicles [START_REF] Vollebregt | FASTSIM2: A second-order accurate frictional rolling contact algorithm[END_REF].

In the FASTSIM algorithm, the ellipse is discretised into strips with coordinate y i , of width δy and length 2a i . The ith strip itself is discretised into elements of length δx i as shown in Figure 3.2. For each ith strip, the tangential tractions p t are deduced iteratively, starting from the leading edge (j = 0) where the tractions vanish to zero. The expression for the tractions p ta , with adhesion being first assumed at element j, is given by:

p tij,a = p tij-1 -δx i ∂p tij,ns ∂x , (3.13)
where the partial derivatives of Equation (3.11) and Equation (3.12) are given as: The magnitude of the traction vector p tij at element [i, j] is limited by the traction bound, which is Coulomb's law applied locally:

∂p xns ∂x = 1 a 3
p tij =        p tij,a , if ||p tij,a || ≤ µp n , p tij,a ||p tij,a || µp n , if ||p tij,a || > µp n , (3.16)
where µ is the friction coefficient, and p n is the normal pressure.

From FASTSIM over an ellipse to FASTSIM over strips

To extend FASTSIM to non-elliptical patches, certain modifications must be introduced in the expressions presented in the previous section. First, the spin term φ associated with the longitudinal creepage ν x in Equation (3.3) vanishes, and the creepages at a given lateral coordinate y i become:

w ti = [ν xi ν yi + xφ i ] , (3.17)
where the subscript i indicates local values for each strip. If the spin creepage is supposed to be purely geometric, the neglected term from Equation (3.3) accounts for the rolling radius variation. This variation is instead considered in the local expression for the longitudinal creepage ν xi . When the contact angle γ i does not vary much over the patch, the contribution of the spin in the longitudinal component of w ti can be approximated as:

y i φ ≈ - y i sin γ i R = δr i R , (3.18)
where R is the nominal rolling radius and δr i is the rolling radius variation. Subsequently, the local expression for ν xi , as shown in Figure 3.3, is given as: The advantage of the above expression is that it enables taking into account any form of the rolling radius variation. Moreover, one does not need to define the origin of the tangent plane explicitly according to a predefined criterion (initial contact point, location of the maximum pressure or barycentre of the normal pressure distribution [START_REF] Vollebregt | Detailed wheel/rail geometry processing using the planar contact approach[END_REF]). The assumptions made in Equation (3.17) and Equation (3.18) still introduce a slight approximation of the linear theory in Equation (3.11) and Equation (3.12) for the nonsaturated shears, since the ν x term is proportional to c 11 and φ to c 23 in the original theory [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF]. The spin creepage is supposed to be purely geometric,

ν xi = ν x - δr i R . (3.19)
φ i = - sin γ i R , (3.20)
and the lateral creepage in each strip is similarly given as,

ν yi = ν y cos γ i . ( 3
p x,ns = - 3 8 Gc 11i ν xi a i a a i -x a i , (3.22) p y,ns = - 3 8 Gc 22i ν yi a i a a i -x a i - 2 π Gc 23i n i m i φ i a i a i a a 2 i -x 2 a 2 i , (3.23)
where m and n denote the Hertzian coefficients for the longitudinal and lateral semi-axes respectively. For non-elliptical patches, the ratio of the longitudinal semi-axes a i /a is replaced by a coefficient k i , between zero and one, whose value is a i /a in the case of an elliptical patch.

For convenience, the original FASTSIM algorithm uses normalised creepages as input, which are defined as the ratio of the non-saturated forces by the Coulomb's limit, both expressed over the ellipse [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF]. The output of the algorithm are the normalised creep forces f t :

f x = F x µN , (3.24)
f y = F y µN , (3.25)
where N denotes the normal force acting over the contact patch.

A similar procedure is used when extending FASTSIM to non-elliptical contact patches, the difference being that the normalised creepages are expressed for each individual strip. In order to derive the expressions of these entries of FASTSIM, the non-saturated forces F ti,ns are deduced from the summation of Equation (3.22) and Equation (3.23) over each strip. Thus,

F xi,ns = - 3 4 Gc 11i a i k i δy i ν xi , (3.26)
F y i,ns = - 3 4 Gc 22i a i k i δy i ν yi - 8 3π Gc 23i a 2 i k i n i m i δy i φ i . (3.27)
The normalised forces f ti,ns , and the associated normalised creepages A xi , A yi , and A φi are finally defined as,

f xi,ns = F xi,ns µ i N i = -A xi , (3.28)
f y i,ns = F y i,ns µ i N i = -A yi -A φi , (3.29)
where N i denotes the normal force acting over each strip.

The final adjustment to the original FASTSIM algorithm concerns the choice of the traction bound. Here, several possibilities may be explored: taking the traction bound either as parabolic or semi-ellipsoidal in the longitudinal x direction. According to Hertz theory, p n is semi-ellipsoidal. However, a parabolic traction bound has been shown to present better results in terms of the division of the stick-slip zones [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF][START_REF] Vollebregt | FASTSIM2: A second-order accurate frictional rolling contact algorithm[END_REF]. Through some mathematical manipulation, the expression of the normal pressure distribution p n,p corresponding to the choice of a parabolic traction bound in terms of normal force per strip is given as,

p n,p = 32 9π k i 3 4
N i a i δy i a 2 i -x 2 a 2 i . (3.30)
Similarly, the expression for the semi-ellipsoidal normal pressure distribution p n,e is given as,

p n,e = 2 π N i a i δy i a 2 i -x 2 a 2 i . (3.31)
Using a parabolic expression can sometimes lead to cases where the shears exceed the elliptic traction limit, which in turn implies that Coulomb's law is violated. In the proposed approach, the normalised output is weighted using the Hertzian expression to be coherent with Coulomb's theory. The weighing process consists in multiplying the shear stresses obtained using FASTSIM with a parabolic traction bound by the ratio p n,e /p n,p . This weighted parabolic traction bound ensures that the condition ||p t || µpn,e ≤ 1 is verified everywhere within the contact patch. The extended method is subsequently referred to as FASTSIM SH , where the subscript denotes the idea of the semi-Hertzian approach.

Results and discussion

In the following sections, the FASTSIM SH algorithm is tested in different contexts. In the first case, numerical experiments are carried out to determine the optimal settings for the traction bound. A design of experiments approach is subsequently used to validate the algorithm by comparing the obtained creep forces to the results obtained using CONTACT. Comparisons are also made with the FaStrip algorithm described in Chapter 1 when using MIM-1D for the normal contact results. Finally, the influence of normal contact modelling on damage prediction is studied briefly.

Numerical experiments to determine the traction bound

To demonstrate the effect of the different settings for the traction bound, we consider a non-Hertzian contact case using worn profiles of the wheelrail pair S1002-UIC60 with an inclination 1:40 [START_REF] Shackleton | Comparison of wheel-rail contact codes for railway vehicle simulation: an introduction to the Manchester contact benchmark and initial results[END_REF]. The nominal rolling radius is taken as 460 mm. The material properties are those of steel, with the modulus of rigidity G = 81890 MPa, and the Poisson's ratio ν = 0.27. The coefficient of friction µ is taken as 0.3. A multibody simulation of the passenger vehicle from the Manchester Benchmark [START_REF] Iwnicki | Manchester benchmarks for rail vehicle simulation[END_REF] running on a curved track is used to obtain the steady state input parameters. The normal contact force is 78.5 kN. The lateral position of the wheel over the rail ∆y = 5.8 mm, where a positive value indicates an outward movement of the wheel. The creepage values are taken as ν x = 0.58 ‰, ν y = 0.061 ‰, and φ = 0.274 m -1 , with the origin located at the point of geometrical contact.

The normal contact results for the considered wheel position and normal contact force are determined using CONTACT. The tangential stresses, their directions and the stick-slip zone division using the various available options are presented in Figure 3.4. The advantage of using a parabolic traction bound as opposed to a semi-ellipsoidal one is clear, which is also consistent with the results obtained previously by Kalker [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF] and Vollebregt [START_REF] Vollebregt | FASTSIM2: A second-order accurate frictional rolling contact algorithm[END_REF]. Nevertheless, the parabolic traction bound sometimes results in cases where the tractions exceed Coulomb's limit. The weighted parabolic traction bound ensures that the condition ||p t || µpn,e ≤ 1 is verified everywhere within the contact patch, as seen in Figure 3.4. The contact stresses and directions using the weighted parabolic setting can be observed to be in a relatively good agreement with the reference results obtained using CONTACT. The presented approach also provides an adequate estimation of the stick and slip zones, which is denoted using the solid line.

Validation of FASTSIM SH using CONTACT

In MBS codes, the creep forces rather than the contact stresses are used during online vehicle dynamics simulations. The original FASTSIM algorithm provides a good, fast estimation of the tangent forces, which is one of the reasons for its popularity. To evaluate the performance of the extended FASTSIM SH algorithm, we consider a combination of all three creepages within a range of realistic values. The longitudinal and lateral creepages are assumed to vary between 0 and 5 ‰, and the spin creepage from 0 to 5 m -1 .

In the statistical studies presented in [START_REF] Vollebregt | FASTSIM2: A second-order accurate frictional rolling contact algorithm[END_REF] and [START_REF] Vollebregt | Assessing the accuracy of different simplified frictional rolling contact algorithms[END_REF] for Hertzian contact patches, the ellipse ratio is considered as an additional varying parameter. Introducing a contact patch parameter in a non-Hertzian case is not so simple, where the normal contact results differ considerably depending on the method that is used [START_REF] Qazi | A semi-analytical numerical method for modelling the normal wheel-rail contact[END_REF]. We restrict ourselves to the five different wheel positions for the nominal S1002-UIC60 (1:40) profiles, as done in Chapter 2. Moreover, the normal contact results from CONTACT are used as the input for the FASTSIM SH algorithm. A design of experiments is constructed with 1000 simulations for each considered contact patch1 , using the FASTSIM SH algorithm with the weighted parabolic traction bound, and CONTACT. A uniform grid with a discretisation of approximately 0.58 mm is used for CONTACT. The same discretisation is used for FASTSIM SH in the lateral direction. In the longitudinal direction, two settings for the discretisation are used, with either 10 or 50 elements being considered for each strip. The absolute error in the normalised creep forces with respect to CONTACT,

|f FASTSIM SH x,y -f CONTACT x,y
|, as given by Equation (3.24), is used as the quantity of interest to be assessed, and is presented in Figure 3.5 as a function of the percentage of the total tested cases.

For the case of the normalised lateral creep forces f y , the maximum absolute error using 10 elements seems relatively high at around 0.18, although this drops significantly to about 0.06 for approximately 98% of the tested cases. Using m x = 50, the absolute error in f y is found to be less than 0.075 for all the tested cases, and similarly drops to below 0.03 for about 98% of the simulations.

In the case of the normalised longitudinal creep forces f x , the absolute error follows a similar trend, with comparatively better results. The error is found to be below 0.085 and 0.045 for all the tested cases using 10 and 50 elements respectively, which drops to below 0.02 for 90% of all the simulations that are carried out. In a previous study for a unique non-Hertzian case, different normal contact methods were used and the error in normalised creep forces was found to be less than approximately 0.1 for all the simulations using FASTSIM SH [START_REF] Qazi | An extension of FASTSIM for steady state non-Hertzian contact[END_REF]. The improvement in the results presented here is unsurprising, as the effect of the normal contact modelling is eliminated through the use of the reference program CONTACT for the tangential input parameters. In general, with a moderate discretisation being used in the longitudinal direction, the results of the FASTSIM SH seems to agree relatively well with the performance of the original FASTSIM algorithm for Hertzian contact cases [START_REF] Vollebregt | FASTSIM2: A second-order accurate frictional rolling contact algorithm[END_REF].

Combination of MIM-1D with a rolling contact model

In the following section, the semi-analytical method MIM-1D presented in Chapter 2 is combined with the FASTSIM SH as well as the FaStrip algorithms. The performance of both algorithms is compared to the fully detailed method through a comparison of the tangential tractions, their directions and the stick-slip zone divisions. Three contact models are considered to focus on the influence of the normal contact modelling on the tangential contact results. We have:

1. The CONTACT software [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF] based on Kalker's complete theory for the contact between two elastic half-spaces [START_REF] Kalker | Rolling contact phenomena: Linear elasticity[END_REF], which is used as the reference;

2. The semi-analytical normal contact method MIM-1D [START_REF] Qazi | A semi-analytical numerical method for modelling the normal wheel-rail contact[END_REF] with FASTSIM SH [START_REF] Qazi | An extension of FASTSIM for steady state non-Hertzian contact[END_REF];

3. The semi-analytical normal contact method MIM-1D [START_REF] Qazi | A semi-analytical numerical method for modelling the normal wheel-rail contact[END_REF] with FaStrip [START_REF] Sichani | An alternative to FASTSIM for tangential solution of the wheel-rail contact[END_REF][START_REF] Sichani | A fast wheel-rail contact model for application to damage analysis in vehicle dynamics simulation[END_REF].

Using the nominal S1002-UIC60 (1:40) profiles, two lateral positions of the wheel with respect to the rail are considered, with ∆y = 0 mm and ∆y = -1 mm respectively, which result in two different non-Hertzian contact cases. The normal contact results for different wheel positions using the various available normal contact methods have already been presented in Here, the normal pressure distribution over the contact area for the two aforementioned wheel positions using CONTACT and MIM-1D are presented in Figure 3.6 and can be noted to agree closely. The normal contact force is 78.5 kN. The creepage values are taken as ν x = 0.5 ‰, ν y = 0 ‰, and φ = 0.05 m -1 , with the origin located at the point of geometrical contact.

The tangential stresses, their directions and the stick-slip zone division for the two cases are presented in Figure 3.7. Both FASTSIM SH and FaStrip provide reasonably good estimations of the tangential contact results when compared to reference results from CONTACT. In general, FASTSIM SH tends to slightly underestimate the tangential tractions, while FaStrip tends to slightly overestimate them for both of the considered wheel positions. This trend has also been observed for other contact results not presented here, and is not surprising when compared to the results already available in the literature [START_REF] Sichani | A fast wheel-rail contact model for application to damage analysis in vehicle dynamics simulation[END_REF][START_REF] Qazi | An extension of FASTSIM for steady state non-Hertzian contact[END_REF]. The stick and slip zone division, represented 3.7 by the solid line in Figure 3.7, also follows the same trend, with a slight overestimation in the case of FaStrip and a slight underestimation with FASTSIM SH . The shear stress distribution is compared by plotting the longitudinal and lateral shear stresses p x and p y for the considered contact cases along the rolling direction at different lateral coordinates. These results are presented in Figure 3.8 and Figure 3.9. FaStrip tends to capture the nonlinear distribution of the shear stress in the stick area, which is one of the prominent advantages of using the strip theory. The shear stress distribution using FASTSIM SH is also non-linear due to the weighting process, however, it is captured less accurately as compared to FaStrip. Concerning the longitudinal and lateral creep forces, both contact models provide similar results, as shown in Table 3.1.

Comparison of rolling contact models used for damage prediction

In most cases, using the detailed integral relations implies that the normal stress distribution and the contact patch obtained using MIM-1D agree more closely with the fully detailed method based on Kalker's complete theory [START_REF] Kalker | Three-dimensional elastic bodies in rolling contact[END_REF] as compared to the simplified analytical approaches. This is an advantage when the interest of the user lies in the prediction of damage or fatigue analysis, which require the contact stresses instead of the total normal forces.

To highlight this advantage, two different approaches are used to evaluate the damage within the wheel-rail contact area. In the first model, an energy index presented in [START_REF] Dirks | The development of a crack propagation model for railway wheels and rails[END_REF] and denoted as EI, is used locally for each element. This permits the introduction of the local spin creepage:

EI = p x (x, y)[ν x -φy] + p y (x, y)[ν y + φx] . (3.32)
The second model used for the damage prediction is based on the idea that energy dissipated at the contact patch could also be the source of damage and crack initiation. This model is characterised by the calculation of the wear number, denoted as T γ , which is an estimator of the dissipated energy per metre of the distance travelled [START_REF] Burstow | A whole life rail model application and development for RSSB -continued development of an RCF damage parameter[END_REF]. The value of T γ helps in identifying regions where material removal through wear would be the dominant process, and those regions where the rolling contact fatigue (RCF) damage would be likely to accumulate, which is highlighted in Figure 3.10.

In the region of mild wear, it can be assumed that the energy would be likely to contribute to damage, whereas in the severe regime, cracks which were about to initiate would be removed by the process of severe wear. Furthermore, fatigue damage can also be extended to incorporate crack propagation prediction1 [START_REF] Dirks | The development of a crack propagation model for railway wheels and rails[END_REF]. The expression for T γ is given as:

T γ = F x ν x + F y ν y , (3.33)
where ν x and ν y are the longitudinal and lateral creepages respectively, while F x and F y are the longitudinal creep force and the lateral creep force respectively. While some authors additionally consider the product of the spin moment with the spin creepage in Equation (3.33), for example, in [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF] or [START_REF] Nielsen | Switch panel design based on simulation of accumulated rail damage in a railway turnout[END_REF], the influence of this term for railway dynamics applications is generally negligible compared to the previous two. Similar to the previous section, three contact models are used to present the tangential contact results. As the goal here is to investigate the influence of the normal contact method on the damage prediction, this time we include the semi-Hertzian method STRIPES. The two positions considered previously are such that the first one (∆y = 0 mm) corresponds to a case where the normal contact results using all three approaches are close to each other, and the second one (∆y = -1 mm) such that the analytical approach STRIPES leads to a comparatively large discrepancy with MIM-1D and the reference results from CONTACT. The normal pressure distribution, the tangential tractions, their directions and the stick-slip division when using STRIPES for the two wheel positions are shown in Figure 3.11, which are to be compared with those presented in Figure 3.6 and Figure 3.7. Here, FASTSIM SH is used as the tangential contact method, although as shown in the previous section, using FaStrip should lead to more or less similar results. The different contact models considered this time are:

1. The CONTACT software;

2. The semi-analytical normal contact method MIM-1D with FASTSIM SH ;

3. The analytical semi-Hertzian method STRIPES with FASTSIM SH .

The energy index distribution using the three contact models for the considered wheel positions are presented in Figure 3.12. At ∆y = 0 mm, the three contact models can be observed to be in relatively good agreement with each other. This is not surprising, as the normal contact results for the three approaches match closely. Nevertheless, the more evenly distributed normal pressure obtained using STRIPES means that the peak of the energy distribution at around y c = -11 mm is underestimated compared to the other two models.

The results at the wheel position ∆y = -1 mm are visibly different, which is linked to the dissimilarities in the normal contact results. MIM-1D and CONTACT are observed to be in good agreement with each other. The peak of the energy index using STRIPES is approximately of the same order, which corresponds to the normal pressure results presented in the previous chapter. However, the underestimation of the contact patch implies that a large zone undergoing the majority of the energy dissipation is entirely missed.

To present the results using the T γ model, steady state input parameters for the different contact models are obtained through a multi-body simulation of the passenger vehicle from the Manchester Benchmark [START_REF] Iwnicki | Manchester benchmarks for rail vehicle simulation[END_REF] running on the curved track shown in Figure 3.13. This track consists of a straight section, followed by a transition zone, and finally by a curve with a constant radius of 245 m. The curve transition is 25 m long with a linearly increasing curvature. The vehicle speed is taken as 100 km/h and As in the case of the EI model, the creep forces are evaluated using the FASTSIM SH tangential contact method, and compared to the reference results from CONTACT. The longitudinal and lateral creep forces as functions of the track position are presented in Figure 3.15. The creep forces can be observed to be in relatively good agreement with the reference results. This is expected as the FASTSIM SH algorithm has been shown to provide a good estimation of the tangential contact forces [START_REF] Qazi | An extension of FASTSIM for steady state non-Hertzian contact[END_REF]. In general, MIM-1D provides values that are relatively closer to the CONTACT results as compared to STRIPES linked to the better estimation of the normal contact results. Figure 3.16 shows the evolution of the wear number T γ as a function of the vehicle's position in the track. Variations begin to appear between the three approaches about halfway through the curve transition as the wheel moves outwards. This zone also corresponds to a high variation in the longitudinal, lateral and spin creepages. The steady state values a few metres after the end of the transition zone seem to be in better agreement with each other. While the T γ model provides a good fast estimate of the damage conditions, it is not easy to convincingly separate MIM-1D and STRIPES based on it alone, which explains the need of detailed approaches such as the EI model when a more thorough investigation is required.

Conclusions

Two non-elliptical tangential contact methods, namely the FASTSIM SH and the FaStrip algorithms, are coupled with the semi-analytical method MIM-1D [START_REF] Qazi | A semi-analytical numerical method for modelling the normal wheel-rail contact[END_REF] presented in Chapter 2. The FASTSIM algorithm based on Kalker's simplified theory [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF] has been adapted to be used in the case of steady state non-Hertzian contact by using the local geometric properties of the interacting bodies. It may be noted that other variants of FASTSIM for non-elliptical patches exist [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF], but they are rather based on mean curvature as opposed to local ones, and the derivation of creepages is carried out differently.

The FASTSIM SH approach is presented in detail and validated here.

The tangential tractions, their directions and the stick-slip zone divisions obtained with both the FASTSIM SH and the FaStrip approaches have been compared for different wheel-rail contact cases, using the results from the program CONTACT [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF] as the reference. The following global conclusions can be made subsequently:

• Numerical experiments are carried out to investigate the influence of different settings used for the traction bound in FASTSIM SH . A new variant is introduced using a weighing process which shows better results than the formulation proposed in the original algorithm in cases where the parabolic traction bound setting violates the local Coulomb's law. At the same time, the stick-slip zone division retains the merits of the parabolic setting used in the original algorithm;

• The FASTSIM SH has been validated through a design of experiments approach, considering different non-Hertzian cases with various combinations of the three creepages. The absolute error in the normalised creep forces is used as the quantity of interest. This error is found to be within an acceptable range, and corresponds well with the results already present in the literature for Hertzian contact [START_REF] Vollebregt | FASTSIM2: A second-order accurate frictional rolling contact algorithm[END_REF][START_REF] Vollebregt | Assessing the accuracy of different simplified frictional rolling contact algorithms[END_REF]. The absolute error in f y is found to be less than 0.075 for all the tested cases, and drops to below 0.03 for about 98% of the simulations. For f x , this error is found to be 0.045 for all the tested cases, which drops to below 0.02 for 90% of all the simulations that are carried out;

• The FASTSIM SH and FaStrip algorithms provide comparable creep forces results for two methods that are based on different theories. FASTSIM SH uses a direct extension of the FASTSIM method, by using parameters defined in every strip. The FaStrip algorithm uses a corrected strip theory [START_REF] Kalker | A strip theory for rolling with slip and spin[END_REF] for the tractions, and the FASTSIM algorithm to determine the stress directions in the slip region.

Although the creep force results are comparable, the non-linear pressure distribution in the strip theory corresponds better with the detailed method used in CONTACT as compared to FASTSIM SH ;

• As the basic principle remains more or less the same, FASTSIM SH is as fast as the original algorithm, and the same is also true for FaStrip. From the point of view of implementation in MBS codes, the FASTSIM SH method provides a relatively straightforward approach, with the global parameters used in the original FASTSIM approach replaced by their local values. FASTSIM is already widely used in the railway industry, and using the presented method should offer a good procedure for its extension to non-Hertzian contact patches;

• The influence of the normal contact modelling on the tangential results is investigated by using two damage prediction models. The normal contact results from MIM-1D and STRIPES are used with FASTSIM SH and compared with CONTACT. Using the wear number T γ does not permit a clear differentiation between the two methods, however, it provides a fast quantitative assessment of the damage conditions in vehicle dynamics simulations. For applications where a detailed investigation of the wear conditions in the contact patch is required, a more thorough indicator such as the energy distribution index may be used.

The new normal contact method MIM-1D coupled with either FASTSIM SH or FaStrip represents a new solution for the whole rolling contact problem, providing better results than the existing approximate methods usually implemented in MBS codes. This is accomplished with less computing effort than the reference method implemented in CONTACT. The ultimate proof is to certify these rolling contact models through their implementation in an MBS software, which would enable wide-ranging studies not limited to the theoretical cases generally investigated in the literature.

Parallel to the development of new contact methods, the existing one in the MBS code VOCO has been benchmarked in an international study dedicated to switches and crossings [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]. The next chapter describes the approach used in VOCO to model the wheel-rail contact as well as the railway dynamics. The results of the benchmark are also presented briefly. Following the benchmark, a study of wheel-rail contact methods in switches and crossings has been initiated with other universities. This study represents another benchmark, exploring in more detail the results provided by STRIPES and the Kik-Piotrowski methods described in Chapter 1. Some results are also compared with the new normal contact method MIM-1D.

Introduction

In continuation of the general overview on railway MBS presented in Chapter 1, this chapter begins by looking at the different steps of modelling and simulation of railways systems when using the MBS code VOCO. This is followed by the application of VOCO to handle two different cases studies in railway vehicle dynamics. The first case deals with the participation in an international benchmark of different commercial and academic MBS codes for simulations concerning switches and crossings (S&C), which require specific considerations. Following the benchmark, the assessment of wheelrail contact modelling has been particularly addressed in a collaborative study with other universities, with special emphasis on non-Hertzian contact and impact loads.

The MBS software VOCO

The multibody systems simulation code VOCO has been developed extensively at the Université Gustave Eiffel over the past three decades under different variants, specialising in wheel-rail contact and dynamics of railway vehicles. It is mainly dedicated to mechanical engineers working in the railway sector, as well as to railway operators, suppliers of rolling stock and track infrastructure, or design engineers and project managers of new systems, working on new or existing railway lines.

Through various ad hoc and proficient algorithms, VOCO is usually in the range of achieving real-time simulation rates, despite being mainly developed to implement advanced contact models in the past few years [START_REF] Chollet | The VOCO multi-body software in the context of real-time simulation[END_REF]. The possibility of carrying out more sophisticated simulations exists nevertheless, primarily for the multibody modelling of the track, mainly in response to the demands of the railway industry. The following sections outline the different steps of modelling and simulation using VOCO, defining the different functions that can be handled by this software, the user settings, and the main software output.

Software features

The software features of VOCO are governed by a great number of options. Several variants are present, as shown in Table 4.1, depending on the intended application and the desired level of complexity for the models to be simulated. For instance, several dynamic conditions may be considered, including constant acceleration, free acceleration, user imposed speed law, or braking/traction scenarios. Different types of vehicles may be considered when carrying out simulations, mainly railway, but also rubber-tyred tramway, rubber-tyred metro, or Guided Light Transit (French: Transport sur Voie Réservée or TVR).

Wheelsets represent special types of bodies as they are connected to the track through the rolling contact model. In addition to the common properties (mass, inertia, centre of gravity, external forces), a given wheelset may have rigidly connected wheels, independent wheels (where the rotations of the two wheels are independent of each other), or connected wheels with torsional elasticity (where the left and right wheel rotations are linked through a linear stiffness). A given wheelset can also be a motored wheelset (always with brakes), a carrier wheelset, a carrier wheelset without brakes, or a motored wheelset with motor on the left/right side (and brakes). One can assign a relative difference in nominal rolling radius between the left and right wheels, which can be important when taking a curve, passing through S&C, or when traversing a particular track defect. Considering flexible bodies (excluding the wheelsets) is possible via flexible body files generated from FE models via a modal reduction approach.

Bodies are linked together by connecting links and VOCO provides about ten types of connection elements, including unidirectional bi-slope springs with play, tridirectional springs and dampers, non-linear unidirectional spring with play, non-linear spring with asymmetric plays, air cushions, bislope dampers with stiffness in series and play, unidirectional or bidirectional elements with controlled dry friction etc.

Theoretical or measured track geometries may be considered. A theoretical track can be divided into a number of horizontal and vertical zones. For each horizontal zone, the zone length, the radius at the end of the zone, the cant at the end of the zone, and the gauge variation at the end of the zone may be defined. Similarly, for each vertical zone, the zone length, the zone start radius, and the zone end radius can be prescribed. Several predefined forms (dips, troughs, bumps etc.) as well as user defined functions can be introduced as track defects at a chosen longitudinal position in the track. In the standard use of VOCO, no dynamic track model is considered and the vibratory behaviour of the track is not taken into account within the time loop. When the track dynamics are considered, options are provided to take into account track dynamics without roll, track dynamics with roll, track dynamics with structure, and track dynamics via the co-simulation approach with an FE model of the track. Each of these models uses a cyclic representation of the track, except for the co-simulation.

The fairly complex modelling, as well as the reliable behaviour of the different contact models, is one of the main strengths of VOCO. The normal contact is handled either through the semi-Hertzian (STRIPES) method [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF] or the multi-Hertzian approach [START_REF] Ayasse | Paramètres caractéristiques du contact roue-rail[END_REF]. The STRIPES method, described in Section 1.4.2.3, is slower compared to the multi-Hertzian approach (in the order of around 5 times), but allows precise calculation of contact stresses with elliptical as well as non-elliptical contact patches. Options are provided to take the roll into account, or the separation of rails when the rail profile consists of multiple bodies (as in the case of S&C) and dynamic track model is considered. The multi-Hertzian approach is faster than real time, and generally sufficient if one is mainly interested in the overall results (accelerations, wheel-rail contact forces, forces in the connecting elements).

The contact patches in this case are, however, restricted to ellipses, although multiple contact can be handled. The tangential contact parameters may be computed using the CHOPAYA method, which is a functional approximation of the FASTSIM SH approach described in Chapter 3. Similar to the Shen-Hedrick-Elkins (SHE) method [START_REF] Shen | A comparison of alternative creep force models for rail vehicle dynamic analysis[END_REF], CHOPAYA evaluates the creep forces using empirical relations for given values of the local normalised creepages per strip. This functional approximation of FASTSIM SH may also be used in a Hertzian context, and it better approximates FASTSIM than the SHE method the when spin is significant. Options are also provided to use the SHE method, the Polach method [START_REF] Polach | A fast wheel-rail forces calculation computer code[END_REF], or the Polach with falling friction method [START_REF] Polach | Creep forces in simulations of traction vehicles running on adhesion limit[END_REF], although these are limited to the multi-Hertzian approach. The material behaviour within the contact patch may be elastic or perfectly plastic [START_REF] Sebès | A fast-simplified wheel-rail contact model consistent with perfect plastic materials[END_REF].

User settings

Vehicle and track modelling

The vehicle model in VOCO is represented as a set of bodies connected by links. The model classically consists of three levels of solids, with carbodies, intermediate solids (which include bogie frames, motors) and wheelsets. Each layer is assigned a number, as shown in Figure 4.1, using the indices (I, J, K) such that (I, 0, 0) always refers to a carbody, (I, J, 0) always represents an intermediate structure, and (I, J, K) refers to a wheelset. For each solid, the mass, inertia, and centre of gravity properties are defined by the user and stored in the vehicle model file. It is also possible to apply an external force on a solid. If the body flexibility is taken into account, the associated flexible body file is linked to the vehicle model. The Matlab interface allows the computation of the eigen frequencies, modal damping, and the visualisation of the eigen modes, which serves as a model verification tool as well as for the computation of the linear critical speed. This critical speed is defined as the speed above which the hunting motion becomes unstable. Saving the model results in the calculation of the static equilibrium of the vehicle, which is necessary as the subsequent dynamic calculations are made from an equilibrated vehicle, with initial static preloads [START_REF] Chollet | User manual: VOCO platform under MATLAB[END_REF].

The track geometry and the track defects (if applicable) are also specified at the input stage. If a dynamic track model is considered, the corresponding track model file is included in the input.

Wheel-rail contact

For the wheel-rail contact, both theoretical and measured wheel and rail profiles may be used. These profiles are smoothed to avoid numerical artefacts, and the degree of smoothing is defined by the user at the input stage. The contact parameters of STRIPES and CHOPAYA are pretabulated. The user defines the interval of lateral movements of the wheel relative to the rail ∆y, the number of strips used to discretise the rail profile, as well as whether a search for gaps or discontinuities must be performed, which is useful when considering S&C profiles.

The strips are not distributed at random: the zones in contact are determined under a large virtual load. Zones where there is no contact are coarsely meshed. Conversely, a finer mesh is used in zones with close and conformal contact. The option for a uniformly distributed mesh is also available. These settings are used for the creation of the contact tables. The For the creation of contact files, the rail cross-sections are represented by cubic splines. To ensure a consistent interpolation between two crosssections, as in the case of S&C, the profiles are automatically divided into bodies. For example, in the crossing shown in Figure 4.3 (not to scale), the wing rail and crossing nose are coloured red and blue, respectively. One spline is used per body. The strip discretization is indicated by the coloured zone in Figure 4.3. In the lateral direction, rail profiles are discretised using strips. The same number of strips is assigned to each body.

The wheel dynamics in VOCO is coupled to the rolling contact model through a combination of springs and dampers in the lateral and vertical directions, with spring stiffness K y,z and the damping coefficients C y,z such that:

Y = -K y (y w -y c ) -C y ( ẏw -ẏc ) , Q = -K z (z w -z c ) -C z ( żw -żc ) , (4.1)
where Y and Q represent the lateral and vertical wheel-rail forces, respectively, y w and z w denote the lateral and vertical positions of the wheel deduced from the dofs of the wheelset, respectively, while z c and y c are the vertical and lateral displacements of the wheel at the mean contact level,

y c = y r ± ∆y , (4.2)
where y r is the lateral position of the rail as prescribed by the track geometry and the defects as well as the track dynamics model if applicable. The ± symbol is related to the adopted convention of a positive ∆y corresponding to an outward movement of the wheel. Due to the addition of Equation (4.1), two auxiliary dofs per wheel are added to the dofs of the wheelset: the lateral position y c and the vertical interpenetration δtz seen in Section 1.4.2.3. z c is defined such that:

z c = z r + t z , (4.3)
where z r represents the vertical motion of the rail as prescribed by the track geometry and the defects as well as the track dynamics model if applicable. t z denotes the wheel vertical kinematics with respect to the rail, as shown in Figure 4.4. In the case of S&C, t z also varies as a function of the track position and plays an important role in the vertical dynamics.

For the crossing shown in Figure 4.3, due to its conicity, the wheel descends as it rolls along the wing rail, which is presented by the red curve The stiffness and damping parameters K y,z and C y,z are user defined, depending on the simulation conditions, for instance, if the wheel passes over a defect, or whether a dynamic track model is included in the simulation. The stiffness and damping coefficients can either be constant for the entire duration of the simulation, or variable as a function of the track position.

Several options may also be used for the friction model. The coefficient of friction can remain constant throughout the simulation, with values defined on the tread and flange regions globally, or locally in each strip. The friction values can also be variable as a function of the track position. The adhesion recovery model described in [START_REF] Toumi | Contribution to the certification of a WSP device in a HIL environment through a new adhesion recovery model[END_REF] is also a possible option, where the coefficient of friction is both space and time dependent. Using Polach's falling friction model [START_REF] Polach | Creep forces in simulations of traction vehicles running on adhesion limit[END_REF], it is also possible to define friction as a function of the creepage and the vehicle speed.

In addition to the initialisation parameters mentioned above, the global dynamic settings defined by the user at the initialisation stage include the vehicle speed, the acceleration (imposed or free), the duration of the simulation, the frequency of output storage, and the time step used for the solution of the equations of motion.

Software output

The software output may be divided into two parts, with global results for the vehicles, and the local results related to the wheel-rail contact, as shown in Table 4.2. The global output includes the wheel-rail contact forces for each wheel, the displacements, the velocities and the accelerations of each body in the MBS, as well as the resulting forces in the connecting links. With regard to the local contact results, options are provided to visualise the normal and tangential contact forces per strip for each wheel, the creepages for each wheel, as well as the contact patches and contact stresses at a given longitudinal position in the track.

The switches and crossings benchmark

The earliest railway networks used sliding rails to move vehicles between tracks until the railway switch was invented in the early 19th century. Switches and crossings (S&C), also known as turnouts, are now integral components of the railway infrastructure. S&C enable a railway vehicle travelling along a given rail to pass over or "switch" to the rail of a track which crosses its path, and thus ensure a fluid railway traffic operation.

A train can traverse the turnout in a straight line without changing its direction, called the through route, or run along the curve, which is called the diverging route. As opposed to the rest of the network, there is usually no cant to compensate curving acceleration in S&C. Moreover, there is no transition zone between the tangent track and the curve of constant radius. This implies that high lateral contact forces are usually observed in the diverging route.

A few important components of a turnout are shown in Figure 4.6. A switching machine or an actuator is used to operate the movable switch rail, which switches the wheels of the train from one rail to another according to the desired direction. The through route and diverging route rails cross each other at the crossing nose (or the crossing frog), which allows the wheels to travel along either of the intersecting paths. To switch the vehicle between different tracks, the running rails must be interrupted at the crossing level, leading to a void. This is to ensure that the wheel flanges of a train taking the other route have space to pass through to the other side. The discontinuities in the running rails may be avoided in the case of high-speed lines by using movable frogs. Check rails or guard rails are provided on either side of the crossing nose to enforce a constraint on the lateral position of the wheel.

The discontinuities resulting from the interruptions to the running rails lead to high impact loads as the wheel jumps onto the rail. Due to the rapid change of the rail geometry in S&C, high frequency responses (around 1 kHz) are obtained at the wheel-rail contact level. The dynamic wheel-rail interaction due to variation of the geometry and stiffness properties along the switch and the crossing panel, and the short transition zones also cause excessive wear, damage and rolling contact fatigue. The speed at which a train runs through the S&C needs to be carefully monitored, as high lateral loads can lead to a risk of derailment in a diverging route. It is easy to understand why turnouts are seen as one of the most maintenance-intensive components of the track. While a number of MBS benchmarks have been carried out by different universities and research groups, investigating a variety of modelling and application aspects such as the vehicle dynamic behaviour [START_REF] Iwnicki | Manchester benchmarks for rail vehicle simulation[END_REF], vehicle-track interaction at high frequencies [START_REF] Grassie | Models of railway track and vehicle/track interaction at high frequencies: Results of benchmark test[END_REF], wheel-rail contact modelling [START_REF] Shackleton | Comparison of wheel-rail contact codes for railway vehicle simulation: an introduction to the Manchester contact benchmark and initial results[END_REF], the pantograph-catenary interaction [START_REF] Bruni | The results of the pantograph-catenary interaction benchmark[END_REF], and longitudinal train dynamics simulations for freight trains [START_REF] Wu | International benchmarking of longitudinal train dynamics simulators: results[END_REF], a benchmark on S&C has been missing so far. It is partly out of this reason that a new benchmark was announced at the International Association for Vehicle System Dynamics (IAVSD) symposium, held in 2019, by researchers from University of Huddersfield (UK) and Chalmers University of Technology (Sweden) [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: modelling description and simulation tasks[END_REF]. The hope is that a deeper understanding of the vehicletrack interaction in S&C could lead to better design, better interpretation of the damage and wear processes, and consequently better cost management of maintenance operations for railway network operators. In order to be handled by a maximum number of participants, the benchmark setting do not correspond to a set of physical measurements. Therefore, the results cannot be compared to an absolute reference. Some comparisons with measured data are still available in the literature, including for VOCO [START_REF] Kassa | Simulation of dynamic interaction between train and railway turnout[END_REF][START_REF] Pålsson | Dynamic vehicle-track interaction in switches and crossings and the influence of rail pad stiffness -field measurements and validation of a simulation model[END_REF][START_REF] Sebès | Adaptation of the semi-Hertzian method to wheel/rail contact in turnouts[END_REF]. The S&C benchmark has received great interest, with 9 individual software participants (including VOCO), and a total of 18 universities/research groups/software developers/consultancy companies. The modelling descriptions and simulation tasks for the S&C benchmark, including the input data set, has been presented in detail in [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: modelling description and simulation tasks[END_REF].

Two different types of S&C configurations are considered: a Swedish 60E1-R760-1:15 and a British 56E1-R245-1:9.25. The denomination of a turnout refers, successively, to the nominal rail profile (60E1 or 56E1), the curve radius (760 m or 245 m) and the value of the tangent between the through and the diverging routes (1:15 or 1:9.25). The two configurations present two different types of curving and load transfer conditions. The passenger vehicle from the Manchester benchmarks is used as the vehicle model [START_REF] Iwnicki | Manchester benchmarks for rail vehicle simulation[END_REF], and the wheel profile is a nominal S1002 [START_REF] Shackleton | Comparison of wheel-rail contact codes for railway vehicle simulation: an introduction to the Manchester contact benchmark and initial results[END_REF]. Each participating institution was asked to submit a method statement for their respective software, detailing mainly the implementation of the variable track geometry, the interpolation procedure used for rail cross-sections, the implementation of the track model, the modelling of the changing curve radius when taking the diverging route, the wheel rail contact modelling and the time integration method used. These method statements are also available as a part of the provided data set, including the one for VOCO [START_REF] Sebès | VOCO statement of methods[END_REF]. A summary and comparative discussion on the modelling and simulation approaches used by the different participants can be found in the paper presenting the benchmark results [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF].

In the following sections, we outline the new developments integrated into VOCO as a direct outcome of the S&C benchmark participation: the development of a co-running track model, as well as some improvements in the modelling of the check rail and the contact stiffness. These developments have also been presented in the available method statement [START_REF] Sebès | VOCO statement of methods[END_REF]. Selected results from the benchmark, with comparisons of the different MBS software are also presented. 4.1, some track models are available in VOCO, but the one prescribed by the benchmark was missing. The S&C benchmark track model consists of a planar co-running system of masses and bushing elements, which is replicated independently for each wheelset [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]. This is a two level sprung mass system. Activating or disabling the different rail masses in the co-running track model for a given simulation enables the representation of the switch and crossing cases, as shown in Figure 4.7. The track mass m t has three degrees of freedom (vertical, lateral and roll), while all the rail masses have two degrees of freedom (vertical and lateral). In the switch panel, the rail 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 proposed co-running track model has been integrated as a new feature in VOCO. The model is composed of lumped masses labelled with double indices. The first index refers either to the left or the right rail, or to other bodies in the track model (which in this case is the track mass). The second index refers to a specific element of the selected given body. For instance, there are two elements for the right rail in the switch model: the stock rail and the switch rail. The bodies are numbered from the inside to the outside of the track, for example, the check rail is always the first element. This convention is adopted in order to match with the same convention used in the wheel-rail contact model. The wheel-rail contact model of the crossing panel is initially divided into two elements, representing the wing rail and the crossing nose. In the initial approach, the crossing mass in the track model was also divided into two masses linked by rigid springs. In order to prevent possible discrepancies with other participants, the contact model of the crossing was subsequently merged into a single element, allowing the use of the proposed track model.

New developments in VOCO

The track model is initially unloaded, and wheel loads are progressively applied on the track model during the first few meters of the simulation. At the entrance of the S&C panels, the track model is already in a steady state, with non-zero displacements. The elementary co-running track model allows very fast simulations, which in turn should ease the development of new contact methods. These new methods can then be validated through comparisons with more computationally expensive approaches, e.g. cosimulation with FE codes. The determination of the contact patches in VOCO is based on the vertical interpenetration of profiles through the semi-Hertzian method described in Section 1.4.2.3. An almost vertical contact angle, as in the case of a check rail, may therefore lead to numerical instabilities. Contact between the wheel and the check rail is taken into account by an equivalent spring and damper joint with a play in the normal direction, as shown in Figure 4.8 [START_REF] Sebès | Adaptation of the semi-Hertzian method to wheel/rail contact in turnouts[END_REF]. The play between the check rail and the wheel is a function of the location on the track and is determined from the cross-sections of the check rail. The joint stiffness is 100 kN/mm and its damping is 20 N.s/mm, and these values are chosen heuristically. While yaw is not an entry of the pre-computed contact table described in Section 4.2.2.2, a first-order approximation of the contact point offset due to the yaw, either backward or forward, is still taken into account. This offset is also considered in the evaluation of the longitudinal creepage ν x (see Section 1.2.4.2). The length of this offset ds with respect to the rail cross-section, as shown in Figure 4.9, is assumed to be:

ds = 2R(z CR + z f lange ) , (4.4)
where R is the nominal wheel radius, z CR is the distance between the top of the rail and the top of the check rail, and z f lange is the distance between the top of the rail and the bottom of the vertical flange of the wheel, as shown in Figure 4.8 and Figure 4.9.

In previous case studies, a total saturation within the contact patch was assumed, with the tangential force acting solely in the transverse direction, with no component in the longitudinal direction [START_REF] Sebès | Adaptation of the semi-Hertzian method to wheel/rail contact in turnouts[END_REF]. However, first

F x = -µN ν x ν 2 x + ν 2 y , (4.5)
F y = -µN ν y ν 2 x + ν 2 y , (4.6)
where F x and F y are the longitudinal and lateral tangential forces, respectively, ν x and ν y are the longitudinal and lateral creepages, respectively, µ is the friction coefficient, and N is the normal force acting between the check rail and the wheel. Taking the longitudinal forces into consideration in this manner seems to give good agreement with the results from the other software developers in the benchmark.

Modelling of the contact stiffness

In simulations without a track model, lateral and vertical stiffness and damping are introduced in VOCO at the contact level, and their respective values are user-defined in the input files. These parameters account for the combined track, contact and wheel stiffness, and their order of magnitude corresponds to typical values used in railway dynamics, with stiffness equal to 50 kN/mm, and the damping set to 200 N.s/mm [START_REF] Iwnicki | Manchester benchmarks for rail vehicle simulation[END_REF]. When a track In the initial tests, when using the values for the stiffness and damping 500 kN/mm and 20 N.s/mm, respectively, the contact stiffness was found to be too small to capture a short length phenomenon such as the combination of the P1/P2 peaks of the normal impact forces that are expected in a crossing [START_REF] Jenkins | The effect of track and vehicle parameters on wheel/rail vertical dynamic forces[END_REF]. Comparison between VOCO and an FE model in a rail joint simulation indicates a better value for the stiffness should be around 1000 kN/mm [START_REF] Nelain | Rail joint dynamic forces computation[END_REF]. Simulations with this order of magnitude present the expected peaks, but also lead to some contact loss.

For the benchmark, the standard model is improved by replacing the constant vertical stiffness K z with a Hertzian stiffness k h , as shown in Figure 4.10. The expressions used for the wheel dynamics in Equation (4.1) are subsequently modified as follows:

Y = -K y (y w -y c ) -C y ( ẏw -ẏc ) , Q = -k h (z w -z c ) 3 2 -C z ( żw -żc ) ,
k h = 2 3 r -2 3 E 1 -ν 2 1 √ A + B , (4.8)
where r is the Hertzian coefficient used to assess the rigid body approach, and A and B are the longitudinal and lateral curvatures, respectively. The lateral stiffness is set equal to 500 kN/mm, and a damping of 20 N.s/mm in both directions is considered sufficient to take the dissipative effects into account. These values of stiffness and damping are assessed heuristically based on previous internal studies. The damping values may still be assessed by assuming a complex contact stiffness such that:

Kz = K z (1 + jη) , (4.9)
with the damping ratio η of around 5%. As the complex stiffness is not applicable in the time domain, one must define an equivalent viscous damping C z , whose value is defined by:

2πf C z = K z η , (4.10)
where f is the characteristic frequency of around 1 kHz in order to capture the P1/P2 peaks. With K z around 1000 kN/mm, this leads to a C z value of approximately 10 N.s/mm. Figure 4.11 shows the comparison of the vertical Q forces obtained in a crossing panel when using different settings for the contact stiffness. When using the standard model with a small stiffness, a single peak is found at approximately 0.57 s. Increasing the contact stiffness value provides a better representation of the Q forces, with both the P1 and P2 peaks visible in approximately the same time interval. However, using these values, a wheel lift is also observed at around 0.6 s. Selecting the new model with the Hertzian contact stiffness enables capturing both the P1/P2 peaks expected in crossings, without the wheel lift observed in the previous case.

Simulation cases

The different simulation cases for the S&C benchmark are listed in Table 4.3, with 9 runs in total, including Run 9 as a baseline comparison of each participant's simulation set-up and vehicle model [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: modelling description and simulation tasks[END_REF]. The switch and the crossing panels are evaluated separately to maintain the integrity of their respective modelling and simulation approaches. The diverging crossing case, however, involves running along the switch curve before reaching the crossing. The varying S&C rail profiles are always located on the righthand side of the track, as shown in Figure 4.12. The markers in Figure 4.12 indicate the start of the S&C for the switch panel, and the intersection point (IP) for the crossing panel, respectively. Only a select few results, in particular for the wheel-rail contact parameters, are presented here. For a comprehensive comparison between the results from different software developers, one may refer to [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF].

Contact methods

Most MBS codes use one of the different VP methods or Hertz' theory for the normal contact, as described in Chapter 1. FASTSIM is the most popular choice for evaluating the tangential contact parameters, and is commonly combined with an equivalent ellipse when the contact patch is non-Hertzian. Table 4.4 lists the various normal and tangential contact models that have been used by the different software 1 developers for the benchmark simulation cases.

S&C benchmark results

Switch panel

If we consider the simulations in the switch panels in the through route, load is transferred from the stock rail to the switch rail, and also briefly shared between the two. Moving along the curved track in the case of the diverging track implies that this brief contact with both rails lasts longer as compared to the through route. The lateral wheelset dynamics leads to amplification of wheel-rail contact forces, which is comparatively more significant in the diverging direction. Figure 4.13 and Figure 4.14 show some results for the Run 1 in the through route with the British S&C. Here, the stock rail is shown on the left, between 0 m and 5 m approximately. The switch rail is shown between 1 m and 6 m approximately, corresponding to the top left configuration 1 Vampire was represented by VDG as a software user in Figure 4.12. Despite taking the through route, the diverging stock rail implies that the contact position follows the curved track, before the sudden jump onto the switch rail, as shown in Figure 4.13. The corresponding load transfer from the stock rail onto the switch rail also leads to a dynamic load amplification, as seen in the plot for the vertical wheel-rail Q forces in Figure 4.14 at approximately 3.4 m. In the majority of cases (including VOCO) the load transfer is not instantaneous but spans over approximately 20 cm. This is due to the flexible track model, if one takes into account the relative displacement of the switch and stock rail, as shown in Figure 4.15. An additional step is added to the contact determination, where the vertical position of each body is offset to take into account the vertical displacement of the flexible track model. In VOCO, the so-called separation of rails option needs to be activated to capture this phenomenon [START_REF] Sebès | Adaptation of the semi-Hertzian method to wheel/rail contact in turnouts[END_REF].

The contact patch location during Run 2, with the diverging route, is shown in Figure 4.16, corresponding to the top right configuration (mirrored) in Figure 4.12. Here, an earlier contact with the switch rail as compared to the through route can be observed. This is shared with the stock rail initially, before the two-point contact continues onto the switch rail once the wheel takes the curve. The contact patch size in Figure 4.17 shows the contact on both the stock and switch rails. An increasing contact size on the diverging stock rail can be observed, followed by a reduction of the patch size as it is shared with the contact on the switch rail. This is followed by a two-point contact on the switch rail, before a large contact patch increase as the load transfers back to a one-point contact. The variation between the different software highlight the differences that are mainly due to the various contact methods used by the different software developers. This is especially visible at the end of the switch panel, where a large variation in the contact patch size can be observed. Some codes predict a double contact, while others show only a single one.

A sharp-edge contact is observed in Run 3 and Run 4 of the S&C benchmark when the vehicle takes the through and the diverging route with the Swedish switch, respectively. As previously, the wheel initially follows the direction of the stock rail before jumping onto the switch rail, with a short period of double contact, as shown in Figure 4.18 and Figure 4.19. Unlike the British S&C, in this case the wheel makes a very brief sharp-edge contact at approximately 8 m as it moves to the switch rail. Even though the order of the contact duration is very small (in the order of milliseconds), the sharp-edge contact presents a situation which leads to damage of the switch blade due to extreme normal contact pressures exceeding the elastic limit of the rail material. The contact patch and the corresponding wheel-

Crossing panel

If we consider the simulations in the crossing panels, there is a very rapid transfer of load between the wing rail onto the crossing nose rail. Passing from the wing rail onto the tapered crossing nose leads to a large vertical dynamic impact-like load1 . Moving up the crossing nose as the wheel advances forward also gives rise to inertial forces, together with other lateral and steering effects due to the rapid changes in effective rolling radius [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF].

In the diverging route, the check rail running parallel to the stock rail on the opposite side of the crossing imposes a constraint on the lateral displacement of the wheelset. This ensures that the wheel does not clash against the crossing nose. The dynamics is thus further complicated due to the redistribution of forces between the opposite running rail, the check rail and crossing wing. the Run 7, considering the through route with the Swedish S&C, which corresponds to the bottom left configuration in Figure 4.12. The contact patch location in Figure 4.21 can be seen to move steadily along the path of the diverging wing rail before the contact transfers over on to the crossing nose. The vertical Q forces in the crossing panel, as shown in Figure 4.22, are characterised by a number of dynamic amplifications, first when the wheel moves along the wing rail, followed by a larger peak when the load transfers onto the crossing nose. A similar behaviour is also shown using VOCO in Figure 4.11. Brief contact losses due to the wheel rebound on the rail can also be observed for certain MBS codes. The P1/P2 peaks are found at a distance of approximately 0.5 m from the IP of the turnout. The contact patch location in the diverging route, as shown in Figure 4.23, follows a similar trend as in the through route in general, although the contact on the crossing occurs much nearer the rail gauge corner. The influence of the check rail on the wheel results in high lateral wheel-rail Y forces throughout the crossing panels, as shown in Figure 4.24. The force on the crossing rail is initially positive which is pulling the axle away from the crossing, before it changes directions as the wheelset is dragged towards the track centre line due to the action of the check rail [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]. The lateral forces reach their peak as we approach the point where load is transferred from the wing rail to the crossing. Compared to the other configurations, the crossing panel in the diverging route exhibits more dispersion between MBS codes due to the delicate modelling of the check rail.

Conclusions

The novel S&C benchmark provided an opportunity to incorporate several new features into the MBS code VOCO, developed at the Université Gustave Eiffel, including the new co-running track model, an improved modelling of the check rail, and improved modelling of the impact-contact process. The overall benchmark results from VOCO perform competitively against other academic and commercial MBS codes. Even the modelling choices where comparatively simpler approaches (e.g. pre-computed contact tables, vertical and lateral bushing elements approximating the contact stiffness) are used as compared to other software (e.g. online 3D optimisation for the contact detection, velocity dependant Hertz-based impact models with coefficient of restitution etc.) provide steady and reliable results [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]. Importantly, the wheel-rail contact forces do not present any anomalies, which further highlights the numerical stability of VOCO in realistic simulation conditions. Some conclusions from the benchmark can be highlighted as follows:

• In the switch panel with the through route, the transfer between the stock rail and the switch rail is relatively well modelled thanks to the separation of rails;

• In the crossing panel with the through route, the expected P1/P2 peaks are obtained using the improved modelling of contact stiffness, as well as a loss of contact especially in Run 5 (results not shown here);

• In the crossing panel with the diverging route, the modelling of the contact between the wheel back flange and the check rail is found to be consistent despite the simplifications.

Although a lot of results have been compared in the S&C benchmark, there was no room to compare the local results such as the contact shape or the normal pressure. This served as a motivation for some benchmark participants to further study specific configurations with a focus on the contact modelling procedures used in their respective MBS codes. This collaborative study is presented in the following section.

Modelling of wheel-rail contact in the presence of switches and crossings

As seen in the previous section, the passage of trains through zones containing switches and crossings may lead to various instances of sharpedge contact on the switch blade, or impact scenarios with the crossing nose. They may also lead to large conformal contact in the diverging route, especially for worn profiles. Conformal contact is not specific to S&C, but has received much attention in the railway dynamics community as its modelling represents a significant challenge using the existing fast methods. Kalker's variational theory [START_REF] Kalker | Rolling contact phenomena: Linear elasticity[END_REF] has been adapted for conformal contact in the latest implementations in the reference program CONTACT [START_REF] Vollebregt | Detailed wheel/rail geometry processing using the planar contact approach[END_REF][START_REF] Vollebregt | Detailed wheel/rail geometry processing with the conformal contact approach[END_REF], although the half-space approximation used in CONTACT is still violated for a sharp-edge contact.

The previous section shows a good agreement between the different software in the modelling and results for the kinematics and wheel-rail forces. Here, different aspects of the contact modelling in three MBS codes, namely, VOCO, MUBODyn, and VI-Rail, are investigated in detail. This is based on a collaborative study between few S&C benchmark participants, including Université Gustave Eiffel (France), Universidade de Lisboa, Universidade do Minho, and Instituto Politécnico de Lisboa (Portugal), and University of Huddersfield (UK) [START_REF] Magalhaes | Wheel-rail contact models in the presence of switches and crossings[END_REF]. While only three software have been chosen for this study, they cover a range of the available choices from the S&C benchmark [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF], with independent/commercial codes with an access to the internal details (MUBODyn/VOCO), as well as a commercially marketed MBS code, VI-Rail, where not as many details are provided. The goal is to study the handling of several contact scenarios that may occur in the negotiation of switches and crossings. Bearing that in mind, three different idealised running scenarios are used for the test cases, to achieve sharpedge or conformal contacts that may be seen in the switch, and the impact cases that take place at the crossing nose. The different aspects of the contact modelling that are considered for the comparison are now outlined for MUBODyn and VI-Rail, and compared to VOCO.

Aspects of contact modelling investigated for MBS software comparisons

Handling of variable rail profiles

MUBODyn and VI-Rail use a similar approach to model the rails with variable cross-section along the track via section breaks, as shown in Figure 4.25b. A section break is recommended when the profile width for two adjacent cross-sections differs by a certain tolerance. To build the section break, the profile with the larger width is trimmed such that the remaining part of the profile matches the profile with the smaller width. The trimmed profile is placed at the same longitudinal coordinate as the original smaller profile, and thus two cross-sections are defined at the same location. One geometric interpolation is done leading up to this location, and a continuing one from the same position as the two cross-sections match one another exactly on the overlapping segment [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]. A cubic interpolation in the lateral as well as the longitudinal direction is used to represent the 3D geometry of the rail.

In VOCO, a linear interpolation is used in the longitudinal direction to represent the 3D geometry of the rail. The rail cross-sections at each longitudinal position in the track are represented by cubic splines. To ensure a consistent interpolation between two consecutive cross-sections, profiles are divided into bodies, as described in Section 4.

Contact detection

MUBODyn uses an online 3D contact detection approach, taking in to account the effects of yaw and roll [START_REF] Marques | A three-dimensional approach for contact detection between realistic wheel and rail surfaces for improved railway dynamic analysis[END_REF]. The contact detection in VI-Rail is not well documented, but the yaw and roll are considered [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]. VOCO searches the contact in solely a 2D plane and the roll is taken into account.

Contact patch determination

In the MUBODyn version used in this study, the penetration function is determined based on the modified Kik-Piotrowski (MKP) model proposed in [START_REF] Sun | A robust non-Hertzian contact method for wheel-rail normal contact analysis[END_REF], which addresses some of the drawbacks of the original Kik-Piotrowski (KP) method. The length of each strip is then obtained according to the expressions given by the KP method [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF]. VI-Rail also uses relations given by the KP method. In VOCO, the semi-Hertzian or the STRIPES method, as described in Section 1.4.2.3, is used.

Evaluation of the creep forces

After the contact patch determination, MUBODyn determines an equivalent elliptical contact patch where the semi-axes a and b are evaluated to match the same contact patch area and same width of the original non-Hertzian contact patch, as described in [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF]. This equivalent elliptical contact patch is then used to evaluate the creep forces using Polach's creep force model [START_REF] Polach | A fast wheel-rail forces calculation computer code[END_REF], together with the global creepages defined at the point of maximum penetration, as shown in

Normal contact/Impact

In VOCO, the impact-contact process is modelled using the lateral and vertical stiffness and damping elements. Within the context of the S&C benchmark, the vertical contact stiffness is replaced by its Hertzian equivalent, as described in Section 4.3.1.3, which enables a better prediction of the P1/P2 peaks of the impact forces expected in a crossing [START_REF] Sebès | VOCO statement of methods[END_REF]. The same model is used in the current study. MUBODyn uses a Hertz-based model that allows some energy dissipation in the contact via a coefficient of restitution, which accounts for the damping and depends on the penetration depth, as follows:

N =      Kδ n c e if δ ≤ -v 0 , Kδ n [c e + (1 -c e )(3r 2 -2r 3 )] if -v 0 < δ < v 0 , Kδ n if δ ≥ v 0 , (4.11)
where r = ( δ + v 0 )/2v 0 is the transition parameter, K is the generalised contact stiffness that depends on the wheel-rail contact conditions, δ is the maximum virtual penetration of a contact patch, δ is its the first time derivative, n defines the degree of non-linearity (1.5 for wheel-rail contact according to Hertz' theory), c e is the restitution coefficient (typically 0.75), and v 0 is the penetration velocity tolerance (typically 0.1 m/s) [START_REF] Marques | A three-dimensional approach for contact detection between realistic wheel and rail surfaces for improved railway dynamic analysis[END_REF]. VI-Rail uses a parameter named as the "Hertzian damping ratio" [START_REF]VI-grade GmbH. VI-grade: VI-rail documentation[END_REF] which indicates that the normal contact force predicted in VI-Rail also depends on the speed of penetration, but unfortunately the implemented model has not been detailed in the software documentation.

Simulation cases

Three different scenarios are set up to explore the performance of the different MBS codes presented in the previous sections. These simulation cases are hereafter referred to as Run 1, Run 2 and Run 3, and are presented schematically in Figure 4.27. The first two runs present the wheel-rail contact conditions which may be observed when the vehicle negotiates a diverging route through a switch. In Run 1, the cross-sections of both wheel and rail are defined to achieve a large and conformal contact, as depicted in Figure 4.27a. This is similar to the type of contact that is obtained as the outer wheel is pushed against the rail when taking the diverging route. In Run 2, the rail cross-section is defined to represent an intermediate profile of the switch rail, which can have a very small radius of curvature, as shown in Figure 4.27b. This "sharp-edge" contact represents a contact condition for which the most common wheel-rail contact models, as described in Chapter 1, are not very realistic. Run 3 is used to investigate the impact at the crossing nose, when the wheel load is transferred from the wing rail to the crossing nose as schematically represented in Figure 4.27c with three consecutive conditions from top to bottom. For all the test cases, a simple bogie model, based on the Manchester benchmarks [START_REF] Iwnicki | Manchester benchmarks for rail vehicle simulation[END_REF], running at 100 km/h is considered. The track model is a co-running model, similar to the one used in the S&C Benchmark [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]. The parameters used for the bogie as well as the track models are given in Appendix D. The theoretical wheel and rail profiles used here are parametrized to be easily reproducible, but don't represent real profiles. The parameters used to describe the different profiles, and the crossing layout, are presented in detail in Appendix E. To ensure that the contact conditions represented in Figure 4.27 effectively occur, a curved track is considered for both Run 1 and Run 2. This curved track consists of a straight section of 5 m, followed by a curve transition of 25 m where the curvature changes linearly, and finally followed by a curve of constant radius 245 m. A tangent track is used in Run 3.

Results and discussion

Run 1: Conformal contact

Run 1 deals with the case of conformal contact as the vehicle takes the diverging route. The largest contact patches are expected at the outer rail where the flange contact occurs, and the wheel conforms to the rail profile, as shown in Figure 4.27a. The results for both the leading and trailing wheelsets are presented here, as opposed to the S&C benchmark, where only the leading wheelset was studied. The lateral coordinates of the contact positions on the outer (right) rail for the leading and trailing wheelsets are shown in Figure 4.29. A value near 0 mm indicates a contact near the tread, while a value near 40 mm indicates a flange contact, as can be seen from Figure 4.28. Figure 4.30 shows the area of the contact patches for the outer wheels of the leading and the trailing wheelsets.

For the leading wheelset, the wheelset dynamics predicted by the three software are more or less similar, as highlighted by their respective contact positions in Figure 4.29a, with one on the wheel tread and the other on the flange. One may therefore conclude that the principal differences in the contact areas are linked to the differences in the contact modelling, rather than the effects of the wheelset dynamics. A good agreement is observed between MUBODyn and VOCO for the smaller contact patch located on the flange (which is close to a Hertzian ellipse). VI-Rail exhibits a slightly larger contact area in this zone. This may be on account of the shape correction in the original KP method [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF], where the contact ellipse after the correction still has the same area as the interpenetration area, which is larger than the Hertzian ellipse [START_REF] Sichani | Comparison of non-elliptic contact models: Towards fast and accurate modelling of wheel-rail contact[END_REF]. For the larger contact patch on the tread, MUBODyn and VI-Rail show more or less similar results, while VOCO predicts a comparatively smaller contact area. With respect to the trailing wheelset in Figure 4.29b, the contact positions in VOCO show lower amplitude variation when compared to those obtained by MUBODyn and VI-Rail. The dynamic results also exhibit more variations between the three MBS codes than in the leading wheelset. The largest contact area is observed here on the wheel tread at the end of the transition curve at approximately 25 m, as shown in Figure 4.30b. In general, VOCO and VI-Rail exhibit significant differences, with higher contact areas obtained with VI-Rail. The results from MUBODyn show close agreement to VI-Rail, notably at the beginning of the curve transition, beyond which the wheelset dynamics have an influence on the contact position, and consequently on the contact areas obtained using these two approaches. This is again not surprising, as both MUBODyn and VI-Rail use adaptations of the KP method [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF][START_REF] Sun | A robust non-Hertzian contact method for wheel-rail normal contact analysis[END_REF], while the semi-Hertzian method STRIPES is used in VOCO.

The dimensions of the contact patch are compared for the outer trailing wheel when it reaches the beginning and end of the curve transition, i.e., at track distance equal to 0 m and 25 m, as shown in Figure 4.31. Note that VI-Rail is not considered here. The available contact patches are compared with the reference results from CONTACT [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF]. The contact patches using the MIM-1D approach described in Chapter 2 are also presented simultaneously. The main input parameters required for CONTACT are the undeformed distance (or the separation function), and the normal contact force. The difference in the lateral position of the wheel with respect to the rail, in particular during and after the curve transition, implies that the separation obtained for the two software are different, and hence the two contact cases are presented separately.

At 0 m, an elliptical contact patch is obtained using both software, as well as with MIM-1D, and can be seen to be in good agreement with CONTACT as well as with each other, as shown in Figure 4.31a and Figure 4.31c. In this case, the contact point is observed at the origin of the wheel, where the contact angle is approximately zero and the curvatures are locally constant. Wider contact patches are observed at the end of the curve transition, as shown in Figure 4.31b and Figure 4.31d, where the contact patches deviate from the elliptical shape. The contact patch obtained with MUBODyn still resembles a Hertzian ellipse, with a slight variation of the curvature on one side, and this shows a good agreement with CONTACT. On the other hand, the patch at 25 m for the contact position predicted by VOCO is strongly non-Hertzian, and slightly underestimates the patch area when compared with CONTACT. MIM-1D predicts a contact patch that is much closer to the reference contact patch obtained using CONTACT. The normal pressure distribution in the contact patches using VOCO, CONTACT and MIM-1D are presented in Figure 4.32, and seem to be consistent with similar non-Hertzian cases presented in the literature [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF][START_REF] Qazi | A semi-analytical numerical method for modelling the normal wheel-rail contact[END_REF]. With the rail and wheel profiles parametrized to provide large conformal configuration, a slight difference in the lateral displacement of the wheel induces a comparatively high difference in the contact position and area obtained using the two MBS codes used here.

The contact forces in the outer wheels of the leading and trailing wheelsets are presented in Figure 4.33, including the normal contact force, and the longitudinal and lateral creep forces. In general, MUBODyn and VI-Rail show good agreement with each other for the normal contact forces, in particular because they both use different variants of the KP method. For the leading wheelset, the normal contact forces are characterised by a sharp variation at the start of the curve transition as the contact area changes suddenly and a second contact patch appears. The longitudinal creep forces in the two contact patches tend to oppose each other due to the variation of the rolling radius between both locations. Similar to the normal contact forces, sharp peaks can also be observed for the lateral creep forces at the beginning of the curve transition as the wheel flange makes a second contact with the rail.

For the trailing wheelset, the normal contact forces present several peaks corresponding to the swaying motion of the contact position seen in Figure 4.29b. While a generally good agreement is observed between all three MBS codes for the creep forces, MUBODyn does provide lower values, especially in the case of the leading wheelset. These discrepancies may be linked to the assumptions made regarding an equivalent ellipse in MUBODyn when using Polach's method, whereas VOCO and VI-Rail both use different adaptations of FASTSIM for non-elliptical contact patches.

Run 2: Sharp-edge contact

Run 2 focuses on the sharp-edge contact that is typically observed when the wheel makes contact with the switch blade. The contact positions on the rail for the outer wheels of the leading and trailing wheelsets are shown in wheelset dynamics for the leading wheelset as evident from contact positions shown in Figure 4.34a. For the trailing wheelset, MUBODyn and VI-Rail predict similar contact positions during most of the simulation, although the variations increase as the vehicle exits the curve transition zone and the wheelset dynamics takes over. The initial contact position in VOCO differs from those obtained with MUBODyn and VI-Rail and these differences also persist for the remainder of the simulation. Similar to the Run 1, contact patches obtained with VOCO, MIM-1D, and MUBODyn, at the beginning and the end of the curve transition, are shown in Figure 4.35 and compared with the results from CONTACT. At the beginning of the curve transition, a half-elliptic shape is obtained with both software, which corresponds to the contact case with the sharp-edged rail profile. At the end of the curve transition, the size of the contact patches increases as the outer wheel pushes against the rail profile. The contact patches obtained using MUBODyn and CONTACT, as shown in Figure 4.35b, correspond closely, while the contact patch provided by VOCO shows a large difference when compared to the CONTACT results.

While the variations in the contact positions obtained using MUBODyn and VOCO, and the corresponding differences in the undeformed distances and normal contact forces, should have some influence on the contact patches, the similarities in the CONTACT results obtained via independent input from VOCO and MUBODyn suggests this may not be the sole reason for the discrepancies. Further investigation reveals that this is rather linked to the smoothing of the B curvature, which is inherent to the STRIPES method [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF]. STRIPES, as implemented in VOCO, involves the evaluation of a contact stiffness per band k i , which is given as:

k i = E 2(1 -ν 2 ) 1 + A i /B i n 3 i δy i , (4.12)
where n i are the local Hertz' coefficients for the lateral semi-axes, and δy i is the width of each band. Using the smoothed curvature shown in Figure 4.36a leads to higher contact stiffness, and consequently to a smaller contact patch.

Using the non-smoothed curvature leads to better results in terms of contact surface, as shown in Figure 4.36b, but the abrupt change of curvature at y = 0 leads to a sharp variation in the contact patch. Using the nonsmoothed curvatures does not have a significant effect on the contact forces. An alternative procedure can be to trim the profile to remove the sharpedge, which may lead to more consistent contact patch results with the other methods. Although the KP method used in MUBODyn is similar to STRIPES in its strip-based approach, it circumvents this issue as the mean curvatures rather than their local values are used for the evaluation of the contact patch dimensions.

A similar issue is also encountered with MIM-1D when using Equation (2.26) or Equation (2.31) to determine the half-length of the contact patch within the iterative loop, which is appropriate and robust for bodies whose profiles are both smooth and continuous. Johnson remarks in [18, p. 150] that the presence of a sharp corner at the edge of the contact, when using such a method, would lead to the pressure rising towards infinity. With the expression of the half-length used in MIM-1D dependent on both the pressure distribution and the local curvature properties, initial results using the standard MIM-1D approach described in Chapter 2 led to unsatisfactory results when treating the sharp-edge contact presented here.

The expression for the half-length a(y) was subsequently modified by using mean curvatures, inspired by the KP method, with Equation (1.63) given by the semi-Hertzian theory. This was found to give better results, as shown in Figure 4.35a and Figure 4.35b. The modified expression of a(y) used with MIM-1D in this case is presented in detail in Appendix C. This is a workaround which works well in the contact case considered here, but provides less accurate results when considering the wheel-rail contact cases treated in Chapter 2. It should also be noted that the sharp-edge contact presents a situation where any method using the half-space assumption, including CONTACT, is not valid any more. Figure 4.37 shows the normal pressure distribution in the contact patches of the outer trailing wheel obtained with VOCO, MIM-1D, and CONTACT at 0 m and 25 m. The extreme pressure values near the sharp-edge (at approximately y = 0) exceed the elastic limit of steel, and plastification is likely to occur after the first few passages. Any result using such methods must therefore be interpreted with extreme care.

The contact forces acting on the outer trailing wheel are plotted in Figure 4.38. Similar to Run 1, the normal contact forces are relatively in good agreement with each other, while the creep forces show differences owing to the different tangential contact method used as well as the different contact patches obtained with each software. MUBODyn seems to show some numerical instabilities, such as the sharp peak at approximately 16 m for the trailing wheelset, which are mainly caused due to the difficulties in the online determination of the contact points near the sharp-edge.

Run 3: Impact

Run 3 focuses on the transfer of wheel load to the crossing nose, where an impact-like situation is observed as a result of the geometry of the crossing panel. The contact location along the crossing panel, as well as the normal contact forces at specific intervals, are shown in Figure 4.39. As also seen in the S&C benchmark, the contact position in the crossing panel initially deviates laterally from its position, moving along the diverging wing rail, before it jumps abruptly onto the crossing nose. This is followed by a brief loss of contact as the wheel rebounds on the crossing before settling down.

The start of the crossing panel is located at 0 m. The first contact with the crossing nose occurs at around 0.92 m, as seen in the plot for interval 1 (from 0.9 m to 1.1 m). During this first impact, all three software are in phase and agree closely. The plot for interval 2 (from 1.5 m to 2 m) shows the second impact with the nose, after the loss of contact at the end of interval 1 when the wheel rebounds on the rail. Here, VI-Rail predicts the second impact slightly sooner than both VOCO and MUBODyn. While the contact search algorithm for VI-Rail is not well documented, this may suggest that VI-Rail searches potential contact points along the entire surface of the wheel, and not just in the 2D plane as done in VOCO, although MUBODyn also uses a three-dimensional approach for the contact detection [START_REF] Marques | A three-dimensional approach for contact detection between realistic wheel and rail surfaces for improved railway dynamic analysis[END_REF] suggesting that this might not be the reason for the discrepancies. Another possible explanation could be that the vertical velocity during the rebound after the first impact is lower in VI-Rail compared to the other software. This is possibly because of the damping used in the normal contact model which causes the wheel to descend quicker and make an earlier second contact with the crossing. The effect of the vehicle modelling in the MBS codes, as well as the influence of the tangential forces is likely to be negligible in Run 3, which is on a tangent track unlike during Run 1 and Run 2, and the normal contact algorithm seems to be the key factor for explaining the differences.

Conclusions

The collaborative study presented here allows investigation of the different aspects of contact modelling used in three MBS codes, namely, VOCO, MUBODyn, and VI-Rail, in detail, focusing on local contact results. With three different test cases simulating different configurations that may be seen in S&C, the study represents another benchmark that can serve as a guide for current and future software developers. The following conclusions can be made:

• In the case of conformal contact, fast contact methods seem to provide sufficiently precise results as compared to the reference results from CONTACT;

• The sharp-edge contact presents a contact configuration where any approach using the half-space approximation is no longer valid. Furthermore, the sharp discontinuity implies that normal contact methods using local curvature properties may provide results that are unsatisfactory. However, without a reliable reference, it is not straightforward to determine the accuracy of contact results in such a case;

• For the impact contact cases observed in the crossing nose, the normal contact modelling plays an important role not only in the representation of the expected P1/P2 peaks, but also on the instant when the second impact is detected by the contact search algorithm.

Searching for the contact only in the vertical plane may also be insufficient in cases where the running rails are interrupted, due to the presence of a discontinuity, such as in a crossing.

This latter study also permits the comparison of the computational efficiency of three MBS software, which wasn't done in the S&C benchmark. In general, the compiling language, optimisation, and the pre-tabulation of the contact parameters make VOCO the fastest code, achieving close to real-time simulations, which is one of the main objectives of the developers. The average CPU time per second of simulation using VOCO is about 15 seconds, 30 seconds for VI-Rail, and 1 hour for MUBODyn. MUBODyn took a few hours per case as the S&C features in MUBODyn are new, and only implemented in a non-compiled language, as opposed to VOCO or VI-Rail.

In general, wheel-rail contact modelling remains a challenging topic that requires further research to improve the capabilities and precision of existing MBS codes. Global results in terms of wheel-rail contact forces obtained in this study are found to be generally in good agreement, as was already case in the S&C benchmark [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]. Still, if one looks closely at the local results such as the contact patches and normal stresses, more variations may emerge which justifies the development of more accurate methods such as MIM-1D or FASTSIM SH .

Conclusions & future work

This PhD work aims to address the need for efficient and accurate modelling of the contact phenomenon for application to railway dynamics simulations. A thorough investigation is carried out on the various rolling contact models in the literature which are available for the prediction of wheel-rail contact results, ranging from the analytical formulae presented by Hertz to Kalker's detailed variational theory. The choice of wheel-rail contact model used in MBS software is foremost governed by the intended application of the study. For the normal contact approximate analytical models based on Hertz' theory, may often be sufficient when numerous simulations must be carried out per second for a train running over several kilometres of track. The real geometry of the wheel and the rail imply contact conditions that are strongly non-Hertzian, and virtual penetration methods are commonly used in MBS codes to determine the wheel-rail contact in these cases. These methods enable better modelling of the contact conditions compared to the Hertzian theory while avoiding the precise but computationally expensive fully detailed methods.

A similar issue is also encountered for the modelling of the tangential contact, with many of the available methods restricted to Hertzian ellipses. This implies that even MBS packages that use non-Hertzian modelling for the normal contact are forced to revert to strategies such as equivalent ellipses for the evaluation of the tangential creep forces. While these approaches work sufficiently well in dynamic simulations, they cannot be used to study the contact stresses or the slip-adhesion characteristics within the contact patch, which are required for investigations of the rail wear or damage. Thus, the research of more accurate solutions for both the normal and the tangential contact problems is still an open field in the context of railway vehicle dynamics.

In the first part, this work proposes the new method MIM-1D to further enhance the precision of the normal contact solution. This is achieved through a semi-analytical methodology, using strip elements instead of a 2D grid to discretise the potential contact area together with the expressions from Boussinesq's theory to evaluate the normal contact pressures accurately. The length of the contact patch is supposed to be a quasi-known quantity, dependent on the form of the normal stress distribution and the profile curvatures, and is evaluated using an expression based on the Hertzian theory applied locally. This ad hoc semi-analytical approach offers better estimations of the contact patch and pressure distribution when compared to the existing virtual penetration methods, using the software CONTACT as the reference. While it is slower than the virtual penetration methods, it is shown to be more computationally efficient than the rigorous solution implemented in CONTACT, thus adding a new method to the spectrum of fast versus detailed methods available in the literature.

The FASTSIM algorithm is currently the most widely used approach in MBS codes to assess the tangential creep forces. However, being restricted to elliptical Hertzian cases means that workarounds such as equivalent ellipses must often be used for the contact patch. The current work addresses this shortcoming by presenting an updated version of the FASTSIM algorithm, the so-called FASTSIM SH algorithm which was first introduced by Ayasse & Chollet. This approach extends the original algorithm to non-Hertzian contact patches by considering the local curvature properties, with the creepages being assessed in each strip of the contact patch. Numerical investigations are carried out to investigate the influence of different possible settings for the traction bound in FASTSIM SH on the slip-stick zone division and the tangential stresses. A new variant is introduced using a weighing process which shows better results than the formulation proposed in the original algorithm in cases where the parabolic traction bound setting violates the local Coulomb's law. The FASTSIM SH algorithm is also validated for different non-Hertzian cases using a statistical design of experiments. In the tested cases, this approach is shown to retain the merits of the original algorithm when evaluating the creep forces, and at the same time provides reasonably good approximations of the tangential contact stresses. The normal contact method MIM-1D is also coupled with the FASTSIM SH algorithm to investigate the influence of the normal contact modelling on the tangential results, particularly for the case of damage prediction. More accurate modelling of the normal contact via MIM-1D combined with a fast non-Hertzian tangential contact method provides a new solution for the whole rolling contact problem, showing better results than the existing approximate methods usually implemented in MBS codes while using less computing effort compared to the reference method implemented in CONTACT.

Alongside the development of new wheel-rail contact methods, this PhD work also addresses case studies in railway dynamics within the context of an international benchmark on switches and crossings, with a specific focus on wheel-rail contact. Several new features have been integrated into the MBS code VOCO including a new track model and improved modelling of the check rail and the contact stiffness. Even the modelling choices where comparatively simpler approaches (e.g. pre-computed contact tables, vertical and lateral bushing elements approximating the contact stiffness) are used as compared to other software (e.g. online 3D optimisation for the contact detection, velocity dependant Hertz-based impact models with a coefficient of restitution etc.) are shown to provide steady and reliable results, thus justifying the choice of the developers to focus on real-time simulation capabilities of VOCO. Following the benchmark, an indepth collaborative study has been carried out to investigate the influence of different aspects of contact modelling on the dynamic behaviour of the wheelset. This study presents the results using three MBS software, including VOCO, with parametrized wheel and rail profiles to reproduce contact conditions often encountered in S&C simulations. Between the different MBS software, global results in terms of wheel-rail contact forces are generally in good agreement, as is also highlighted in the S&C benchmark. The local results such as the contact patches and normal stresses may exhibit more variations, which in turn highlights the need for more accurate modelling of the normal and tangential contact.

The new normal contact method MIM-1D in its standard version has been shown to improve the contact patch and normal pressure computations as compared to existing fast approaches. While using a semi-analytical methodology permits a reduction in the computation time, re-evaluating the matrix of influence coefficients at each iteration implies that up to two-thirds of the CPU time can be spent on the Gaussian quadrature. As a further perspective, it is proposed to replace the numerical integration in MIM-1D with heuristic expressions to further improve computational efficiency compared to CONTACT. Initial results are found to be promising and may merit further investigations to improve precision and robustness.

In the present work, MIM-1D is only tested for a limited number of contact cases, commonly used to validate new wheel-rail contact methods. The future research direction should be to investigate and incorporate improvements using different wheel-rail profiles, kinematic conditions, axle loads etc. Even though Kalker's variational theory is often used as the reference for any new rolling contact model, the shortcomings of the half-space approach are again highlighted in the simulation cases tested in the collaborative study. The calculation of the contact patch half-length in MIM-1D has already been shown to be a delicate subject when treating unique contact configurations such as a sharp-edge contact and warrants further studies. A future research aspect may also be to carry out detailed studies of the contact scenarios encountered in S&C simulations using FE tools, which can then serve as the reference for the development of fast methods adapted to treat these cases. The ultimate validation of a complete rolling contact model, such as MIM-1D + FASTSIM SH or MIM-1D + FaStrip, can only be through its implementation in an MBS software. This would enable wide-ranging studies such as online damage analysis and wear predictions and thus, not be limited to the theoretical cases generally investigated in the literature.

110 4 . 1 57 2. 1 98 xv 3 . 8 107 4. 1

 4157198381071 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.2 The MBS software VOCO . . . . . . . . . . . . . . . . . . . . 111 xii List of Figures 1.1 Example of the inertial and body coordinate systems . . . . . 1.2 Rotations between different coordinate systems . . . . . . . . 1.3 Origin of the mobile track coordinate system and the cant angle ψ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 The hunting phenomenon . . . . . . . . . . . . . . . . . . . . 1.5 Wheel climb and the contact forces . . . . . . . . . . . . . . . 1.6 A standard vehicle model . . . . . . . . . . . . . . . . . . . . 1.7 A standard ballasted track . . . . . . . . . . . . . . . . . . . . 1.8 Description of the different creepages for railway applications 1.9 Example of a Hertzian contact ellipse with semi-axes a and b 1.10 Comparison of the normal contact results obtained for worn wheel-rail profiles using a FE software, CONTACT and Hertz' theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 The half-space approximation . . . . . . . . . . . . . . . . . .1.12 The general contact problem . . . . . . . . . . . . . . . . . . 1.13 Contact between two spheres with radii R 1 and R 2 . . . . . 1.14 Interpenetration region and definition of the origin in Kik-Piotrowski's method . . . . . . . . . . . . . . . . . . . . . . . 1.15 The semi-Hertzian approach when considering an elliptical contact patch . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.16 Determination of the virtual interpenetration g(y) for wheelrail profiles in STRIPES . . . . . . . . . . . . . . . . . . . . . 1.17 Distribution of tangential tractions in Carter's theory . . . . 1.18 Creep curve from Carter's theory . . . . . . . . . . . . . . . . 1.19 Stick-slip zone division according to Johnson, and Vermeulen and Johnson theories . . . . . . . . . . . . . . . . . . . . . . . 1.20 Stick-slip zone division and the traction distribution over each strip according to strip theory . . . . . . . . . . . . . . . . . . 1.21 Stick-slip zone division and the traction distribution over each strip according to the simplified theory, with φ = 0 . . . . . . xiv 1.22 Stick-slip zone division and the traction distribution over each strip according to Polach's theory without spin . . . . . . . . Two elastic bodies in contact . . . . . . . . . . . . . . . . . . 62 2.2 Contact area C divided into strips, and the normal pressure distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.3 MIM-1D algorithm with a given rigid body approach δ . . . . 70 2.4 Separation profiles of the two contacting surfaces . . . . . . . 71 2.5 Results for a Hertzian profile using existing methods and MIM-1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 2.6 Results for a non-Hertzian profile using existing methods and MIM-1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 2.7 Convergence of the proposed algorithm for different potential contact areas . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 2.8 Nominal profiles defined in the tangent plane (top) and the curvatures (bottom) for the wheel-rail pair S1002-UIC60 1:40 75 2.9 Contact patch and the maximum pressure distribution for wheel-rail contact cases using MIM-1D and other methods . . 76 2.10 Comparison of the relative error in the contact area with reference results from CONTACT . . . . . . . . . . . . . . . . 77 2.11 Comparison of the contact patch, and the maximum pressure distribution using different alternatives for the pressure in the rolling direction in MIM-1D . . . . . . . . . . . . . . . . . . . 81 2.12 Comparison of the contact patch area, and the maximum normal pressure for a range of wheel positions using different alternatives for the pressure in the rolling direction in MIM-1D 82 3.1 Elastic body represented as a set of independent springs . . . 87 3.2 The discretisation in FASTSIM . . . . . . . . . . . . . . . . . 90 3.3 Description of the longitudinal creepage ν xi for a strip located in y i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3.4 Tangential tractions, their directions, and the stick-slip zone division using different settings for the traction bound . . . . 94 3.5 Absolute error in the normalised contact forces computed with the FASTSIM SH algorithm w.r.t. CONTACT . . . . . . 95 3.6 The normal pressure distribution over the contact patch using CONTACT and MIM-1D for given wheel positions . . . . . . 97 3.7 Tangential tractions, their directions, and the stick-slip zone division using CONTACT, FASTSIM SH , and FaStrip for given wheel positions . . . . . . . . . . . . . . . . . . . . . . . Longitudinal and lateral shear stresses for the wheel position ∆y = -1 mm at the lateral coordinate y = 0 mm . . . . . . . 100 3.9 Longitudinal and lateral shear stresses for the wheel position ∆y = 0 mm at the lateral coordinate y = -11.59 mm . . . . . 100 3.10 Derived RCF damage index as a function of T γ . Adapted from [129] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 3.11 Normal and tangential contact results using STRIPES with FASTSIM SH for the given wheel positions . . . . . . . . . . . 102 3.12 The energy index EI distribution over the contact patch using CONTACT, MIM-1D, and STRIPES for the wheel positions: (a) ∆y = -1 mm, and (b) ∆y = 0 mm . . . . . . . . . . . . . 104 3.13 Horizontal curvature for the left-curved track, wheel position over the rail ∆y, and the normal contact forces N as functions of the track position . . . . . . . . . . . . . . . . . . . . . . . 105 3.14 The longitudinal and lateral creepages, and the spin creepage as functions of the track position . . . . . . . . . . . . . . . . 106 3.15 The longitudinal and lateral creep forces as functions of the track position . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 3.16 Evolution of the wear number T γ as a function of the track position . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Definition of a vehicle model in VOCO . . . . . . . . . . . . . 114 4.2 Description of the different properties pre-tabulated in the contact tables for a given position of the wheel ∆y . . . . . . 116 4.3 Example of profile discretization and separation into different bodies for a crossing in VOCO . . . . . . . . . . . . . . . . . 117 4.4 Wheel vertical kinematics with respect to the rail . . . . . . . 118 4.5 Variation of t z as a function of the distance in the case of crossing, for ∆y = 0 . . . . . . . . . . . . . . . . . . . . . . . 119 4.6 Layout, components and nomenclature for a standard righthand side S&C . . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.7 Topology of the co-running track model used in the S&C benchmark with masses and bushing elements . . . . . . . . . 123 4.8 Modelling of the normal contact between the check rail and the wheel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.9 Modelling of the normal contact between the check rail and the wheel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 4.10 Modelling of the contact stiffness in VOCO . . . . . . . . . . 126 xvi 4.11 Comparison of the vertical wheel-rail Q forces obtained in a crossing panel when using different settings for the contact stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 4.12 Locations of the varying rail geometry for the benchmark simulation cases . . . . . . . . . . . . . . . . . . . . . . . . . . 129 4.13 Contact position for the Run 1, considering the through route with the British S&C . . . . . . . . . . . . . . . . . . . . . . . 130 4.14 Vertical wheel-rail Q forces for the Run 1, considering the through route with the British S&C . . . . . . . . . . . . . . 132 4.15 Vertical offset of bodies in determination of contact points due to the flexible track model . . . . . . . . . . . . . . . . . 132 4.16 Contact position for the Run 2, considering the diverging route with the British S&C . . . . . . . . . . . . . . . . . . . 133 4.17 Contact position for the Run 2, considering the diverging route with the British S&C . . . . . . . . . . . . . . . . . . . 133 4.18 Contact position for the Run 3, considering the through route with the Swedish S&C . . . . . . . . . . . . . . . . . . . . . . 134 4.19 Contact position for the Run 4, considering the diverging route with the Swedish S&C . . . . . . . . . . . . . . . . . . . 134 4.20 Sharp-edge contact between the switch rail and wheel in VOCO for Run 3, when taking the through route with Swedish S&C . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 4.21 Contact position for the Run 7, considering the through route with the Swedish S&C . . . . . . . . . . . . . . . . . . . . . . 136 4.22 Vertical wheel-rail Q forces for the Run 7, considering the through route with the Swedish S&C . . . . . . . . . . . . . . 136 4.23 Contact position for the Run 8, considering the diverging route with the Swedish S&C . . . . . . . . . . . . . . . . . . . 137 4.24 Lateral wheel-rail Y forces for the Run 8, considering the diverging route with the Swedish S&C . . . . . . . . . . . . . 137 4.25 3D geometry interpolation for variable rail profiles . . . . . . 141 4.26 Creepage assessed at: (a) the centre of the contact patch, (b) each cell of a discretized a contact patch and (c) each strip of a contact patch . . . . . . . . . . . . . . . . . . . . . . . . . . 142 4.27 Schematic representation of the different run cases . . . . . . 143 4.28 Rail profiles used in Run 1 and Run 2 . . . . . . . . . . . . . 145 4.29 Positions of the contact patch on the rail for Run 1 . . . . . . 145 4.30 Area of the contact patches observed in the outer wheels for Run 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 xvii 4.31 Contact patches of the outer trailing wheel for Run 1 obtained with VOCO, MIM-1D, MUBODyn and CONTACT for different track positions . . . . . . . . . . . . . . . . . . . 147 4.32 Normal pressure distribution in the contact patches of the outer trailing wheel for Run 1 obtained with VOCO, MIM-1D, and CONTACT for different track positions . . . . . . . 147 4.33 Contact forces in the outer tailing wheels for Run 1 . . . . . . 149 4.34 Positions of the contact patch on the rail for Run 2 . . . . . . 150 4.35 Contact patches of the outer trailing wheel for Run 2 obtained with VOCO, MIM-1D, MUBODyn and CONTACT for different track positions . . . . . . . . . . . . . . . . . . . 151 4.36 Run 2 results for VOCO at 25 m: (a) lateral B curvature with and without smoothing, and (b) corresponding contact patches including CONTACT results . . . . . . . . . . . . . . 151 4.37 Normal pressure distribution in the contact patches of the outer trailing wheel for Run 2 obtained with VOCO, MIM-1D, and CONTACT for different track positions . . . . . . . 152 4.38 Contact forces in the outer tailing wheels for Run 2 . . . . . . 153 4.39 Normal impact forces in the outer wheel of the leading wheelset at different intervals for Run 3 . . . . . . . . . . . . 155 D.1 The bogie model . . . . . . . . . . . . . . . . . . . . . . . . . 188 D.2 The co-running track model . . . . . . . . . . . . . . . . . . . 188 E.1 Schematic representation of the wheel profile parametrization 192 E.2 Schematic representation of the rail profile parametrization . 193 E.3 Schematic representation of the crossing layout parametrization194 E.4 Cross-sections of the parametrized crossing panel . . . . . . . 195 xviii List of Tables

Figure 1 . 1 :

 11 Figure 1.1: Example of the inertial and body coordinate systems

Figure 1 . 2 :

 12 Figure 1.2: Rotations between different coordinate systems

Figure 1 . 4 :

 14 Figure 1.4: The hunting phenomenon

Figure 1 . 5 :

 15 Figure 1.5: Wheel climb and the contact forces

Figure 1 . 6 :

 16 Figure 1.6: A standard vehicle model

Figure 1 . 8 :

 18 Figure 1.8: Description of the different creepages for railway applications: (a) longitudinal creepage, (b) lateral creepage, and (c) spin creepage

Figure 1 . 9 :

 19 Figure 1.9: Example of a Hertzian contact ellipse with semi-axes a and b

Figure 1 . 10 :

 110 Figure 1.10: Comparison of the normal contact results obtained for worn wheel-rail profiles using a FE software, CONTACT and Hertz' theory.From[START_REF] Telliskivi | Contact mechanics analysis of measured wheel-rail profiles using the finite element method[END_REF]

Figure 1 . 13 :

 113 Figure 1.13: Contact between two spheres with radii R 1 and R 2

Figure 1 . 14 :

 114 Figure 1.14: Interpenetration region and definition of the origin in Kik-Piotrowski's method

Figure 1 .

 1 Figure 1.15: The semi-Hertzian approach when considering an elliptical contact patch

Figure 1 . 17 :

 117 Figure 1.17: Distribution of tangential tractions in Carter's theory

Figure 1 . 18 :

 118 Figure 1.18: Creep curve from Carter's theory

Figure 1 . 20 :

 120 Figure 1.20: Stick-slip zone division and the traction distribution over each strip according to strip theory

Figure 1 . 22 :

 122 Figure 1.22: Stick-slip zone division and the traction distribution over each strip according to Polach's theory without spin

Figure 2 . 1 :⋆ 1 (

 211 Figure 2.1: Two elastic bodies in contact

Figure 2 . 2 :

 22 Figure 2.2: Contact area C divided into strips, and the normal pressuredistribution. Adapted from[START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with nonelliptical contact area[END_REF]

Figure 2 . 4 :

 24 Figure 2.4: Separation profiles of the two contacting surfaces

Figure 2 . 5 :Figure 2 . 6 :

 2526 Figure 2.5: Results for a Hertzian profile using existing methods and MIM-1D

Figure 2 . 7 :

 27 Figure 2.7: Convergence of the proposed algorithm for different potential contact areas

Figure 2 .

 2

Figure 2 . 8 :

 28 Figure 2.8: Nominal profiles defined in the tangent plane (top) and the curvatures (bottom) for the wheel-rail pair S1002-UIC60 1:40

Figure 2 . 9 :Figure 2 . 9 :Figure 2 . 10 :

 2929210 Figure 2.9: Contact patch (left), and the maximum pressure distribution p 0 (y) (right) for wheel-rail contact cases using MIM-1D and other methods (from top to bottom): (a) ∆y = -1 mm, (b) ∆y = 0 mm, (c) ∆y = 1 mm

Figure 2 . 11 :Figure 2 .

 2112 Figure 2.11: Comparison of the contact patch (left), and the maximum pressure distribution p 0 (y) (right) for wheel-rail contact cases (from top to bottom): (a) ∆y = -1 mm, (b) ∆y = 0 mm, (c) ∆y = 1 mm, using different alternatives for the pressure in the rolling direction in MIM-1D

Figure 2

 2 Figure 2.12: Comparison of the contact patch area (left), and the maximum normal pressure (right) for a range of wheel positions using different alternatives for the pressure in the rolling direction in MIM-1D

Figure 3 . 1 :

 31 Figure 3.1: Elastic body represented as a set of independent springs.Adapted from[21, p. 80]

3 2 c

 2 23 φ , (3.9) where G is the modulus of rigidity, and c 11 , c 22 , and c 23 are Kalker's coefficients, which are functions of the Poisson's ratio ν and the ellipse ratio b/a. The linear theory is exact for infinitesimal values of creepages. By multiplying the Equation (3.8) and Equation (3.9) for the creep forces by the flexibility parameter L, and by comparing them with the Equation (3.1), Equation (3.6), and Equation (3.7), we obtain, not one, but instead three expressions for L:

Figure 3 . 2 :

 32 Figure 3.2: The discretisation in FASTSIM

. 21 )Figure 3 . 3 :

 2133 Figure 3.3: Description of the longitudinal creepage ν xi for a strip located in y i . From [76]

Figure 3 . 4 :

 34 Figure 3.4: Tangential tractions, their directions, and the stick-slip zone division using different settings for the traction bound: (a) CONTACT, (b) FASTSIM SH with semi-ellipsoidal traction bound, (c) FASTSIM SH with parabolic traction bound, and (d) FASTSIM SH with weighted parabolic traction bound

Figure 3 . 5 :

 35 Figure 3.5: Absolute error in the normalised contact forces computed with the FASTSIM SH algorithm w.r.t. CONTACT

Figure 3 . 6 :

 36 Figure 3.6: The normal pressure distribution over the contact patch using CONTACT and MIM-1D for given wheel positions: (a) ∆y = -1 mm, and (b) ∆y = 0 mm

Figure 3 . 7 :

 37 Figure 3.7: Tangential tractions, their directions, and the stick-slip zone division using CONTACT, FASTSIM SH , and FaStrip for given wheel positions: (a) ∆y = -1 mm, and (b) ∆y = 0 mm

FaStripFigure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: Longitudinal and lateral shear stresses for the wheel position ∆y = -1 mm at the lateral coordinate y = 0 mm

Figure 3 . 11 : 102 3. 3 .

 3111023 Figure 3.11: Normal and tangential contact results using STRIPES with FASTSIM SH for the wheel positions: (a) ∆y = -1 mm, and (b) ∆y = 0 mm

Figure 3 . 12 :Figure 3 . 13 :

 312313 Figure 3.12: The energy index EI distribution over the contact patch using CONTACT, MIM-1D, and STRIPES for the wheel positions: (a) ∆y = -1 mm, and (b) ∆y = 0 mm

Figure 3 . 14 :Figure 3 . 15 :Figure 3 . 16 :

 314315316 Figure 3.14: The longitudinal and lateral creepages, and the spin creepage as functions of the track position

Figure 4 . 1 :

 41 Figure 4.1: Definition of a vehicle model in VOCO

Figure 4 . 2 :

 42 Figure 4.2: Description of the different properties pre-tabulated in the contact tables for a given position of the wheel ∆y

Figure 4 . 3 :

 43 Figure 4.3: Example of profile discretization and separation into different bodies for a crossing in VOCO

Figure 4 . 4 :

 44 Figure 4.4: Wheel vertical kinematics with respect to the rail

Figure 4 . 5 :

 45 Figure 4.5: Variation of t z as a function of the distance in the case of crossing, for ∆y = 0

Figure 4 . 6 :

 46 Figure 4.6: Layout, components and nomenclature for a standard right-hand side S&C. From [135]

4. 3 . 1 . 1

 311 Implementation of the co-running track modelAs shown in Table

Figure 4 . 7 :

 47 Figure 4.7: Topology of the co-running track model used in the S&Cbenchmark with masses and bushing elements. From[START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]

Figure 4 . 8 :

 48 Figure 4.8: Modelling of the normal contact between the check rail and the back flange of the wheel. From [141]

Figure 4 . 9 :

 49 Figure 4.9: Modelling of the normal contact between the check rail and the wheel

Figure 4 . 10 :

 410 Figure 4.10: Modelling of the contact stiffness in VOCO: (a) standard model, (b) new model with Hertzian stiffness

( 4 . 7 ) 61 Time

 4761 The stiffness is set to a typical value of approximately 0.7E11 N/m 0.57 0.575 0.58 0.585 0.59 0.595 0.6 0.605 0.K z = 500 kN/mm K z = 1000 kN/mm k h = 0.7 10 11 N/m 3/2

Figure 4 . 11 :

 411 Figure 4.11: Comparison of the vertical wheel-rail Q forces obtained in a crossing panel when using different settings for the contact stiffness

Figure 4 . 12 :

 412 Figure 4.12: Locations of the varying rail geometry for the benchmark simulation cases in Table 4.3: Runs 1 & 3 (top left), Runs 2 & 4 (top right), Runs 5 & 7 (bottom left) and Runs 6 & 8 (bottom right). From [102]

Figure 4 . 14 :

 414 Figure 4.14: Vertical wheel-rail Q forces for the Run 1, considering the through route with the British S&C. From [102]

Figure 4 . 15 :

 415 Figure 4.15: Vertical offset of bodies in determination of contact points due to the flexible track model

Figure 4 . 16 :

 416 Figure 4.16: Contact position for the Run 2, considering the diverging route with the British S&C. From [102]

Figure 4 . 17 :

 417 Figure 4.17: Contact patch size for the Run 2, considering the diverging route with the British S&C. From [102]

Figure 4 . 18 :

 418 Figure 4.18: Contact position for the Run 3, considering the through routewith the Swedish S&C. From[START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]

Figure 4 . 19 :

 419 Figure 4.19: Contact position for the Run 4, considering the diverging route with the Swedish S&C. From [102]

Figure 4 . 20 :

 420 Figure 4.20: Sharp-edge contact between the switch rail and wheel in VOCO for Run 3, when taking the through route with Swedish S&C

Figure 4 .

 4 21 and Figure 4.22 show some results for

Figure 4 . 21 :

 421 Figure 4.21: Contact position for the Run 7, considering the through routewith the Swedish S&C. From[START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]

Figure 4 . 22 :

 422 Figure 4.22: Vertical wheel-rail Q forces for the Run 7, considering the through route with the Swedish S&C. From [102]

Figure 4 . 23 :

 423 Figure 4.23: Contact position for the Run 8, considering the diverging route with the Swedish S&C. From [102]

Figure 4 . 24 :

 424 Figure 4.24: Lateral wheel-rail Y forces for the Run 8, considering the diverging route with the Swedish S&C (top: check rail, bottom: crossing rail). From[START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]

  2.2.2, and shown in Figure 4.25a.

Figure 4 .

 4 Figure 4.25: 3D geometry interpolation for variable rail profiles: (a) separation of profiles into bodies, and (b) using section breaks. Adapted from [102]

Figure 4 .

 4 26a. VOCO uses the CHOPAYA function, described in Section 4.2.2.2, to evaluate the creep forces via the local

Figure 4 . 26 :

 426 Figure 4.26: Creepage assessed at: (a) the centre of the contact patch, (b) each cell of a discretized a contact patch and, (c) each strip of a contact patch. From [145]

Figure 4 . 27 :

 427 Figure 4.27: Schematic representation of the different run cases: (a) conformal contact, (b) sharp-edge contact, and (c) impact load on the crossing nose. From [145]

Figure 4 . 28 :

 428 Figure 4.28: Rail profiles used in Run 1 and Run 2

Figure 4 . 29 :

 429 Figure 4.29: Positions of the contact patch on the rail for Run 1

Figure 4 . 30 :

 430 Figure 4.30: Area of the contact patches observed in the outer wheels for Run 1

Figure 4 . 34 .

 434 Similar to Run 1, a value near 0 mm indicates a contact near the tread, while a value near 40 mm indicates a flange contact, as seen in Figure4.28. As in the previous case, all three software predict similar

Figure 4 . 33 :

 433 Figure 4.33: Contact forces in the outer tailing wheels for Run 1

Figure 4 . 34 :

 434 Figure 4.34: Positions of the contact patch on the rail for Run 2

Figure 4 . 35 :Figure 4 . 36 :Figure 4 . 37 :

 435436437 Figure 4.35: Contact patches of the outer trailing wheel for Run 2 obtained with VOCO, MIM-1D, MUBODyn and CONTACT for different track positions: (a, c) 0 m, and (b, d) 25 m

Figure 4 . 38 :

 438 Figure 4.38: Contact forces in the outer tailing wheels for Run 2

Figure 4 . 39 :

 439 Figure 4.39: Normal impact forces in the outer wheel of the leading wheelset at different intervals for Run 3

  The rolling contact problem with Kalker's complete theory .1.3.1 The elastic half-space and quasi-identity . . . . . . . . 1.3.2 The general contact problem . . . . . . . . . . . . . . 1.4 Normal contact problem using fast approaches . . . . . . . . 1.4.1 Hertz' theory . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Virtual penetration methods for non-Hertzian contact 1.4.3 ANALYN . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Tangential contact problem in Hertzian conditions . . . . . . 1.5.1 Carter's theory . . . . . . . . . . . . . . . . . . . . . . 1.5.2 Johnson's theories . . . . . . . . . . . . . . . . . . . . 1.5.3 Kalker's theories . . . . . . . . . . . . . . . . . . . . . Heuristic solutions considering spin . . . . . . . . . . . 55 1.5.6 Lookup tables . . . . . . . . . . . . . . . . . . . . . . . 58 1.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 A new semi-analytical method for modelling the normal wheel-rail contact 60 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 2.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 2.2.1 Theoretical background . . . . . . . . . . . . . . . . . 61 2.2.2 Towards a reduced approach . . . . . . . . . . . . . . 63 2.2.3 Discrete problem . . . . . . . . . . . . . . . . . . . . . 64 2.2.4 Estimation of the contact patch half-length . . . . . . 67 2.2.5 Iterative resolution . . . . . . . . . . . . . . . . . . . . 69 2.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 71 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.2 Extension of FASTSIM for steady state non-Hertzian contact 86 3.2.1 The simplified theory . . . . . . . . . . . . . . . . . . 86 3.2.2 From FASTSIM over an ellipse to FASTSIM over strips 89 3.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 92 3.3.1 Numerical experiments to determine the traction bound 93 3.3.2 Validation of FASTSIM SH using CONTACT . . . . . 94 3.3.3 Combination of MIM-1D with a rolling contact model 96 3.3.4 Comparison of rolling contact models used for damage prediction . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

	Contents
	1.5.5 On the use of semi-analytical contact models in rolling
	contact conditions	85
	3.1 Case studies in railway dynamics
	xi

Introduction 1 Overview of wheel-rail contact models used with the multibody systems approach 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 The multibody systems approach for railway applications . . 1.2.1 Coordinate systems and associated degrees of freedom 1.2.2 Railway vehicle dynamics and motion scenarios . . . . 1.2.3 Vehicle and track modelling . . . . . . . . . . . . . . . 1.2.4 Wheel-rail contact modelling in the context of multibody systems . . . . . . . . . . . . . . . . . . . . 1.3 1.5.4 The FaStrip algorithm . . . . . . . . . . . . . . . . . . 2.3.1 Theoretical profiles . . . . . . . . . . . . . . . . . . . . 71 2.3.2 Wheel-rail contact . . . . . . . . . . . . . . . . . . . . 73 2.3.3 Computational cost . . . . . . . . . . . . . . . . . . . 75 2.4 Alternative pressure distributions in the rolling contact direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 2.4.1 Uniform pressure distribution . . . . . . . . . . . . . . 78 2.4.2 Parabolic pressure distribution . . . . . . . . . . . . . 79 2.5

Table 1 .

 1

	θ •	0	5	10	30	60	90 120	150	170	175	180
	A/B 0	0.0019 0.0077 0.0717 0.3333 1	3.0	13.93	130.6	524.6	∞
	n/m 0	0.0212 0.0470 0.1806 0.4826 1	2.0720 5.5380 21.26	47.20	∞
	m	∞ 11.238 6.612	2.731	1.486	1	0.7171 0.4931 0.311	0.2381 0
	r	0	0.2969 0.4280 0.7263 0.9376 1	0.9376 0.7263 0.4280 0.2969 0

1: Hertz' coefficients as functions of θ. From

[START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF]

Table 2 .

 2 1: Comparison of the CPU time with the different existing approaches ( * the total number of elements for CONTACT is n × n)

			Time (s)
	Mesh size n * STRIPES ANALYN MIM-1D CONTACT
	0.2 mm	100	1.5E-5	6.2E-5	2.0E-3	< 0.1
	0.08 mm 250	3.1E-5	9.4E-5	2.9E-2	0.2
	0.04 mm 500	9.3E-5	1.3E-4	9.2E-2	1.3

Table 3 .

 3 1: Comparison of absolute normalised longitudinal and lateral creep forces using different rolling contact models presented in Figure

	3.3. Results and discussion				Chapter 3
			FASTSIM SH FaStrip CONTACT
	∆y = -1 mm	|f x | |f y |	0.344 0.052	0.321 0.059	0.337 0.063
	∆y = 0 mm	|f x | |f y |	0.409 0.051	0.397 0.059	0.407 0.066

Table 4 .

 4 1: Overview of the different features in VOCO

Table 4 .

 4 2: Overview of the different output provided by VOCO

	Scope	Output quantities
	Global	Wheel-rail contact forces
		Displacements, velocities, and accelerations of each body
		Forces in the connecting links
	Local	Normal and tangential contact forces per strip for each wheel
		Longitudinal, lateral and spin creepages per strip for each
		wheel
		Contact patches and contact stresses

Table 4 .

 4 3: S&C benchmark simulation cases. Adapted from[START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF]

	Run S&C	Panel	Route	Speed (km/h)
	1	British	Switch	Through	100
	2	British	Switch	Diverging	43
	3	Swedish Switch	Through	160
	4	Swedish Switch	Diverging	80
	5	British	Crossing Through	100
	6	British	Crossing Diverging	43
	7	Swedish Crossing Through	160
	8	Swedish Crossing Diverging	80
	9	Identical to Run 2 but with a constant 56E1 profile
		replacing the stock rail, and no switch rail

Excluding the dofs linked to the rolling contact model

In practice, the expression for νx also contains an additional term that takes into account the advance of the contact point as a function of the yaw angle α

1.3. The rolling contact problem with Kalker's complete theoryChapter 1

δtz may be considered equivalent to ϵδ when the contact angle is negligible

As long as contact patches are assumed to be elliptical

2bi is the same as δyi used in STRIPES in Chapter 1

The Aij terms here correspond to the normal influence coefficients Azz ij

5000 simulations in total

This is only valid for small cracks, as it ignores the characteristic that crack growth rate is influenced by crack size

, which

This is also highlighted in Figure4.5 showing the variation of the wheel vertical kinematics tz with respect to the rail in the crossing panel

Conclusions & future work

Chapter 2

A new semi-analytical method for modelling the normal wheel-rail contact On the use of semi-analytical contact models in rolling contact conditions Chapter 4 Case studies in railway dynamics

Appendices

Appendix A

Theory of Boussinesq and Cerruti

The normal stress distribution and normal displacements in an elastic halfspace due to a point loading was studied in detail by Boussinesq in 1885 [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF]. A similar study had been done three years earlier by Cerruti who investigated the action of tangential surface tractions [START_REF] Cerruti | Richerche intorno all equilibrio de corpi elastici isotropi[END_REF] on the stress distribution and displacements. Boussinesq and Cerruti's analytical solutions use the theory of potentials to give the surface displacements at any point of the semiinfinite solid due to the action of a concentrated load acting on the semiinfinite solid. In the surface mechanical form, the relation between the displacements u and the stresses p can be given as:

where x = [x, y, z], x' = [x ′ , y ′ , z ′ ], p = [p n , p x , p y ], and A(x,x') is the influence function matrix describing the displacement difference at x arising due to a unit load in x'. The influence function depends on the geometry of the body. Closed form solutions of the Boussinesq's problem have been determined for various simple forms such as in the works of Love [START_REF] Love | Boussinesq's problem for a rigid cone[END_REF] and Sneddon [START_REF] Sneddon | Boussinesq's problem for a rigid cone[END_REF] for a rigid cone, and Sneddon for flat-ended cylinders [START_REF] Sneddon | Boussinesq's problem for a flat-ended cylinder[END_REF]. More importantly, the influence functions for the semi-infinite half-space can be calculated analytically [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF][START_REF] Love | The stress produced in a semi-infinite solid by pressure on part of the boundary[END_REF]. One may refer to the book by Johnson [START_REF] Johnson | Contact mechanics[END_REF] for a synopsis of the results available in literature.

For the normal and tangential contact problem, the influence functions are defined on the contacting surface, with x = [x, y, 0]. Under the influence of a purely normal load p n acting on the contact surface C, the displacements at all the points can be given as:

where the subscripts 1 and 2 denote the two bodies in contact. The relative displacements due to the normal load can subsequently be expressed by using the relevant expressions of the influence functions from the theory of Boussinesq and Cerruti:

x -x ′ r 2 p n (x')dC , (A.5)

, G is the combined modulus of rigidity, and K is the difference parameter. Following the same methodology, the relative displacements on the surface due to pure longitudinal are calculated using the corresponding influence functions:

and for the pure lateral loading as:

The total surface displacement differences are finally obtained through the superimposition of each of the individual displacement differences:

From the above expressions, the normal and tangential contact can be observed to be coupled through the parameter K. When the two bodies have the same material properties, K is zero and the bodies are said to be quasi-identical. In these conditions, the normal and tangential problems become decoupled and the above equations can be written in the following form:

Analytical expressions for the integrals in the above relations can be found in [START_REF] Johnson | Contact mechanics[END_REF] or the appendix in [25, p. 150].

Contact patch half-length using ANALYN

The expression for the contact patch boundaries 1 in ANALYN [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF] is:

where

Here, g i is the interpenetration, and the term β i h i takes the surface deformation into account analytically, as opposed to methods based on the virtual penetration where the deformation is neglected. The coefficients α i and β i are defined as:

3)

where A i and B i are the relative longitudinal and lateral curvatures, respectively, while m i , n i , and r i are non-dimensional Hertzian coefficients calculated using the local curvatures. The maximum pressure value p 0 i is:

Squaring Equation (B.1) and taking into account the expression for α i , we

1 ai represents the discretised value of the continuous value a(yi)

have:

Combining Equation (B.5) and Equation (B.6), the half-length of the contact patch a i can finally be written as:

which is the same as Equation (2.26).

Appendix C

Alternate expression for the contact patch half-length used in MIM-1D

The assessment of a i is critical in order to get consistent results in MIM-1D. The chosen expressions Equation (2.26) or Equation (2.31) make use of the local curvatures as well as the pressure (or normal force) distribution obtained within the contact patch, and provide best results in the studied cases as compared to other tested expressions. Still, this choice may considerably deviate from the CONTACT results in profiles showing a large gradient of surface curvature, as in the case of sharp-edge contact. This is a shortcoming of the method, even though such profiles should be considered with care, either due to the high values of pressures (leading to plastification) or due to the use of methods based on the half-space assumption.

Another expression of the half-length is proposed here, which is not based on the local values of the curvature, but rather on the assessment of the curvatures of an equivalent ellipse, as done in the Kik-Piotrowski method [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF]. The half-length in the KP method is first approximated, as described in Section 1.4.2.1, using the interpenetration region:

In MIM-1D, ϵ is not computed, but the limits of the contact patch in the plane x = 0 is known. Let y l and y r be the limits of the contact patch, and y m be the location of the minimum of -h(y i ) between y l and y r . Then, g(y i ) may be approximated using:

where ο is a very small positive value to avoid a zero value for a i in Equation (C.1), which would make the influence coefficient matrix C singular. Comparing Equation (C.2) with Equation (C.4), it is clear that ϵδ can be evaluated as:

As mentioned in Chapter 1, in every VP method, a correction must be made in order to get the aspect ratio predicted by the Hertzian theory. In the KP method [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF], the correction is based on the estimation of the aspect ratio of an equivalent ellipse. The mean lateral B curvature is evaluated using the width of the contact patch W = y r -y l , such that:

The longitudinal curvature is corrected using the semi-Hertzian theory:

where m and n are the Hertzian coefficients given by Table 1.1 for the mean curvature ratio A/B. A c then replaces the A value in Equation (C.1) to evaluate the half-lengths a i .

Vehicle and track models

The following sections present the vehicle and track models used in Section 4.4, and are adapted from [START_REF] Magalhaes | Wheel-rail contact models in the presence of switches and crossings[END_REF].

D.1 Vehicle model

A simple bogie is considered as the vehicle model, which consists of a bogie frame supported by two wheelsets through the primary suspension as depicted in Figure D.1. The parameters used for this model are listed in Table D.1 and are based on the passenger vehicle used in the Manchester Benchmarks [START_REF] Iwnicki | Manchester benchmarks for rail vehicle simulation[END_REF], with the exception that only half of the car body mass is considered, and the series stiffness of the primary suspension elements are removed.

D.2 Track model

The track model is a co-running model, similar to the one used in the S&C benchmark [START_REF] Bezin | Multibody simulation benchmark for dynamic vehicle-track interaction in switches and crossings: results and method statements[END_REF], and is shown in Appendix E

Parametrized wheel and rail profiles

The following sections present the parametrization of the wheel and rail profiles, as well as the crossing layout used in Section 4.4, and are adapted from [START_REF] Magalhaes | Wheel-rail contact models in the presence of switches and crossings[END_REF].

E.1 Wheel profile

The wheel profile is parametrized using the curved sections of three ellipses, as depicted in Figure E.1. Each ellipse is defined by the two semi-axes a i and b i , where i corresponds to the respective part of the wheel profile, with 'T' for the tread, 'C' for the concave part, and 'F' for the flange. The union between ellipses is such that it always ensures the continuity of the first derivative. The different parameters used to define the wheel profile, and their corresponding values, are given in Table E.1.

E.2 Rail profile

The rail profiles are parametrized via the union of two ellipses that represent the head or the top of the rail, and two straight lines that represent the lateral faces, as shown in Figure E.2. Similar to the wheel parametrization, each ellipse is defined by the two semi-axes a i and b i , where i corresponds to the respective side of the rail profile, with 'L' and 'R' representing the left and tight sides respectively. The continuity of the first derivative is ensured at all points. The different parameters used to define the rail profile, and their corresponding values, are given in Table E.2.