Comparison of numerical and computational aspects between two constraint-based contact methods in the description of wheel/rail contacts

The numerical and computation aspects of the Knife-edge Equivalent Contact (KEC) constraint and lookup table (LUT) methods are compared in this paper. The LUT method implementation uses a penetration-based elastic contact model for the flange and a constraint-based formulation at the wheel tread. For the KEC method, where an infinitely narrow rail contacts an equivalent wheel, regularization of the tread-flange transition is adopted to simultaneously account for tread and flange contacts using constraints. A comparison between the two methods is carried out using well-known numerical integrators to show the applicability and limitations of both methods. Two fixed-step-size integrators, the explicit Runge–Kutta (RK4) and the predictor–corrector Adam–Bashforth–Moulton (ABM) methods, and two variable-step-size Matlab built-in function integrators, the explicit ode45\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ode45$\end{document} and implicit ode15s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ode15s$\end{document}, were applied to get the numerical solutions to the dynamic problems and study the relative numerical performance of the two contact description methods. To complete the railway vehicle model, both contact methods were implemented for the multibody model of a benchmark railway vehicle (the Manchester wagon 1). Numerical results were obtained for different railway tracks with and without irregularities. Profiles of the S1002 wheel and LB-140-Area rail, which demonstrate the two-point contact phenomenon, were considered. Both methods were implemented in Matlab and validated against commercial simulation software. The kinematic results for both approaches show good agreement, but the KEC method was up to 20% more efficient than the LUT method regardless of integrator used.

Vector of generalized wheel-rail normal contact forces, generalized wheel-rail normal contact forces at tread and flange Q tang Vector of generalized tangential tread and flange forces R Radius of curved track [m] R t Position vector of a track frame with respect to the global frame R t x , R t y , R t z Position components of the position vector R t in X, Y , and Z direction R bti Position vector of the body track frame with respect to the global frame R ijt Relative distance vector from body track frames it to jt R P Absolute position vector of an point P in global framē r ir Position vectors of rail irregularity vectors of rail profile frames with respect to ideal rail profile frames r rp Position vectors of ideal railhead frames with respect to track frame r i P Position vector of point P that belongs to body i with respect to the global framē r i Position vector of body frame with respect to its the body track framē r wi Position vector of wheelset frame with respect to the wheelset track framē

Introduction
A railway vehicle is a complex mechanical system that consists of a large number of bodies, wheel/rail contacts, and complicated suspension elements. Because of the high computational loads, the real-time simulation of the railway vehicle is challenging [1], especially with respect to the wheel-rail contact [2]. Computationally efficient and physically accurate simulations can be used for a wide variety of applications, such as optimization problems and embedded models. Moreover, to analyze the validity and efficiency of different numerical methods, benchmark problems can be found in the literature. Examples include the Manchester benchmark [3] and the switches and crossings benchmarks [4].
In computational railway dynamics, a number of formalisms have been introduced to analyze wheel/rail contact. Among these, two well-known approaches are commonly used in multibody railway simulations. The first is the elastic approach [5][6][7], where interpenetration between wheel/rail surfaces is allowed, and normal contact forces are computed based on interpenetration and Hertzian parameters. However, the large magnitude Hertzian parameters that come about when applying the elastic approach can result in a stiff system of ordinary differential equations. Therefore, implicit integrators with small time steps are used, which increases computational cost. In [8], Liu and Bruni compared Hertzian and non-Hertzian wheel/rail contact models for multibody simulation. In [6], the Hertzian and non-Hertzian models are studied and compared to give the insight of proper selection of the parameters for better computational accuracy and efficiency. The second approach is the constraint method, which is also called Lagrange multiplier method. In this approach, the contact between wheel and rail is described with a set of kinematic constraint equations [9] Fig. 1 (a) Real S1002 wheel profile and LB.140-AREA rail profile geometry, (b) contact points on real wheel/rail profiles and corresponding KEC-equivalent wheel/rail profiles, the real wheel/rail profiles correspond to the profiles in the dashed boxes on the left to ensure both surfaces are in contact without penetration or separation. In this way, the normal contact forces can be computed using the Lagrange multipliers that are associated with the constraint equations.
Two constraint-based formulations for the wheel/rail contact in a multibody simulation have been studied here. The first is the LUT method. The LUT method is an offline constraint approach where the contact detection part of the analysis is carried out in a preprocessing stage and the solution is stored in a lookup table. In the subsequent dynamic simulations, the contact points are obtained by interpolating between the stored data points [10][11][12]. In the work of [10], an offline contact lookup table is used to predict tread contact point and an online elastic approach is used to predict flange contact. This algorithm was later extended to cover the combination of nodal and nonconformal contact detection to solve contact point jumps that can occur in turnouts [11]. The authors of [12] propose a constraint contact LUT method that accounts for track irregularities with two entries, lateral displacement, and track gauge variation. Using the Kalker book of tables for non-Hertzian contact (KBTNH) to build the lookup table is proposed [13] to generate the fast creep force. In that study, a regularization of the nonelliptical contact patch is introduced to solve the wheel/rail normal contact problem with non-Hertzian methods. In addition, the KBTNH is extended to include the full symmetry relations for creep forces and moment [14]. A comprehensive analysis to improve accuracy and efficiency when interpolating the lookup tables with multiple input parameters is presented in [15].
The second approach is the KEC method. The KEC method used here is an online constraint approach that was first proposed [16] to model wheel/rail contact. In the KEC method, the wheel/rail profile combination (see Fig. 1(a)) is established using an equivalent wheel profile (see solid lines in Fig. 1(b)) in contact with an infinitely narrow rail, which yields the equivalent allowable relative motion. This equivalent profile combination produces the same wheel/rail contact kinematics as the real wheel/rail profiles. As shown in Fig. 1(b), the contact points for a set of discrete values of the wheelset lateral displacement are located at the real wheel/rail profiles (dotted-dashed lines in Fig. 1(b)). The wheelset with wheels S1002 profile and rails LB.140-AREA profile are considered in the figure, and the geometry parameters for both wheel/ rail profiles are listed in Table 1. Accordingly, the corresponding contact points can be found on the KEC-equivalent profiles. See the solid lines in Fig. 1(b). The exact position of the contact points can be determined from the online solution of the KEC constraints. In addition, [17] proposed a regularization method to apply a continuous transition from tread to flange between each wheel/rail pair for the dynamic simulation of a two-point contact scenario. This regularization method can help avoid the finite contact point jumps that occur between tread and flange. Active research to develop computationally efficient models for multibody simulation of railroad vehicles is ongoing. Focus areas include linearizing the equations of motion [18], establishing an appropriate time integration method [18], efficiently detecting wheel/rail contact, and formalizing corresponding contact computation [10,19,20]. An explicit predictorcorrector integration method has been developed to solve the large-scale coupling dynamics of a railway vehicle and track system using the explicit method as the predictor and the implicit Newmark method as a corrector [21]. A quasielastic contact model proposed in [19] employs a two-dimensional spline approximation of contact solutions and stores the spline coefficients as a table in a pre-processing stage. The method is implemented in the wheel/rail module of the commercial software SIMPACK. In addition, a modal substructuring approach is proposed in [22] to improve the computational efficiency of a vehicle-track simulation by using a reduced number of modal coordinates. A hybrid contact detection algorithm is developed in [10] to eliminate the online iterative search for the second contact point. Moreover, a differential contact model is presented in [20] to reduce the computation cost by integrating the differential modeling of the wheel/rail contact problem with multibody modeling of the railway dynamics. The work reported in [18] proves that a simplified model based on weakly coupled lateral and vertical dynamics solves about four times faster than the full three-dimensional (3D) coupled dynamic model.
Several works have been published to study the pros and cons of the LUT and KEC methods with respect to computational efficiency. In the case of vehicle ML95 (Lisbon subway company) negotiating a curved track without irregularities in [12], a numerical solution is reached six times faster using the LUT method instead of directly applying contact constraints using the Matlab integrator ode45, ode113, and ode4. The authors of [16] conclude that the KEC method reduces computational cost up to 20% compared to the LUT method using Matlab implicit integrator ode15s. And in [17], the KEC method proved to be 12.6 times more efficient than the 3D elastic contact simulation with the Matlab implicit integrator ode15s. Expanding on this preceding research, the objective of this work is to study the numerical and computational aspects of both methods using the multibody model of the Manchester wagon 1 proposed in [3]. Since efficient and stable time integration is a key issue [23] for the computational multibody dynamics of railroad vehicles, different well-known numerical integrators have been studied and compared in this paper. Two fixed-step-size integrators, RK4 and a predictor-corrector based ABM, and two variable-step-size Matlab integrators, ode15s and ode45, are applied to solve railway dynamic problems. This will reveal numerical and computational aspects of the constraint-based contact approaches and will look at ways an accurate physics-based wheel/rail vehicle system could be analyzed in real-time.
The following Sect. 2 of this paper briefly introduces the multibody modeling of railway vehicles. Computer simulation of railway vehicles is presented in Sect. 3. Different integration methods are briefly introduced in Sect. 4. Sect. 5 compares the kinematics solutions of the modeled Manchester Wagon 1 with those obtained using standard commercial software and investigates the numerical and computational aspects of the KEC and LUT methods in detail. Finally, Sect. 6 provides a summary and offers conclusions.

Multibody modeling of railway vehicles
This section introduces details of the multibody modeling of railway vehicles. To this end, track kinematics vehicle kinematics and wheel/rail contact will be discussed. More insight into the kinematics and contact mechanism can be found in [12,16,24].

Track geometry and kinematics
In the multibody simulation of railway vehicles, the use of moving coordinate systems that follow the vehicles motion along the track centerline has been widely used in the literature [25,26]. This method has advantages in the interpretation of the results because system coordinates are referred to track geometry. This work uses the so-called 'track frame' O t ; X t , Y t , Z t for each vehicle body. The 'track frame is defined along the track centerline and has the X t axis parallel to the tangential line of the track centerline, the Y t axis connecting the centerline of the two rails and the Z t axis perpendicular to the plane of the rails. See Fig. 2.
The absolute position vector and rotation matrix from a track frame to the global frame can be expressed as the functions of the arc-length s as follows: where R t x , R t y , R t z are the components of the position vectors in the X, Y , and Z direction, ψ t (azimut or heading angle), θ t (vertical slope, positive when downwards in the forward direction), and ϕ t (cant or superelevation angle) are Euler angles which define the track centerline orientation.
With the description of the track frame, the geometry of the left and right rails centerlines can be obtained by defining two additional profile frames (left rail profile frame  O lrp ; X lrp , Y lrp , Z lrp and right rail profile frame O rrp ; X rrp , Y rrp , Z rrp ) at the railheads. Both rail profiles are separated a distance 2L r and rotated an angle β in the ideal track configuration, as it is shown in Fig. 3. In addition, this ideal geometry description allows a straightforward definition of the irregular track by means of the left and right irregularity vectorsr lir andr rir and the linearized irregularity angle δ defined as δ = (z lir − z rir )/2L r . See Fig. 3. Note that the components of the irregularity vectors can be extracted from the following well-known centerlines' irregularities measured in the railway industry: All of these components can be implemented in a track preprocessor subroutine to parameterize the three-dimensional (3D) ideal track centerline and track irregularities as continuous functions (cubic splines in this work) of the longitudinal arc-length trajectory coordinate s. A schematic representation of the procedure used for the track preprocessor is presented in Fig. 4. for 1000 m track. As for track irregularities, the track preprocessor is discretized with s = 0.1 m, which results in a Matlab storage size of around 300 kB (matrix size 10000 × 5) for 1000 m track.

Wheel/rail contact
As shown in Fig. 2, a body track frame O ti ; X ti , Y ti , Z ti can be associated with body i at any time-instant, when a body moves along the track with arc length s i . Therefore, the position and orientation of the body track frame can be computed according to Eq. (1), as: Similarly, a wheelset track frame O wti ; X wti , Y wti , Z wti is associated with the wheelset i and with time. See Fig. 5. The position and orientation of the wheelset track frame can be calculated by replacing s i with s wi in Eq. (2).
The set of relative wheelset track frame coordinates of a rigid wheelset i (superscript wi) is given by: where y wi is the lateral displacement and z wi is the vertical displacement of the wheelset frame with respect to the wheel track frame, respectively, and the X-component along the track centerline is zero; ϕ wi , θ wi , and ψ wi are Euler angles which define the orientation of relative wheel frame with respect to wheelset track frame.
The position vectors of right and left contact points, defined in the right and left railheads as shown in Fig. 3, are written with respect to wheelset track frame: where lr and rr represent left rail and right rail,r rrp andr lrp are position vectors of ideal left and right railhead frames with respect to wheelset track frame,r lir andr rir are rail irregularity vectors of rail profile frames with respect to ideal left and right rail profile frames, A t,lrp and A t,rrp are rotation matrices from the railhead frames to the wheelset track frame,û rrp rc andû lrp lc contain the components of the position vector of contact points in the rail profile frames as shown in Fig. 3. The vectors and matrices are written as The position vectors of contact pointsû rrp rc andû lrp lc in the rail profiles from Eq. (4) are parametrized according to the railhead profile geometry as shown in Fig. 6: where s lr 2 and s rr 2 are left and right transverse rail surface parameters, respectively, and h lr and h rr are the functions that define the railhead profiles.
Also, special treatment is given to the railway vehicle wheelset bodies in a railway vehicle. An additional frame called 'wheelset intermediate frame' (wI) is defined as shown in Fig. 5. The wI-frame is defined as a body frame but shows no pitch rotation. As will be shown in the following section, it is a convenient frame to define contact point position vectors, because they are not influenced by the wheelset pitch rotation, which is particularly high in these bodies. More details about body and wheelset kinematics can be found in [12,16,24]. The position vector of contact points on the surface of the left and right wheel profile in Fig. 5 can be obtained in wheelset track frame, such as r wi lc =r wi + A wti,wIiûwIi lc , where vectorr wi = 0 y wi z wi T is the position vector of the wheelset frame with respect to the wheelset track frame,û wIi lc andû wIi rc are the position vectors of the contact points at the left and right wheels with respect to the wheelset intermediate frame, and A wti,wIi is the where A wti,wIi is the rotation matrix of the wheelset intermediate frame with respect to the wheelset track frame, A wIi,wi is the rotation matrix of the wheelset frame with respect to the wheelset intermediate frame. These are given by The position vectors of contact pointsû wIi rc andû wIi lc in the wheel profiles from Eq. (7) are parametrized according to the wheel profile geometry as shown in Fig. 6: where s w 1 and s w 2 are the transverse and angular wheel surface parameters, respectively, and h lw and h rw are the functions that define the wheel profiles.
Based on the parameterization of the wheel and rail surfaces, contact can be modeled using the already defined elastic or constraint approach. For both approaches, 3D surfaceto-surface contact can be reduced to a planar two-dimensional (2D) curve-to-curve contact assuming that the contact points lie in the Y wti -Z wti plane of the wheelset track frame. This decreases the number of surface parameters from four to two, since the contact points at the left and right rails have the same arc-length parameter as the wheelset (s r 1 = s wi ), and the angular wheel parameters of the contact points are constant (s w 2 = π/2) in the wI-frame. Therefore, the simplified nomenclature s w = s w 1 and s r = s r 2 can be and is used in the rest of this paper. However, the main implication of this planar contact is that the wheelset yaw angle is neglected for the location of the contact points, which avoids the consideration of the so-called lead-lag contacts (flange contact points that are longitudinally displaced with respect to tread contacts). These contacts have an effect when the vehicle is negotiating very narrow curves, which usually occurs at very low velocities. Nevertheless, the planar contact approach is sufficiently accurate for most practical applications as shown in [12].
The use of 2D wheel/rail contact using the constraint approach is considered in the following. To this end, the following subsections briefly describe the basis of the following two planar constraint contact methods: contact lookup tables and the KEC method.

Creation of constraint contact lookup table
The procedure to build a constraint contact lookup table starts by solving the contact constraints for different values of the wheelset relative position with respect to the track. The algorithm is summarized as follows: 1. Two sets of discrete numerical values are assigned to the lateral displacement of the wheelset y wi and gauge variation ξ g : where y wi and ξ g are increments.
2. The wheelset coordinates s wi , θ wi , and ψ wi are set to zero (contact points are assumed to lie in the Y wti -Z wti plane of the wheelset track frame). 3. For a one-point contact scenario (when two wheel/tread contacts occur at both wheels), the six contact constraint equations of Eq. (12) are solved for different values of y wi and the nominal distance 2L * r = 2L r + ξ g in two for-loops to obtain six outputs: the wheelset kinematics z wi , ϕ wi , and the vector of tread surface parameters s, where rc represents right wheel and lc left wheel for one wheelset, s = [s lw s lr s rw s rr ] T is the vector of the surface parameters,t wi c is the unit tangential vector, andn rp c is the normal vector defined in the wheelset track frame, which can be expressed as follows: wheret wIi c is the unit tangential vector with respect to the wheel surface at the contact point andn rp c is the unit normal vector with respect to the railhead surface at the contact point [16,25]. In Eq. (12), given the value of lateral displacement y wi , six contact constraints for one wheelset are solved to get six outputs: the wheelset coordinates z wi and ϕ wi , and the surface parameters s of the contact points at the tread. 4 where clt refers contact lookup table. Figure 7(a) shows the wheelset vertical displacement coordinate z wi with respect to the track centerline within a range of track gauge variations ξ g and Fig. 7(b) illustrates the contact lookup table.

Contact detection for two-point contact scenario
For the two-point contact scenario, contact detection at the flange is addressed using the maximum relative-indentation condition [5]. Therefore, two nonlinear equations are solved to obtain two flange surface parameters s fla = [s w fla s r fla ] T as follows: Contact point detection at both the tread and flange for the LUT method is determined in a preprocessing stage that builds the contact lookup tables. During the dynamic simulations, the contact lookup table is used to interpolate the stored discrete data to get the contact point location.

Contact force computation for two-point contact scenario
For the LUT method, constraint contact is not appropriate when the number of contact points varies (i.e., from one contact point to the two-contact point scenario). For this reason, when flange contact occurs in the two-point contact scenario (before wheel climbing), a hybrid method is usually adopted [12]. In this approach, wheel/tread contact with the rail is considered using the constraint method and the wheel/flange contact with the rail is considered using the elastic method. The equations of motion for a railway vehicle subject to contact lookup table constraints can be written in the form of the following differentialalgebraic equation (DAE): where M is the vehicle generalized mass matrix, C LUT includes the LUT constraint equations C wi,LUT for the wheelset and forward velocity holonomic constraints, C LUT q is the corresponding Jacobian matrix, λ is the array of Lagrange multipliers, Q tang is the vector of generalized tangential tread and flange forces, Q susp is the vector of generalized suspension forces, and Q includes the vectors of generalized applied forces and generalized quadraticvelocity inertia forces. The wheelset LUT constraint equations contain two independent nonlinear equations using the lookup tables for vertical displacement and roll angle as presented in [12], whereȳ wi ,z wi , andφ wi are wheelset kinematics with respect to the irregular track centerline, which have the relationship with respect to those related to the ideal track centerline: The normal contact forces in the tread are computed as reaction forces associated with contact constraints from Eq. (16). Therefore, those forces are computed using Lagrange multipliers where C wi,LUT q is the LUT constraint Jacobian matrix associated with wheelsets. In addition, the normal contact forces for the flange are included in the right-hand side of the equations of motion as the applied force from Eq. (16) and it is computed using a Hertzian model based on interpenetration and interpenetration rate: where K hertz and D damp are Hertzian contact stiffness and damping coefficients, the Hertzian stiffness is computed based on the curvatures at the wheel/rail surface and material properties [25,27] under the assumption that the contacting bodies behave like infinite semispaces. The variables δ wi andδ wi from Eq. (20) are the wheel/rail penetrations at the flange and rate at the flange contact. These are given by When the two-point contact scenario occurs, the wheel-flange contact event suddenly appears and is treated as an impact. Therefore, the elastic flange contact event may require small time step size for the time integration [28] and the flange contact stiffness controls numerical simulation performance [24].

KEC method
The KEC method is a computationally efficient online constraint-based method, in which an equivalent wheel profile that contacts an infinitely narrow rail provides the same wheelset kinematics as the real wheel/rail profiles. See Fig. 8.

Construction of KEC-lookup table
To use the KEC method, the equivalent wheel profiles must be obtained first. To this end, the results of Eq. (14) are needed to relate the location of the contact points in the KEC and real profiles by building a KEC-lookup table [16]. This algorithm is summarized as follows: 1. The track irregularities y lir , z lir y rir , and z rir are set to zero. 2. For a one-point contact scenario, the six contact constraint equations of Eq. (12) are solved in one "for-loop" to obtain six outputs: the two wheelset coordinates and the four wheel surface parameters.
• When two wheel/tread contacts occur at both wheels, lateral displacement y wi is used as the independent coordinate to find the wheelset kinematic values z wi , ϕ wi , and the tread surface parameters s = [s lw s lr s rw s rr ] T .
where klt represents the KEC-lookup table. These functions are used to find the position of the contact points in the real profiles once the position in the KEC profile has been obtained. The surface parameter's s lw location of contact points in the real left-wheel profile is plotted as a function of the lateral position s lk of the contact point in the KEC profile in Fig. 9.
The computation of KEC-equivalent wheel profiles requires solving Eq. (22) including no track irregularity. The authors of [16] state that the area of allowable motion for the wheel with null irregularity cannot be guaranteed to be the same as for wheels with irregularities. Later, in [24], the KEC-equivalent wheel profiles were proved to be acceptable for use in irregular tracks with accuracy.

Contact detection for two-point contact scenario
As shown in Fig. 9, with a specific value of the equivalent parameter s lk f , two simultaneous contact points in the real profile s lw f and s lw t appear, which corresponds to a two-point contact scenario (before wheel climbing). However, this discontinuity in Fig. 9 might generate an unstable numerical integration problem. Therefore, the two-point contact scenario can be simulated by adopting a smoothed transition from tread to flange with the help of a regularization approach [17]. Using the transition lengths [17] of the regularization for the tread-flange transition to account for the two-point contact scenario is critical for simulation stability and accuracy. In this work, the transition lengths used are ŝ = 0.5 mm and s = 1 mm. If the transition lengths are too small, the equation of motion matrix will be illconditioned, and the simulation will stall. Alternatively, if the transition lengths are larger, the simulation results may be less accurate.

Contact force computation for two-point contact scenario
One advantage of the KEC method is that contact forces on the tread and flange can be treated equally as reaction forces. The equations of motion of a railway vehicle subject to KEC constraints can also be written in the following DAE form: where C KEC includes the KEC constraint equations for wheelset C wi,KEC and forward velocity holonomic constraints, the matrix N includes the holonomic constraint Jacobian matrix and the matrix N wi , which represents the direction of the reaction forces and is associated with all the wheelsets in the vehicle. The wheelset KEC constraint equations C wi,KEC in the dynamic simulation are the same as in Eq. (22). In this case, lateral displacement is treated as an independent coordinate, and the outputs of the nonlinear constraint equations from Eq. (22) become the vector of wheelset dependent coordinates x wi = [z wi ϕ wi s lk s rk ] T [16].
The normal contact forces for the wheels are computed as the reaction forces as follows: where N wi represents the direction of the reaction forces and torques and is associated with the wheelset. Due to the nonlinearity, the lateral position vector s k cannot be eliminated from Eq. (22). Therefore, the reaction forces for the wheelset Q nor,wi cannot be obtained using the KEC constraint Jacobian matrix Q nor,wi = −C wi,KEC q λ in Eq. (19) [17,24]. As shown in the left and right of Fig. 10, points A and B represent the wheel/tread and wheel/flange contacts when a two-point contact scenario occurs. The transition length between points A and B in the left of Fig. 10 corresponds to the forbidden area between points A and B on the left wheel profile in the right. However, in the numerical simulation, the lateral position of point C in the transition length s lk of the KEC-equivalent profile (see left of Fig. 10) might be obtained by solving Eq. (22). In this case, the contact forces will be directly applied at the contact point C in the forbidden area in the right of Fig. 10. Using contact force Q nor,wi at point C directly, the vehicle dynamic behavior might not be accurate. Therefore, the normal contact force Q nor,wi is transformed into two normal contact forces at the tread-end point Q nor,wi

Computation of tangential contact forces
In the dynamic simulation of the rail vehicle, the tangential forces and spin moment are generated between wheels and rails due to the relative motion of rolling and sliding [25].
Because it offers good computational efficiency and accuracy, the tangential contact force here was computed using Polach creep contact theory [29,30] and treated as applied forces in the equations of motion Eqs. (16) and (24). In general, the computation of tangential contact forces requires the inputs of: (a) normal contact forces acting on the contact points, (b) relative contact velocity, (c) Kalker's constants, and (d) coefficients of friction. In this work, the coefficient of friction was considered constant. The relative contact velocity and Kalker's constants were computed based on the rail and wheel surface curvatures, which can be stored in a lookup table. The normal contact forces used here depends on two different situations: 1. If the normal contact force is computed as the reaction forces, such as the normal tread contact force with LUT method (see Eq. (19)) and normal tread or/and flange contact force with KEC method (see Eq. (25)), the normal forces obtained last time-step is used to find tangential contact force. This assumption can work is simply due to the small time-steps used in railway dynamic simulation. 2. If the normal contact force is computed as the applied forces, such as the normal flange contact force with LUT method (see Eq. (20)), current time-step normal contact forces are implemented for tangential contact force computation.

Computation of suspension forces
This work considers that each modeled body is accompanied by a track frame along the track centerline such that the body coordinates are defined with respect to this frame. See Fig. 11. However, when computing suspension forces, it is more convenient to employ a master track frame. Therefore, the relative kinematics for all moving bodies can be projected on the same track frame. This way, the absolute position vector of a point P and Q in the global frame (see Fig. 11) are given by: Fig. 11 Kinematics projected to master track frame: R ijt is the relative distance vector from body track frames it to jt,û i P and u j Q are the position vectors of the points P and Q with respect to their body frames,r i P andr i Q are the position vectors of point P and Q, which is projected into master track frame i where R ijt = R btj − R bti is the relative distance vector from body track frames it to jt, vectors R bti and R btj are the position vectors of body track frames with respect to global frame, matrices A bti and A btj are rotation matrices of body track frames,r i andr j are position vectors of body i and j with their track frames,û i P andû j Q are the position vectors of the points P and Q with respect to their body frames, and matrices A bti,i and A btj,j are the rotation matrices of body frames with respect to the body track frames. If the position vectors of points P and Q are referenced with respect to the master track frame with their components expressed with respect to this frame, Eq. (27) can be expressed as Therefore, the suspension force of the spring-damper in parallel in Fig. 11 is computed as follows: where q i is the set of generalized coordinates for body i, K is the spring stiffness, D is the damping coefficient, and the deformation length of the spring d is defined as: where l 0 is the undeformed length of the spring. The time derivativeḋ is expressed as [31] d = ∂d ∂q iq i .
The calculations of suspension forces with different suspension elements in the Manchester wagon 1 are given in Appendix A.

Computer simulation
The dynamic simulation of a railway vehicle integrates Eq. (16) or (24) forward in time. This procedure, developed in the framework of multibody dynamics, is described here as a function of the two wheel/rail contact methods used in this work and illustrated in Fig. 12 as follows: 1. Set initial conditions for the system generalized coordinates and velocities. 2. If using the contact LUT method, interpolate stored discrete data from the contact lookup table with two entries,ȳ wi = y wi − ξ a and ξ g , and obtain the wheelset coordinatesz wi = z wi + ξ v ,φ wi = ϕ wi + ξ c /2L r , and surface parameters s lw , s lr , s rw , s rr . Normal contact forces in the tread are computed as reaction forces associated with contact constraints, while in this work, the normal contact forces for the flange are computed as a Hertzian penetration-based model. 3. If using the KEC method, solve the KEC constraint equations from Eq. (22) to get the wheelset coordinates z wi , ϕ wi , and KEC-profile positions s lk , s rk with inputs y wi and track irregularities y lir , y rir , z lir , z rir . Interpolate stored discrete data from the KEC-lookup table to get surface parameters s lw and s lr with the entry of s lk and s rw and s rr with the entry of s rk . Normal contact forces on the tread and the flange are treated equally as the reaction forces. The two-point contact scenario can be simulated with a smooth transition of the normal contact forces from tread to flange. 4. Include the contact force vectors at each wheel, suspension force vectors (primary and secondary), Newton-Euler generalized forces and quadratic-velocity generalized inertia forces into the vector of external forces of the system. Solve the equations of motion of the railway vehicle to get the generalized accelerationsq. 5. Return the new generalized positions q and velocitiesq at time t + t with the numerical integration subroutine. 6. Go to step 2 or 3 until the analysis is finished.

Integration methods
The numerical integration method can be divided into explicit and implicit integration. The decision to use an explicit or implicit method depends on the studied system. In the explicit method, the calculation of the variables at each time step requires the solution of the linear equations that are functions of the values of the variables in previous time steps. In the implicit method, the calculation of the variables at each time step requires the solution of the nonlinear equations that are functions of the values of the variables in the current and the previous time steps. Therefore, the implicit method uses iteration to solve nonlinear equations at each time step.
This study uses a fourth order Runge-Kutta method as the fixed-step-size explicit integrator [32]. It is a single step method, and four function evaluations are required in each integration step. As an alternative time-integration scheme, the predictor-corrector ABM method is used as the fixed-step-size in this work. It allows determining the errors in each iteration step (local truncation error) and a correction term can be included. The Adams-Bashforth-Moulton scheme can be written as follows: (Predictor) y n+1 = y n + t (55f(t n , y n ) − 59f(t n−1 , y n−1 ) + 37f(t n−2 , y n−2 ) − 9f(t n−3 , y n−3 ))/24, (Corrector) y n+1 = y n + t (9f(t n+1 , y n+1 + 19f(t n , y n ) − 5f(t n−1 , y n−1 ))/24, (32) where t n−0..3 and y n−0...3 are the three previous time steps and state variables.
The ABM is a multistep method that has two function evaluations in each integration step, the predictor and the corrector. Compared to the RK4 with four function evaluations, the ABM should be twice as fast.
The Matlab built-in function ode45 is typically selected as the variable-step-size explicit integrator, as it can provide an acceptable solution in most cases. It is a single-step explicit integrator based on the Runge-Kutta (4,5) formula (Dormand-Prince method). For stiff systems, the solution time with ode45 increases, and other Matlab built-in time integrators schemes such as ode15s can perform better; ode15s is a variable order and variable-stepsize implicit solver based on numerical differentiation formulas (NDFs) of orders 1 to 5. Optionally, the maximum order (5) can be reduced or a backward difference formulation can be used [33]. The railway vehicle models usually contain stiff springs for suspension components or contact forces modeling. The explicit integration methods, such as RK4, are not suitable for a stiff problem. In many cases, the explicit integration method is conditionally stable and it need least effort for one time integration step. However, the method requires a small time-step size to solve the stiff problem and makes the simulation less efficient [34]. Furthermore, the predictor-corrector ABM method includes the previous step information to prevent the numerical instability. The Matlab ode45 is based on Runge-Kutta algorithm and has similar features as the RK4 for a stiff problem. The Matlab ode15s is based on NDF formulation, where the time-step and solution order can be changed. These features are desired in the stiff problem solution, as dynamically adjusting these parameters allows obtaining the solution efficiently. However, it requires the iterative solution for the nonlinear equations and the update of the Jacobian matrix at each time step. As the number of unknowns gets larger, the numerical effort for the update of Jacobian matrices starts to be more expensive [23].

Vehicle model description
This section describes the vehicle model used in the numerical simulations of Sect. 3, which is based on the Manchester wagon used in Iwnick [3] as a benchmark. Figure 13 illustrates the various components of the vehicle formed by one carbody, two bogie frames and four wheelsets with primary and secondary suspension elements.
The mass and inertia properties of the different bodies are presented in Table 2. The wheel and rail profiles used for the contact lookup table and computation of the KEC profiles are the S1002 wheel and LB140-Area rail profiles shown in Fig. 14, which present a unique two-point contact scenario in the tread-flange transition. Wheel and rail profile positioning with respect to the track centerline (L w , L r ), rail inclination (β) and wheel nominal radius (r 0 ) are given in Table 1. The coefficient of friction for the wheel rail tangential contact force is μ = 0.2364 and the Hertzian stiffness and damping constant for flange contact with the LUT method are listed in Table 3.
The primary and secondary suspension shown in Fig. 13 involves a total of 48 elements classified into springs and dampers in series, spring and dampers in parallel, and bumpstops, with their stiffness and damping properties given in Table 4.  Damping D damp (N·s/ m 2 ) 7 .075 · 10 11 1 · 10 10 1 · 10 9 1 · 10 8 1 · 10 7 Track irregularities are generated using analytical expressions of the power spectral density functions (PSD) [35]. The alignment, vertical profile, gauge variation, and cross level are shown in Fig. 15.

Validation
In this section, the 3D multibody model of the railway vehicle is implemented in commercial simulation software, Universal Mechanism (UM) [36] to validate that implemented in Matlab software. For both the KEC and LUT methods in Matlab software and UM commer-

Single wheelset on a curved track without track irregularities
In the first case, the single wheelset is simulated at a constant forward velocity of V = 10 m/s (36 km/h) on a 120 m track without irregularities constructed of the following three  Fig. 16. Figure 17 shows the comparison of the lateral displacement and yaw angle of the wheelset with the KEC and LUT methods and UM. The achievement of the steady curving can be observed when the vehicle enters into transition curve (see grey area). In addition, with different values of curve radius, the traveled distance where steady curving starts is different. The steady curving results are quantified in Table 5. The results of lateral displacement and yaw angle show good agreement between both of the examined methods and UM, except for the yaw angle, which is about 2 mrad bigger when using UM on the radius R = 300 m curve track. That is mainly because the planar contact approach is implemented into the KEC and LUT methods, but in UM, the 3D elastic contact approach is utilized for the tread and flange contacts.

Manchester wagon on a curved track without track irregularities
In the second case, the Manchester wagon is simulated at a constant forward velocity of V = 20 m/s (72 km/h), on a 500 m track without irregularities formed by the following three segments: a 100 m tangent, a 50 m transition, and a 350 m left curve of R = 500 m radius segment (see Fig. 18).   Figures 19 and 20 show a comparison of lateral displacements and yaw angles for the vehicle wheelsets, bogies, and car body, when using the KEC and LUT methods and UM. Table 6 lists the quantities of the steady curving results. Again, the results of lateral displacement and yaw angle show good agreement. A steady curving motion is achieved when the vehicle enters into the transition curve. According to the lateral displacement of the four wheelsets, the front wheelsets of both bogies can experience flange contact when negotiating the curve (blue area in the figure), and the rear wheelsets do not.

Manchester wagon on a curved track with track irregularities
In the last case, the Manchester wagon is simulated at a constant forward velocity of V = 20 m/s (72 km/h) on a 500 m track (see Fig. 18) with irregularities. Figures 21 and 22 show a comparison of the lateral displacement and yaw angle for the vehicle wheelsets, bogies, and car body when using the KEC and LUT methods and UM. The results of lateral displacement and yaw angle are almost identical with the results on the irregular track. According to the lateral displacement of the two wheelsets, steady curving motion is not achieved because of the track irregularities.

Computational efficiency
In this case, the Manchester wagon is simulated at a constant forward velocity of V = 20 m/s (72 km/h), on the same 500 m track as shown in Fig. 18 with irregularities. The computational efficiency of the LUT method is compared to that of the KEC method using different integration methods. The comparison was carried out on an Intel Core i5 laptop with a 2.5 GHz CPU and Matlab 2018b. The Mex functions written in C or Fortran were not used to improve computational efficiency in Matlab. To determine the computational time, each case is simulated five times, and the calculated simulation time was as the average value of the five runs. The computation times for the five runs for each case are listed in Appendix B. Two fixed-step-size integrators, explicit RK4 and predictor-corrector ABM, and two variable-step-size Matlab built-in function integrators, explicit ode45 and implicit ode15s, were applied to compare computational efficiency.
According to the Hertz contact theory [25], the flange Hertzian stiffness for the wheelrail profile combination shown in Fig. 14 is K hertz = 7.075 · 10 13 N/m 1.5 with the Poisson's ratio ν = 0.28 and Young's modulus E = 2.1 · 10 11 N/m 2 . Since the location of the flange contact points remains the same when the two-point contact scenario occurs, the constant flange contact stiffness parameters are considered for the hybrid method in the con-  Table 7. With step size t = 1 ms, RK4 and ABM failed to converge at same time as when flange occurs. In addition, the variable-step-size explicit integrator ode45 only succeeds when the absolute and relative tolerance is smaller than 1 · 10 −8 and the variablestep-size implicit integrator ode15s only succeeds when the absolute and relative tolerance is smaller than 1 · 10 −7 . This is because having the high-magnitude Hertzian parameter, K hertz = 7.075 · 10 13 N/m 1.5 , results in a stiff system of ordinary differential equations, and some variations of the equation significantly impact the solution. If the time step and/or tolerances are too large, the first impact between the wheel flange and rail head may result in a substantial penetration depth, which leads to a numerically instability and failure of the time integration. Therefore, a large number of time steps are needed to solve the problem, which results in high computational cost. Under the assumption that the contacting bodies behave like infinite semispaces, the Hertzian stiffness of K hertz = 7.075 · 10 13 N/m 1.5 might be higher than the real value. In fact, reducing the value of Hertzian stiffness is realistic because the influence of flange flexibility is not considered. In what follows, different values of flange contact stiffness will be presented. Computational cost and function evaluations of the vehicle simulation with different values of flange Hertzian stiffness for LUT method are compared against KEC method in Fig. 23. The computation time of five runs for each case are listed in Tables 12,13,14,15,16. When using variable-step-size integrators, ode15s and ode45, the KEC method is much superior to the LUT method with higher flange contact stiffness (such as K hertz = 1 · 10 11 N/m 1.5 and 1 · 10 12 N/m 1.5 ). This is due to the large value of stiffness that leads to a stiff system of ordinary differential equations and increases the computational effort during the simulation. It is true that stiff integrators (ode15s) are more expensive. But in case of stiff systems, they lead to a stable solution that is not sometimes possible or take even more time in the case of nonstiff integrators (ABM and RK4). When working with lower flange stiffness K hertz = 1 · 10 10 N/m 1.5 and 1 · 10 9 N/m 1.5 , the LUT method shows better performance in computational efficiency. Table 8 lists a comparison of maximum and average errors of lateral displacement between LUT and KEC methods when using different integrators, with the Manchester wagon simulated at a constant forward velocity of V = 20 m/s (72 km/h) on a 500 m track with irregularities. See Fig. 18. Different values of flange stiffness are considered at the wheel flange. The fixed step sizes are t = 1 ms for fixed-step-size integrators RK4 and ABM. The maximum step size is t = 1 ms for Matlab integrators ode45 and ode15s. The absolute and relative tolerance for ode45 and ode15s is 1 · 10 −6 . The error between both methods is computed as

Accuracy
where y wh KEC is the lateral displacement for the KEC method within a 25 s simulation, and y wh LUT is for the LUT method. Some conclusions can be made based on the numerical results from Table 8: • For all integrators, the lower flange stiffness for the LUT method (especially K hertz = 1 · 10 9 N/m 1.5 ) results in higher maximum and average absolute errors. This is due to the low flange stiffness that results in large flange to rail head indentations. • Compared to the other integrators, the ABM errors are much higher when selecting a large value of the flange contact stiffness K hertz = 1 · 10 12 N/m 1.5 . This is because constant timestep integrators are not appropriate for impact problems. As shown in Fig. 24, the flange indentation and normal contact forces become artificially large for the ABM integrator.
As shown in Fig. 25, the use of a low stiffness (K hertz = 1 · 10 9 N/m 1.5 ) helps to simulate the vehicle motion. However, low stiffness results in high indentations (around 0.75 mm) that can be considered as physically inadmissible. Even though the flange indentations are different for different flange contact stiffnesses K hertz in Fig. 25, the resulting flange normal contact forces F hertz are on average almost the same. However, the higher the stiffness, the higher the force oscillations.  24 Comparison of lateral displacement, flange indentation, and normal contact force with the LUT method when using ABM, RK4, ode45 and ode15s. The fixed step sizes are t = 1 ms for fixed-stepsize integrators RK4 and ABM. The maximum step size is t = 1 ms for the Matlab integrators ode45 and ode15s. The absolute and relative tolerance for ode45 and ode15s is 1 · 10 −6 . Flange contact stiffness K hertz = 1 · 10 12 N/m 1.5

Fig. 25
Comparison of flange indentation and normal contact force with LUT method when using ABM with step size t = 0.1 ms. Different flange contact stiffnesses are considered a set of nodal points. Later in the dynamic simulations, the required parameters or variables for the dynamic analysis are obtained by using linear interpolation between the stored data. The KEC table has only one entry, the lateral wheel parameter s lk or s rk . In addition, the LUT table has two entries, gauge variation irregularity ξ g and wheelset lateral displacement y w . In this section, each KEC table is refined by increasing the number of  points from 465 × 1 to 5001 × 1 as listed in Table 9. Furthermore, each LUT table is refined by increasing the number of points from 361 × 11 to 721 × 11 as listed in Table 10.
The refinement of lookup tables tends to increase the precalculation time and the computer storage while it is assumed to improve the accuracy of the interpolated values [15]. To verify if the refinement of lookup tables can improve the computational efficiency, Fig. 26 compares the computation cost by using KEC and LUT methods with the flange Hertzian stiffness 1 · 10 10 N/m 1.5 . The computation time of five runs for each refinement case are listed in Tables 17-18. The results indicate that the refinement of both KEC and LUT tables can slightly improve the computational efficiency in most particular cases (14 out of 18). Table 11 lists the computational costs and function evaluations of the vehicle simulation with flange Hertzian stiffness K hertz = 7.075 · 10 13 N/m 1.5 for the LUT method after refinement. Compared to the result before refinement from Table 7, it can be concluded that the refinement only helps ode15s and ode45 to smooth the numerical integration that increases the maximum tolerance to 10 times bigger.

Conclusion
Two constraint-based formulations for the wheel/rail contact simulation were studied in terms of their accuracy and computational efficiency based on the Manchester wagon 1 model. Two fixed-step-size integrators, RK4 and predictor-corrector-based ABM, and two  variable-step-size Matlab integrators, ode15s and ode45, were used to carry out the comparison.
The studied contact models are the KEC and LUT methods. The differences between each approach can be summarized as follows: 1. Contact detection for the LUT method is based on an offline approach where the contact points are determined by interpolating a predefined lookup table. The KEC method is an online contact detection approach where the exact position of the contact points can be determined from the online solution of the KEC constraints. 2. The second difference is the computation of normal contact forces at the flange. For the LUT method, the normal contact forces are defined in terms of stiffness and damping ratio using Hertz's contact theory. This implies that computational efficiency and numerical stability will be affected by the large magnitude values of the stiffness parameters. However, with using the KEC method, the normal contact forces are determined through Lagrange multipliers that avoid the numerical disadvantages in the case of stability using high values of contact stiffnesses.
Numerical simulations of vehicle dynamics have been carried out in different case studies with and without track irregularities. Comparison between the kinematic results obtained by Matlab and the commercial simulation software UM has allowed a reliable validation of both the LUT and KEC methods. The performance of both approaches is compared in terms of efficiency and accuracy. The KEC method is faster (more computationally efficient) than the LUT method when using all the studied integrators (up to 20%). Also, as expected for the LUT method, results are more accurate when Hertzian stiffness is lower. The integrators ode15s and ode45 showed better efficiency compared to the integrators RK4 and ABM, only when flange contact stiffness K hertz is lower than 1 · 10 10 N/m 1.5 . As for the KEC method, the integrators ode15s and ode45 are faster than the fixed-step-size case under some specific tolerances.