Wheel–rail contact model for railway vehicle–structure interaction applications: development and validation

An enhancement in the wheel–rail contact model used in a nonlinear vehicle–structure interaction (VSI) methodology for railway applications is presented, in which the detection of the contact points between wheel and rail in the concave region of the thread–flange transition is implemented in a simplified way. After presenting the enhanced formulation, the model is validated with two numerical applications (namely, the Manchester Benchmarks and a hunting stability problem of a suspended wheelset), and one experimental test performed in a test rig from the Railway Technical Research Institute (RTRI) in Japan. Given its finite element (FE) nature, and contrary to most of the vehicle multibody dynamic commercial software that cannot account for the infrastructure flexibility, the proposed VSI model can be easily used in the study of train–bridge systems with any degree of complexity. The validation presented in this work proves the accuracy of the proposed model, making it a suitable tool for dealing with different railway dynamic applications, such as the study of bridge dynamics, train running safety under different scenarios (namely, earthquakes and crosswinds, among others), and passenger riding comfort.


Introduction
The dynamic effect on structures, such as bridges, caused by moving vehicles is a topic that has been attracting researchers and engineering practitioners for a long time. These effects can be assessed through transient moving load models [1][2][3] or by using more realistic vehicle-structure interaction (VSI) models [4][5][6][7]. While the first method is restricted to the analysis of the structural response, since the vehicle is represented as a set of moving loads corresponding to its static axle loads, the latter can also be used to assess the vehicle's behavior, including its dynamic response and the contact forces that arise from the contact interface. These models can be used in both roadway [8][9][10] and railway [4][5][6][7] applications, with differences in the contact interface, namely, between tire-surface and wheel-rail contact mechanisms. Since the present work focuses only on railway, only the latter will be addressed hereinafter.
VSI models in railways can be used in a wide range of applications, such as the study of bridge dynamics [11][12][13][14], evaluation of the train running safety [15][16][17], and assessment of passenger riding comfort [18,19], among others. Depending on the objective, the VSI models may have different degrees of sophistication. The simpler ones are those that deal only with the vertical dynamic behavior, which have been widely developed since the 1980s [20][21][22] due to its simplicity on dealing with the coupling between the wheel and rails. However, they are restricted to the analysis of the vertical vibrations on both the train and the structure, thus not being able to deal with phenomena caused by lateral excitations, such as wind, earthquakes, or even lateral track deviations. Naturally, such drawback makes these models incompatible with any type of analysis that strongly depends on this type of excitations, namely, running safety analyses or riding comfort evaluation.
To overcome the aforementioned limitations, several researchers start to focus their work on developing complex contact models that could accurately simulate the behavior 1 3 Railway Engineering Science (2023) 31 (3):181-206 on the wheel-rail contact interface, thus allowing the evaluation of the contact forces in all the directions, including the lateral one. Within this topic, special attention has been given by mechanical engineers in the development of wheel-rail contact models in multibody environments, where the track flexibility is not considered or is considered in simplified forms. Examples of such applications can be found in several publications. Shabana and Rathod [23] made a comparison between planar and spatial contact detection methods using nonlinear geometric equations, in which the contact points could be in a 2D plane or a 3D curve (contact points in different longitudinal positions), respectively. The authors concluded that the two options do not affect the vehicle's critical speed, but some deviations may occur in the estimation of the contact forces due to small differences in the contact point locations. Sugiyama and Suda [24] proposed a hybrid strategy to detect the contact points, in which an offline approach is used to detect the tread contact point, and an online method is applied to detect the contact point in the flange. While the former is based on precalculated lookup tables that are interpolated during the analysis, the later consists of detecting the contact point in real time during the dynamic analysis. Bozzone et al. [25] proposed an alternative to the nonlinear geometric equations used by the previous authors. In this work, the position of the contact point is also determined online, but based on intersection volumes, i.e., the contact point location is determined based on the farthest points belonging to the intersection of the contacting surfaces that define the wheel and rail profiles. More recently, Magalhaes et al. [26] proposed a complex and comprehensive online contact model able to detect multiple contact points in switches and crossings scenarios, while Sun et al. [27] used an alternative offline contact approach based on lookup tables to efficiently detect contact points in hunting stability scenarios. These models, however, are restricted to multibody simulations and cannot address the response of the structure. Multibody software, such as SIM-PACK® [28] or VAMPIRE® [29], also use sophisticated wheel-rail contact models but they are also unable to model the flexibility of a structure, such as a bridge. Thus, these models are mainly used to study the behavior of the vehicle itself and to perform detailed analysis of the contact interface in switches and crossings scenarios or to study wear of the wheel-rail interface. Enhancements in this regard have been made recently by Antunes et al. [30], which developed a co-simulation process to interconnect the multibody capabilities of an in-house railway dynamics software [31] with a track model developed within a finite element (FE) framework. However, these important developments continue to be mainly focused on the behavior of the vehicle from the mechanical engineering point of view, thus neglecting the detailed analysis of the infrastructure dynamics or the modeling of more complex structures, such as bridges.
To study the coupling behavior of the train-bridge systems, several authors brought together the knowledge of mechanical and civil engineers, to develop tools that perform the dynamic train-bridge coupling with appropriate wheel-rail contact models and consider the full flexibility of the infrastructure through FE models with any degree of complexity. A detailed state-of-the-art review has been published recently by Zhai et al. [32], where the authors demonstrate the evolution of these models that took place in the last years, as well as their main applications. Some VSI models based on FE, however, do not take into consideration the track model [19]. By considering the flexibility of the bridge, but neglecting that of the track, the computation of the wheel-rail contact forces, which are crucial for the assessment of the train running safety based on derailment criteria, becomes less accurate. Therefore, several authors started to incorporate the track in the bridge models through elaborate flexible components that simulate the flexibility of the different components of the track, namely, ballast, pads, fasteners, etc. Zhai et al. [6,33,34] presented a theoretical framework of a train-track-bridge interaction tool, as well as several experiments conducted with it. Montenegro et al. [35] also developed the tool "VSI-Vehicle-Structure Interaction Analysis," which can model structures with any degree of complexity and that has been used in several train-bridge applications, such as running safety [17,36,37], evaluation of the train stability depending on the track condition [38], as well as in the evaluation of riding comfort considering different indicators in different directions [18]. Recently, other research groups also developed their own VSI tools to study the train-bridge coupling behavior under different conditions [15,16,[39][40][41][42], as stated in the literature reviews published by Montenegro et al. [43] and Li et al. [44]. However, most of these models that can be found in the literature do not thoroughly present the validation process undertaken by their developers, making it hard to check which examples and tests were used to guarantee the correct correspondence between the model performance and reality.
The present work aims to show not only an enhancement in the wheel-rail contact model presented before in [35], in which the wheel profile can now be parameterized without any simplification in the concave regions, but also to demonstrate its validity with three different validation examples, namely, the comparison of the results obtained with the proposed model with those published in the wheel-rail contact Manchester Benchmarks [45], the simulation of a hunting stability problem with semi-analytical solution, and the validation with an experimental test that took place in the Railway Technical Research Institute in Japan. The model proposed here has been implemented in a [46] framework, which couples the vehicle and structure subsystems modeled with an FE package, thus not restricting it to multibody analyses of the vehicle, but also to the analysis of the railway infrastructure, such as bridges or plain track, modeled with any degree of complexity.

Wheel-rail contact formulation
In a previous publication from the author [35], the contact model did not contemplate the possibility of finding wheel-rail contact pairs in the concave region of the wheel (transition zone between the tread and flange) due to the multiple solutions that could arise from the geometrical nonlinear system of equations that will be presented later in this section. In this work, this drawback has been overcome with an alternative algorithm for detecting contact pairs in the concave regions of the wheel, thus avoiding the simplification of the wheel surface geometry adopted in [35]. The full contact algorithm is described in the following sub-sections.

Geometric contact problem
The present wheel-rail contact formulation is based on profile surfaces that are parameterized using cubic splines defined from a group of points that represent, as rigorous as possible, the geometry of the contact surfaces. The wheel is parameterized by two functions (see Fig. 1), one for the tread f w,t and another for the flange f w,f , making the detection of the contact points in the two regions of the wheel profile fully independent (see [35] for details about the profiles parameterization) and allowing the detection of two contact points simultaneously, as shown in Fig. 1b. This option is particularly important in scenarios with important lateral loads, such as in the presence of earthquakes [37] or crosswinds [36], where double contact with the tread flange contact occurs often. In Fig. 1, the contact point and the point with maximum curvature, which defines the division between the two surfaces, are denoted by P c and C m , respectively, while f r is referred to the rail profile function.
After defining the surfaces of the contacting bodies, the next step of the geometric problem consists of determining the position of the contact points between the wheel and the rail. In the present work, and contrary to the author's previous work [35], two algorithms for the detection of contact points are implemented. The first algorithm is used to detect the position of contact points lying on convex regions of the surfaces (see Sect. 2.1.1), while the second is applied when the contact point is located on concave regions (see Sect. 2.1.2), as illustrated in Fig. 2.

Convex contact search
The location of the potential wheel-rail contact points is determined based on the following system of nonlinear equations: where t t r,y and t t w,y are the tangent vectors to the rail and wheel surfaces at the potential contact point, respectively; n t r is the normal vector to the rail surface at the potential  Transition zone contact contact point; d t wr is the vector that defines the relative position of the point of the wheel with respect to the point of the rail; and the superscript "t" refers to the fact that all the aforementioned vectors are defined with respect to the target element coordinate system (x t , y t , z t ). The x t axis has the direction of the longitudinal axis of the rail profile, the y t axis is parallel to the track plane, and the z t axis completes the right-handed system (see [35] for details). By simultaneously guaranteeing that the tangent vector to the rail is perpendicular to the vector that defines the relative position between the contact pair (first condition) and that the normal vector to the rail is perpendicular to the tangent vector of the wheel (second condition), it is possible to define the contact point position.
The system of Eq. (1) may have multiple solutions if one of the contact surfaces is not convex, i.e., if the potential contact point lies in the transition zone between the flange and the rail. Therefore, after solving Eq. (1), the algorithm checks the convexity sign of the wheel surface at the potential contact point location by computing its curvature w,y along the lateral direction y. According to Garg and Dukkipati [47], the radius of curvature of a surface is considered to be positive if the corresponding center of curvature is within the body, i.e., if the surface is convex. Thus, the potential contact point lies on a convex region if the following condition is fulfilled: otherwise, the potential contact point lies on a concave region, and the solution obtained with Eq. (1) is discarded. When this situation occurs, the concave contact search algorithm, presented in Sect. 2.1.2, is called to determine the actual position of the contact point.
In the scenarios in which the algorithm finds a unique solution for Eq. (1), i.e., the potential contact pair lies in a convex region, the existence of an actual contact point (2) w,y > 0; is not yet mathematical guaranteed. Hence, the two bodies are actually in contact only if the vectors d t wr and n t r point in opposite directions, which can be mathematically described through the following equation: Finally, the penetration d between the two bodies in contact is given by

Concave contact search
As mentioned before, one of the main contributions of the present paper consists of the development of an algorithm to detect the contact point position not only in convex regions, as in [35], but also in concave zones. The concave contact approach consists of determining the location of the contact points in the regions where the convex contact approach cannot find a single solution, i.e., in concave surfaces (see condition expressed in Eq. (1)). Unlike the algorithm used in the convex contact approach, the accuracy of this algorithm depends on the degree of discretization of the profiles. Therefore, although this approach may also be used to locate the contact points in convex regions, the higher computational cost that is required to achieve a good solution makes it computationally less attractive. As a result, the concave contact approach is used only if the convex approach finds a solution that lies in a concave region.
In the concave contact search approach, the rail and wheel surfaces are discretized in n r and n w points, respectively, by interpolating the profile functions f r and f w shown in Fig. 1. Hence, the evaluation of the potential contact between the two surfaces consists of determining if any of these points lie inside the opposite surface, forming the so-called intersection volume, illustrated in Fig. 3.
To determine which points belong to the intersection volume, the points belonging to the rail surface are projected into the wheel surface and vice versa. Then, the vertical distances between the points of each surface and the respective projection on the other surface, h r,i and h w,j , are computed as follows: where u t r,i and u t w,j are the position vectors of the ith rail point into the wheel surface and of the jth wheel point into the rail surface (see Fig. 3), respectively, defined with respect to the target element coordinate system; u If there are no points belonging to the intersection volume, the bodies are not in contact, and the potential contact point is discarded from further considerations. On the other hand, when contact is detected, each point of one of the surfaces belonging to the intersection volume has a potential contact pair in the other surface. Thus, the potential contact pair of a given point of the rail surface belonging to the intersection volume is the closest point of the wheel surface, which also belongs to the intersection volume, and vice versa. The distance d i between the ith rail point belonging to the intersection volume and the point of the wheel surface j that forms the potential contact pair is given by where n IV r and n IV w are, respectively, the number of points of the rail and wheel surfaces which belong to the intersection volume.
Finally, out of all the pairs giving the maximum distance between the rail point and the correspondent wheel , i = 1, 2, … , n IV r and j = 1, 2, … , n IV w , point, the pair where contact occurs is the one that leads to the maximum penetration d, given by A schematic representation of the selection of the contact pair ij is depicted in Fig. 4.

Normal contact problem
When two non-conforming bodies are compressed against each other, they will deform around the first point of contact and form a contact area. In the present model, the nonlinear Hertz contact theory [48] is used to analyze the normal contact problem. The normal contact force F n between the wheel and rail is given by where d is the penetration defined in Eq. (4) or (10), and K h is a generalized stiffness coefficient that depends on the material properties of the bodies in contact, such as the Young modulus and the Poisson ratio, and on the curvatures of the surfaces at the contact point [49]. Note that, in a wheel-rail contact problem, the assumptions of the Hertz theory are not met, since the surfaces of the contacting bodies are not totally frictionless and may be conforming. Moreover, the wheel and rail profiles may have non-constant curvatures in the contact area, and plastic deformations may occur in the contact zone. Hence, more complex and realistic contact shapes [50,51] may be necessary for analyzing local problems, such as wear or switch and crossings problems. Nevertheless, the Hertz theory is still widely used by  [52] due to its sufficient accuracy in predicting the wheel-rail normal contact forces in the majority of the railway applications, even when the Hertzian assumptions are not fully met. For this reason, the nonlinear Hertz theory is used in the present application in both the convex and concave contact regions. For the concave regions, the curvature ratio B is limited to a minimal positive value in order to allow the calculation of the Hertzian ellipse semiaxes, as proposed in [50]. In this work, this minimal value adopted was B = 1 × 10 -6 m −1 to ensure that this parameter is close to zero but does not pose any numerical problem to the algorithm.

Tangential contact problem
When two bodies are allowed to roll over each other, some points on the contact area may slip while others may adhere. The difference between the tangential strains of the bodies in the adhesion area leads to a small apparent slip, called creep, which is crucial for the determination of the tangential forces that develop in the contact area. Hence, the creep may be defined as a combined elastic and frictional behavior in which two elastic bodies that roll over each other share a contact area where both slip and adhesion occur simultaneously. The tangential contact forces, also called creep forces, are computed through non-dimensional quantities named creepages, which consist of the relative velocities between wheel and rail at the contact point normalized to the vehicle forward speed.
In this work, a lookup table based on Kalker's USETAB [53] has been used. Like in the original version, the values from the table are normalized and calculated according to the following criteria: 1) the combined shear modulus of the wheel and rail materials G is 1; 2) the Coulomb's friction limit given by F n is 1; and 3) the square root of the ellipse's semi-axes product √ ab is 1. While the table inputs are the semi-axes ratio of the contact ellipse and the normalized creepages, the outputs consist in the normalized longitudinal and lateral creep forces and the normalized spin creep moment, which are linearly interpolated during the dynamic analysis. To build the table, the software CONTACT [54], which is based on Kalker's exact three-dimensional rolling contact theory [55], has been used. The normalized creepages and semi-axes ratios were discretized in two intervals as in the original USETAB, namely 0 ≤ x ≤ 1 and 1 ≤ x < ∞ , where x is the input of the table. However, a linear and a logarithmic distribution of ten values were used for the discretization of the first and second intervals, respectively, instead of the original linear intervals with seven values. Adopting an 40 × 40 element discretization of the contact ellipse, and by considering all possible combinations of the creepages and semi-axes ratios, a total of 320,000 calculations were performed using the software CONTACT.
An important point is the consideration of an upper limit for the table to avoid inaccurate extrapolations. Furthermore, according to [53], the linear interpolations with semi-axes ratios close to zero or infinite should be avoided, since the creepage coefficients are singular in these cases. Therefore, an upper limit of 10 3 and a lower limit of 10 -3 are used for the semi-axes ratios, and an upper limit of 10 3 is adopted for the normalized creepages. If a combination of creepages and semi-axes ratio falls outside these intervals, the Polach [56] method is used to solve the mentioned singularities.

Vehicle-structure coupling system
The vehicle-structure coupling system is solved using the direct method originally proposed by Neves et al. [57] for vertical dynamics and later extended by Montenegro et al. [35] to deal with lateral behavior, in which the governing equilibrium equations of the vehicle and structure (see Sect. 3.1) are complemented with additional constraint equations (see Sect. 3.2) that relate the displacements of the contact nodes of the vehicle with the corresponding nodal displacements of the structure. A brief explanation of these two types of equations is given in the next subsections, but the readers should refer to [35] for a detailed description of the coupling formulation.

Equilibrium equations
Considering the α method [58], the equations of motion of the vehicle-structure system can be written as follows: where M is the mass matrix, R is the internal nodal force vector, F is the vector of externally applied nodal loads, a is the nodal displacements, α is the algorithmic dissipation factor proposed by Hilber et al. [58] for the α-method (taken as 0 here to correspond to the Newmark method), and the superscripts t and t + Δt indicate the previous and current time steps, respectively. After a series of mathematical operations in detail in Sect. 5.1 from the authors' previous publication [35], the final system of equilibrium equations is given by where K is the current effective stiffness matrix; D is a matrix that relates the contact forces, defined with respect to the target element coordinate system, with the nodal forces defined in the global coordinate system; Δa and ΔX are the incremental nodal displacements and contact forces, respectively; is the residual force vector that arises from the fact that the wheel-rail contact problem is nonlinear; and the superscripts i and i + 1 indicate the previous and current iteration, respectively.

Constraint equations
When contact occurs, the contact element and the target element are coupled in the three directions. Thus, the following constraint equations must be imposed: where v c and v t are the displacements of the contact (wheel) and target (rail) elements, respectively, and r is the vector of irregularities between the contact and target elements in the vertical and lateral directions. After a series of mathematical manipulations described in detail in Sect. 5.2 from the authors' previous publication [35], the final constraint equation is given by where H is the transformation matrix that transforms the displacements of the internal nodes of the contact and target elements from the global to the local coordinate system of the contact pair; and r is the modified irregularity vector that is consulted in Sect. 5.2 from [35].

Hybrid system of equations to make the couple between vehicle and structure
The equilibrium Eq. (13) and the constrain Eq. (15) form, therefore, a single hybrid system, with displacements and contact forces as unknowns, that is solved directly. Hence, the vehicle-structure interaction problem can be expressed as follows: Since the efficiency of the algorithm used for solving the system of equations is critical, a block factorization algorithm is used to solve the system of Eq. (16) that takes into account the specific properties of each block, namely, symmetry, positive definiteness, and bandwidth. This algorithm is described in detail in Appendix A of the authors' previous publication [35].
The vehicle-structure interaction numerical tool described above named "VSI-Vehicle-Structure Interaction Analysis" has been implemented in MATLAB® [46], which imports the structural matrices from both subsystems previously modeled in a finite element (FE) software, which in this work was ANSYS® [59]. As mentioned before, due to space limitations, only a brief description of the vehicle-structure coupling model is presented in this article, but a complete description of the mathematical formulation may be consulted in the author's previous publication [35].

Initial considerations
In the present section, the proposed vehicle-structure interaction model is validated with two numerical applications. First, the Manchester Benchmark organized by Shackleton and Iwnicki [45] is revisited to validate the wheel-rail contact model. The benchmark consisted of a series of tests simulated with ten different software with the aim of allowing an informed choice when selecting a contact model for a particular railway vehicle simulation scenario. The second numerical application consists of a hunting stability analysis of a suspended wheelset, in which its lateral displacements and yaw rotations are compared with those obtained with a semi-analytical model described by Wickens [60].

Benchmark description
Shackleton and Iwnicki [45] proposed a benchmark with the aim of allowing an informed choice when selecting a contact model for a particular railway vehicle simulation scenario. There is a wide range of wheel-rail contact models in the vehicle simulation software, and, to achieve acceptable computational times, all of them make simplifying assumptions. As a result, each model has a limit of its validity and restrictions to its applications that are not always apparent to the user. Thus, the Manchester Metropolitan University conducted a series of tests with ten railway vehicle simulation software and compared the results. These software tools vary in the way they establish the position of the contact point between the wheel and the rail; in this way, they predict the size and shape of the contact area and in terms of the methods used to simulate the forces that are generated in the contact interface. Table 1 summarizes the formulations adopted by each software to solve the contact problem.
The tests consisted of prescribing, both statically and dynamically, lateral displacements and yaw rotations to a single wheelset to analyze its behavior. Four case studies were conducted during the benchmark with real wheel and rail profiles, S1002 wheel and UIC60 rail with a 1:40 inclination, as depicted in Fig. 5, and a vertical load of 20 kN applied at the center of mass of the wheelset. Note that the thread-flange transition is considered without any simplification, so that the concave contact search algorithm proposed in the present article is fully tested with this application. These case studies are the following: • Case A1.1: The wheelset is subjected to a prescribed lateral displacement from 0 to 10 mm with 0.5 mm increments. A static analysis is performed in each position, and the normal contact is evaluated. • Case A1.2: The wheelset is subjected to the previously described lateral displacements combined with a yaw rotation from 0 to 24 mrad with 1.2 mrad increments. A static analysis is performed in each position, and the normal contact is evaluated.
• Case A2.1: Forward speed of 2 m·s -1 is given to the wheelset on straight track while it is subjected to the previously described lateral displacements. A dynamic analysis is performed, and both the normal and tangential contacts are evaluated. • Case A2.2: The wheelset is subjected to the combinations of lateral displacements and yaw rotations described in the case A1.2. The dynamic conditions are the same as for the case A2.1, and both the normal and tangential contacts are evaluated.
The results obtained with the proposed model described in Sects. 2 and 3 are compared with those obtained with the tested software shown in Table 1. The local coordinate systems considered in the benchmark, as well as the adopted conventions, are described in Appendix B of [45]. Each tested code has been assigned a line/marker style (see Fig. 6), being the results obtained with the proposed model superimposed over those obtained with the tested software.

Rolling radius difference
The rolling radius difference ΔR between the left and right wheels obtained in the test case A1.1 is plotted in Fig. 9.  Again, a good agreement can be observed between the proposed model and the tested software. For lateral displacements between 4.5 and 6 mm, the contact occurs in the tread/ flange transition of the right wheel, causing a small increase in the rolling radius. After 6 mm, however, the contact point jumps to the flange, leading to an abrupt increase in the right rolling radius and, consequently, in the rolling radius difference between the two wheels.  obtained with VOCOLIN derive from the non-consideration of the roll rotation of the wheelset to locate the contact point, while the trend followed by LaGer and CONTACT PC92 is justified by the fact that the output given by these codes is related to the wheelset coordinate system rather than to the track centerline coordinate system (see Appendix B of [45] for the definition of these coordinate systems).

Longitudinal creepages
The longitudinal creepages in the left υ ξ,lft and right υ ξ,rht contact interfaces obtained in the dynamic test case A2.1 are plotted in Fig. 11. According to [45], the lack of conformity between the results predicted by the several codes is due to differences in the way the total longitudinal creepage is distributed between the left and right contact interfaces. Therefore, in absolute terms, the output of the proposed model is in agreement with the outputs obtained with the software tested during the benchmark, except CONPOL, which follows an isolated trend.

Lateral creepages
The lateral creepages in the right contact interface υ η,rht calculated in the test case A2.2 are presented in Fig. 12 (the lateral creepages obtained in the left contact interface have not been published in [45]; therefore, only the results of the right interface are presented). The results given by the proposed model accompany the main trend followed by all codes, except CONPOL, which again shows a different output. These differences are justified by the fact that CON-POL neglects the effects of the yaw angle of the wheelset in the calculation of the creepages. This important limitation also affects the spin creepages, as will be shown in the next section.

Spin creepages
Finally, the spin creepages in the left υ ϕ,lft and right υ ϕ,rht contact interfaces obtained in the test case A2.2 are plotted in Fig. 13. The spin creepages follow the same trend as the contact angle (see Fig. 10), since they depend directly on it. Therefore, the different trends observed in the left side are justified by the same reasons presented in Sect. 4.2.4, while the discrepancies obtained with CONPOL in the right side are, once more, due to the non-consideration of the yaw angle effects in the calculation of the creepages. Regarding the proposed model, a good agreement is observed between the results obtained with it in the right contact interface, and those obtained with the software tested during the benchmark, with exception of CONPOL for the reasons stated above. In the left side, the proposed model follows again the same trend as GENSYS and VAMPIRE.

Final remarks
Although a general agreement between the tested software and the proposed model is observed, there are, in some cases, notable discrepancies. However, the main discrepancies are mainly justified by limitations of the contact models adopted by some of the tested software, especially CONPOL and VOCOLIN [45], rather than by limitations of the proposed model. Moreover, the results obtained with the proposed model are, in most cases, in an excellent agreement with those obtained with GENSYS, NUCARS, and VAMPIRE, which are widely used in dynamic simulations of railway vehicles. Therefore, it can be concluded that the wheel-rail contact model developed in this work is suitable for being used in railway dynamics applications.

Hunting phenomenon
Due to the specific conic shape of the train wheels, when a wheelset is running on a straight track and is subjected to a lateral perturbation, the rolling radii of the left and right wheel differ from each other. Hence, since both wheels have the same angular velocity if the wheelset is running with a constant speed, the wheel with larger radius will experiment a higher velocity than the opposite wheel. This phenomenon will force the wheelset to yaw and go back to the centered position, making the rolling radius of the opposite wheel to become larger. This process, called hunting motion, tends to continue indefinitely in an unsuspended wheelset making it unstable [47,60]. However, the creep forces that arise in the contact interface act as damping forces that dissipate energy and ensure the existence of a certain range of speeds where the wheelset is stable. The speed above which the wheelset become unstable is called critical speed. In addition to the creep forces, the critical speed of a wheelset also depends on the wheel conicity, wheelset mechanical properties, and suspensions. The last one is particularly important to ensure that the wheelset instability occurs only at higher ranges of speeds. Klingel [61] derived the following expression based on purely kinematic relationships, which describes the wavelength λ Klingel of the hunting motion a non-suspended wheelset: where R 0 the initial rolling radius of the wheel, L cp is the half lateral distance between contact points, and γ 0 is the conicity of the wheels. However, as aforementioned, this expression only considers the kinematic components of the movement, ignoring the inertial effects due to the mass of the wheelset, the influence of the suspensions' flexibility, the creep forces that arise in the contact interface, and the real shape of the wheels (the wheels are perfectly conical in the Klingel's model). Hence, to allow a reliable validation of the proposed model, a fully dynamic model of the wheelset is adopted in this second validation study.

Numerical models
Three different numerical models have been considered in the present study, namely (see generic model in Fig. 14): i) model A, suspended wheelset without suspensions and with conic wheels; ii) model B, suspended wheelset connected to a moving frame by the primary (lateral and longitudinal) suspensions and with conic wheels; and iii) model C, equal to model B but with a realistic wheel with thread and flange (Shinkansen wheel with diameter ϕ of 860 mm and JIS60 rail profile, as shown later in Fig. 20, and used in the experimental validation example from Sect. 5).
The geometrical and mechanical properties of the model are presented in Table 2. Note that the contact ellipse semi-axes a and b, as well as the Kalker creepage coefficients c 11 and c 22 , are calculated for a static position of the wheelset centered with the track and maintained constant throughout the analysis.

Governing equations of motion of the semi-analytical model
The results obtained with the proposed model are compared with those obtained with a semi-analytical model described by Wickens [60] and based on the following simplifying assumptions: a) The wheelset is rigid and is connected to a reference moving frame by lateral and longitudinal suspensions. b) The running speed of the wheelset is constant. c) The wheelset movement is characterized exclusively by two degrees of freedom: the lateral displacement y ws and the yaw rotation ψ ws (see Fig. 14). are valid, and the dimensions of contact area remain constant throughout the analysis. f) The slip inside the contact area is neglected, being the tangential contact problem solved with the Kalker's linear theory [62]. g) The secondary effects, such as gravitational stiffness, gyroscopic effects, and spin creep, are neglected.
Based on these assumptions, the linear equations of motions that govern the dynamics of the system can be written as follows: where V is the forward speed of the wheelset, f x = Gc 11 ab and f y = Gc 22 ab . The remaining variables present in Eq. (18) are described in Table 2. The system of linear differential Eq. (18) can be solved using a direct integration method. Note that for model C, in which the wheel is formed not only by the thread, but also by the flange, the governing equations of motion described in Eq. (18) are only valid until the moment the flange hits the rail.
The speed above which the wheelset become unstable, called critical speed V crit , can be determined from a stability study described in detail in [63] and is given by For running speeds below the critical value, the wheelset experiences a sinusoidal lateral motion that tends to damp out if no further disturbances occur. However, if the critical speed is exceeded, the wheelset undergoes an increasing oscillatory motion that makes it unstable.
The hunting wavelength is also a characteristic of the hunting motion of the wheelset, since is independent from the running speed. By performing a quasi-static analysis of the dynamic equations of motion, the theoretical hunting wavelength λ theory is found to be [63] Notice that the theoretical hunting wavelength expressed in Eq. (20) becomes equal to the wavelength proposed by Klingel, defined in Eq. (17), if the dynamic terms are neglected.

Bases of the analysis
The results obtained with the numerical integration of the system of Eq. (18) are compared with those obtained with the proposed vehicle-structure interaction method proposed in the present article. The time step used in all the analysis with both the numerical and semi-analytical models is Δt = 0.001 s. The Newmark integration scheme with integration parameters α = 0, β = 0.25, and γ = 0.5 is used to solve the equations of motion. At the beginning of the dynamic analysis, a lateral impulsive load of 10 kN is applied at the center of mass of the wheelset to drive the system away from its equilibrium position, causing it to oscillate over the track centerline. Since the equations that govern the dynamic behavior of the analytic model assume that there is no slip inside the contact area (see Eq. (18)), the Kalker's linear model [62] is adopted to compute the creep forces in the simulations performed in the present section. The results obtained in the dynamic analyses using the three models described above are presented in the next section.

Validation results
To investigate the influence of the suspensions in the wheelset stability, a dynamic analysis using model A has been performed. Figure 15 shows the lateral response of the wheelset using model A for two different speeds. Without the inclusion of suspensions, the wheelset is always unstable, even for low speeds, such as V = 20 km/h. Moreover, the numerical hunting wavelength λ num analyzed for all the speeds correctly matches the theoretical value calculated with Eq. (20) (λ num = 22.639 m vs. λ theory = 22.745 m). In terms of the response, a very good agreement between both formulations is also observed. On the other hand, when the suspensions are considered in the model, the wheelset's behavior clearly improves in respect to its hunting stability, as it is shown in Fig. 16. This figure shows the wheelset response for a running speed of V = 100 km/h in terms of lateral displacement and yaw rotation when considering model B. It is clear that the wheelset is running below its critical speed, since the energy dissipation due to the creep forces and the stability provided by the primary suspensions lead to a decrement of the hunting amplitude. Such conclusion is supported by the analytical evaluation of the critical speed using Eq. (19) which, for this particular case, is V crit = 234.4 km/h. An excellent agreement between numerical and semi-analytical results is observed, proving the accuracy of the numerical tool.
It is important, therefore, to numerically analyze the behavior of the wheelset using model B for a running speed above the critical value. The results obtained with such analysis are plotted in Fig. 17, in which the dynamic response of the wheelset for V = 250 km/h is presented. As expected, the movement is now unstable, and the system does not return to its centered position. Again, for analysis performed for speeds above the critical value, the wheelset's behavior is clearly captured by the numerical tool with slight differences in relation to the analytical model that has been previously justified.
Regarding the hunting wavelength, which is independent from the running speed, the numerical values obtained for the analysis below and above the critical speed are, respectively, λ 100 = 22.583 m and λ 250 = 22.708 m. These values are in a good agreement with the theoretical value calculated with Eq. (20), which is found to be λ theory = 22.754 m.
The lateral displacement of the wheelset for three different running speeds, including the critical speed, is plotted in Fig. 18a. As the speed increases, the oscillation decay rate tends to decrease, reaching a null value at the critical speed, proving that the proposed numerical tool clearly captures the dynamic behavior predicted with the analytical expression for evaluating the critical speed defined in Eq. (19). After that, the hunting motion grows indefinitely, and the behavior of the wheelset becomes unstable. The critical speed is, therefore, a transition in the dynamic behavior of the wheelset that can be analyzed with the logarithmic decrement factor δ n , given by where y ws (t) and y ws (t + nT) are two peak displacements separated by n consecutive cycles with period T. The logarithmic decrement related to the responses of the wheelset for speeds ranging from 50 to 300 km/h is depicted in Fig. 18b. As expected, the logarithmic decrement for the lower speeds is positive, but starts to decrease as the speed increases. Once the critical speed is reached, the decrement becomes null, since the hunting motion maintains the amplitude throughout the analysis. Then, once the speed exceeds the critical value, the decrement turns negative, and the wheelset experiences an unstable behavior.    To validate the two-point contact scenario described in Fig. 1, the same situation shown in Fig. 17 is simulated using model C and compared with the analytical model in Fig. 19. It can be observed in Fig. 19a that, until t = 2.9 s, the numerical model follows the semi-analytical solution. However, after this point, the gap of approximately 5.7 mm that exists between the right wheel flange and the lateral side of the rail closes, and the wheelset lateral movement becomes limited by the flange. Naturally, afterwards, the numerical solution drives away from the analytical one. By observing Fig. 19b, where the left and lateral contact forces are superimposed on the relative lateral displacements, it is clear that the peaks observed in the lateral contact forces are due to the wheel flange impacts with the lateral side of the rail, since they occur at the same moment as the gap between the flange and rail closes. Figure 19c and d depicts the relative position between the left wheel and rail in two distinct moments indicated in Fig. 19b when the flange is not touching the rail and when it hits the rail, respectively. It is clear that the high-frequency lateral contact forces are fully related with the lateral impacts between the flange and the inner face of the rail. Such result proves the efficiency and accuracy of the model to capture not only the tread contact, but also the flange contact that occurs for higher values of relative lateral displacement between wheel and rail.
Finally, with the objective of comparing the previous model presented in [35] with the one presented in this article that considers the detection of contact points in concave regions, Fig. 20a and b depicts the lateral response and lateral contact forces, respectively, obtained for the same example analyzed before (model C) but considering both contact models. The profiles used by both models are depicted in Fig. 20c, where it is possible to observe the simplification in the wheel profile geometry. As it can be observed, when the concave region is neglected, the moment where the first lateral impact occurs is different, showing that when simplifying the wheel profile geometry, some differences may occur in the response. With the simplified profile, since the flange geometry is not so smooth, the stronger lateral impact occurs more suddenly and, therefore, sooner. Naturally, after the first impact, the response deviates from each other, leading to different results. Such comparison shows that using simplified approaches might lead to less accurate and realistic responses.

Initial considerations
In addition to the numerical and analytical validations shown in the previous section, the present article also presents a validation based on an experimental test that took place in the Railway Technical Research Institute (RTRI) in Japan. This test, which consisted of the analysis of the dynamic behavior of a full-scale railway vehicle mounted over a test rig that imposes vertical and lateral rail deviations while the train is running at different speeds, was used to draft the Displacement Limit Standard for Railway Structures [64] currently in use in Japan. This code provides specifications regarding the maximum deformations of the railway track that should be guaranteed to ensure the stability of railway vehicles running over bridges during ordinary operating and seismic conditions. A revision of the main criteria stipulated in this normative document can be accessed in [43].
In this section, and in addition to the numerical validation described in Sect. 4, the proposed model is validated based on experimental results, namely, the accelerations measured in the carbody above the rear bogie. While in one of the authors' previous publication [35], only one single example has been used to validate the model (lateral accelerations measured for a scenario where the train runs at 300 km/h), in this article, a much more comprehensive validation is presented by comparing the experimental and numerical responses of the vehicle in the vertical and lateral directions caused by track deviations imposed by the actuators on both directions and for train speeds ranging between 100 and 400 km/h.

Track
Since the test has been performed on a rolling stock test plant, in which the railway vehicle runs over four wheel-shaped rails connected to the actuators, the track is modeled as rigid to be consistent with the test rig's characteristics, while the rail deviations imposed by the actuators are considered as track irregularities and included in the dynamic analysis through the r vector described in Sect. 3. Therefore, the results presented later in Sect. 5.3 only refer to the vehicle's response. It is important to highlight that, although this example considers a rigid structure, the model proposed here is based on a FE environment (see Sect. 3), which allows the analysis of flexible structures with any degree of complexity, such as those studied in several authors' previous publications [17,38,65,66].
During the tests, the actuators imposed two types of rail deviation geometries that aimed to simulate the deflection of consecutive bridge spans as rigid bodies. These geometries, which were used in the dynamic vehicle-structure interaction analyses that gave rise to the displacements limits imposed by the Japanese code [64] to ensure traffic stability under ordinary and seismic conditions, are the bending shape (BS) and the parallel shift (PS), as shown in Fig. 21. While the first simulates the relative rotation between two adjacent spans, the latter reflects the movement of a single span rotation, keeping the adjacent spans parallel. In the authors' previous publication [35], only a lateral deviation corresponding to span lengths L of 20 and 40 m was tested, while, in this work, the model validation considered a much wider range of scenarios, namely, four  (10,20,40, and 60 m) and lateral and vertical BS and PS rail deviations with maximum amplitude δ of 8 mm (see Fig. 21, which shows the deviation geometry for the lateral direction for exemplification purposes).
To guarantee a smooth geometry of the rail in the decks' extremities and avoid unrealistic impacts in the numerical calculation, transition zones have been added to the deviation geometries mentioned above, as depicted in Fig. 22. Here, L t denotes half-length of the transition, while θ t is the span rotation, and x t is the distance from the start of the transition. According to [35], the vertical and lateral deviation geometry of the transition zone r t imputed as a track irregularity in the numerical model is given by in which β t is the relative bending stiffness of the rails and pads in the vertical or lateral directions, given by where k p is the pad stiffness, E is the Young modulus of the steel, and I r is the moment of inertia of the rail on both directions (see Table 3).

Vehicle
The vehicle used in the experiments carried out in the test rig consisted of a full-scale passenger carriage, whose properties were provided by the RTRI. 1 The numerical model has developed in ANSYS® [59], using rigid beams to consider the rigid body movements of the vehicle connected through spring-dashpots elements to simulate the primary and secondary suspensions. Mass point elements were included in the center of gravity of each rigid body of the vehicle, namely one at the carbody and one in each bogie and each (22) wheelset, to take into consideration both their mass and inertial effects. The wheel profile is a conic and arc profile wheel with diameter ϕ of 860 mm, same as that used in the Shinkansen trains, while the rail corresponds to JIS60 profile, as shown previously in Fig. 20.

Validation results
The results obtained with the vehicle-structure interaction tool presented in Sects. 2 and 3 are compared with those obtained in the experimental test described before. The validation is carried out based on the comparison between the numerical and experimental accelerations obtained in the carbody above the rear bogie. All the results presented in this section regarding the vertical accelerations are obtained exclusively with the imposition of vertical deflections to the track, while the results relative to the lateral accelerations are obtained with the imposition of transversal deflections.
The time-histories of the vertical accelerations measured in the rear part of the carbody relative to a scenario in which the vehicle is running at 300 km/h over vertical rail deviations corresponding to span lengths of 20 and 40 m are plotted in Fig. 23, while the analogous results but for the lateral direction are depicted in Fig. 24. A general good agreement can be observed for both analyzed directions between the measured data, and the numerical results using both the previous model presented in [35], and the one presented in this article that considers the detection of contact points in concave regions with realistic wheel profiles. The differences observed with the experimental data may be caused by the fact that the numerical model of the vehicle used in this work does not account for any kind of structural flexibility that may exist in the carbody. Nevertheless, the global behavior of the responses is well captured by the proposed model, proving its capacity to simulate generic railway dynamic scenarios. Regarding the comparison between both numerical models, the differences are negligible, since the consideration of the realistic wheel profile does not significantly impact in the carbody's response due to the filtering effect provided by two levels of suspensions. Figure 25 shows a comparison between the two numerical models but at the wheel-rail contact interface level, namely, in terms of wheel-rail contact forces (due to space limitations, only the BS scenarios are represented, but the conclusions to the remaining ones are the same). As it can be observed, the only differences occur on the lateral contact forces corresponding to the scenario with 20-m span with lateral irregularity (see Fig. 25c). This is due to the fact that it is the only scenario where lateral flange impact takes place, which highlights the differences between the two models (note that in the scenario with 40-m span depicted in Fig. 25d, the lateral contact forces are much smaller because there are no flange impacts). This conclusion agrees with that drawn in the application presented in Sect. 4.3. Naturally, when the irregularity is imposed only in the vertical direction (see Fig. 25a and b), the responses are the same because there is only contact in the tread region, which does not differ on both models. Hence, depending on the level of the lateral impacts, these differences may be higher or smaller. However, it is important to stress that in scenarios in which the lateral contact forces pose a significant role in the train running safety, such as in the presence of crosswinds or earthquakes, the consideration of a more realistic contact model leads to more accurate results and, consequently, to more accurate predictions of the safety performance.
To further enhance the model validation, Figs. 27 and 28 present, respectively, the maximum vertical and lateral accelerations in the rear part of the carbody obtained in the experimental tests and with the proposed model with realistic wheel profiles. The results refer to a wider range of span lengths, namely, 10, 20, 40, and 60 m and vehicle speeds ranging between 100 km/h and 400 km/h with steps of 50 km/h (the experimental and numerical results for each span length have been assigned with a specific marker and line, respectively, as shown in Fig. 26). Again, the numerical results show a good agreement with the experimental results, with some discrepancies which may be justified by the inaccuracy of the vehicle's numerical model to reproduce the effects that arise from the flexibility of the carbody. However, as stated before, the global trend is perfectly captured by the numerical model, which reinforces once again its ability to be used in railway dynamic applications, such as the study of train running safety [17,36,37] or passenger riding comfort [18].

Conclusions
In a previous work from the authors, the contact point detection has been performed through the solution of a pair of nonlinear geometrical equations. However, this approach does not give a single solution if one of the contact surfaces is concave. Therefore, this paper proposes an additional algorithm, called the concave contact search, to deal with contact points located in the tread-flange concave transition zone. This algorithm consists of discretizing the wheel and rail profile surfaces in points to evaluate the potential contact between them by determining if any of these points lie inside the opposite surface, forming the so-called intersection volume. The validation process is of the utmost importance for any new developed tool; however, most of the existing models in the literature related with VSI analysis do not thoroughly present it, making it hard to check which examples and tests were used to guarantee the correct correspondence between the model performance and reality. Thus, the VSI model developed in the present article is validated with three applications, two numerical and one experimental, from which the following conclusions were drawn: • In the first application, the Manchester Benchmark is revisited to validate the concave contact search proposed in this article. The benchmark comprised a series of tests that consisted of prescribing, both statically and dynamically, lateral displacements and yaw rotations to a single wheelset to analyze its behavior. Several contact characteristics were analyzed during the benchmark, namely, the contact point positions on both wheels of the wheelset, the rolling radius difference between wheels, the contact angles, and the creepages. The results obtained with the proposed model for all the analyzed quantities showed an excellent agreement with those obtained with other railway vehicle dynamics multibody software, such as GENSYS, NUCARS, and VAMPIRE. The few discrepancies observed are mainly justified by limitations of the contact models adopted by some of the tested software, especially CONPOL and VOCOLIN, rather than by limitations of the proposed model. Marker and line styles assigned for the experimental and numerical results, respectively, with respect to each span length • The second validation example consists of evaluating the lateral stability of a single wheelset running at several speeds. The dynamic response of the wheelset calculated with the proposed model is compared with that obtained using a semi-analytical model with two degrees of freedom available in the literature. A good agreement between the responses obtained with the proposed model and those obtained by the integration of the equations of motion of the semi-analytical model is observed. As expected, for speeds below the critical limit, both the lateral displacement and the yaw rotation of the wheelset tend to damp out after being driven away by a lateral disturbance. This is due to the energy dissipation provided by the creep forces and to the stability provided by the suspensions. However, when the speed exceeds the critical value, the behavior of the wheelset becomes unstable, leading to a hunting motion that grows indefinitely. The critical speed predicted by the proposed formulation using a logarithmic decrement factor is also in a perfect agreement with the theoretical value determined from a stability study described in the literature. Moreover, the two-point contact scenario (flange and tread contact points) was also validated through an example in which a flange was introduced to the wheel profile. The results showed that, exactly after the moment when the gap between flange and rail closes, the wheelset lateral movement became limited by the flange, and an abrupt increment in lateral contact forces occurred, demonstrating the capabilities of the model on capturing the flange impacts with the rail. Finally, differences were detected in the magnitude of the contact forces and in the moment when the first flange impact occurs between the simplified model presented in one of the authors' previous publications and the model presented here that considers the concave region of the wheel profile. • Finally, an experimental test conducted in a test rig from the RTRI, in which a full-scale railway vehicle runs over a track subjected to vertical and lateral deviations, is reproduced numerically. Unlike in one of the authors' previous publication, where only a single example has been used to validate the model (lateral accelerations measured for a scenario where the train runs at 300 km/h), a much more comprehensive validation is presented in this paper, consisting of the comparison between the experimental and numerical responses of the vehicle in the vertical and lateral directions caused by track deviations imposed by the actuators on both directions and for train speeds ranging between 100 and 400 km/h. The results showed a good agreement between all the examples, proving the capacity of the model to deal with generic railway dynamic applications. Regarding the differences between the simplified and the realistic model presented here, they are only notorious during the occurrence of lateral flange impacts, because it is in these moments that the contact points may be located in the concave transition zone.