Simulation of wheel and rail profile wear: a review of numerical models

The development of numerical models able to compute the wheel and rail profile wear is essential to improve the scheduling of maintenance operations required to restore the original profile shapes. This work surveys the main numerical models in the literature for the evaluation of the uniform wear of wheel and rail profiles. The standard structure of these tools includes a multibody simulation of the wheel–track coupled dynamics and a wear module implementing an experimental wear law. Therefore, the models are classified according to the strategy adopted for the worn profile update, ranging from models performing a single computation to models based on an online communication between the dynamic and wear modules. Nevertheless, the most common strategy nowadays relies on an iteration of dynamic simulations in which the profiles are left unchanged, with co-simulation techniques often adopted to increase the computational performances. Work is still needed to improve the accuracy of the current models. New experimental campaigns should be carried out to obtain refined wear coefficients and models, while strategies for the evaluation of both longitudinal and transversal wear, also considering the effects of tread braking, should be implemented to obtain accurate damage models.

close to each other, and the shape of the resulting contact patch is not elliptic [3]. Conformal contact instead refers to a contact condition in which the contact patch area is not planar, and its dimensions are not negligible with respect to the dimensions of the contacting bodies [4]. In such condition, the hypothesis of elastic half-space, which holds in the Hertzian theory and in many other contact theories, is hence violated. An example of conformal contact in the railway field is the contact between the wheel tread and the brake shoe during tread braking operations.
When the wheel and rail profiles tend to become locally more conformal with respect to the original contact conditions due to the material removal, the wheel tread can feature a hollow shape. In such condition, the contact point location becomes an irregular function of the lateral displacement, i.e. the contact point position "jumps" as the wheelset shifts laterally, and this inevitably produces vibrations and rolling noise [5,6].
Shifting focus to irregular wear, surface shape changes in longitudinal direction cause wheel out-of-roundness [7,8] and rail corrugation [9,10], which can lead to impact loads and rolling noise as well as to possible damages to structural components of the wheelset and the track.
To avoid the drawbacks related to changes in the wheel and rail profiles in both tangential and longitudinal directions, the profiles of wheels and rails are periodically restored through maintenance operations, namely wheel turning and rail grinding. Wheel turning operations are usually scheduled at fixed running distances, while rail grinding intervals are typically related to a certain amount of line traffic, expressed in gross ton mile (GTM). Nevertheless, these operations are expensive and cause loss of performance for the entire railway system, since wheel maintenance prevents the normal vehicle circulation, while rail grinding clearly leads to temporary traffic interruptions.
Therefore, a deep knowledge of wheel and rail damage mechanisms is the key to improve the wheel-rail interface tribological behaviour, by minimizing wear and RCF negative effects, as well as to optimize maintenance scheduling. Recently, a massive experimental campaign was performed to investigate transversal and longitudinal wear on wheel and rails on Chinese high-speed lines [11]. Nevertheless, such experimental tests demand a huge amount of data as well as long times and high costs. Hence, the high computational power of modern computers allows to address these goals by means of numerical simulations. Different strategies, even based on co-simulation and code parallelization, can be implemented to develop effective numerical tools.
Due to the huge extent of the topic, the present paper is mainly focused on a review of modelling strategies for the calculation of the wear of wheel and rail transversal profiles, thus neglecting the RCF phenomenon and the models for the calculation of wheel out-of-roundness and rail corrugation. In fact, for the sake of computational speed, usually models intended to compute the change of the transversal profile spread the worn material uniformly along the running direction, while models aiming to compute the changes in the longitudinal direction assume a constant wear along the transversal direction. Furthermore, also the wear of wheel and rail profiles in railway turnouts is outside the scope of the present paper, as the changes of the rail geometry and stiffness along the running direction require the introduction of special solutions in the view of a proper modelling. For rail profiles, the paper thus only focuses on models which consider constant rail profiles for the whole length of the considered track section, hence variations of the rail profile along the longitudinal direction [12] and rail wear phenomena due to local irregularities [13] are neglected.
The wear of wheels and rails is a strongly nonlinear phenomenon, depending on several variables, namely operating conditions, initial profiles, materials involved in the process, contact conditions, etc. Although several wear models exist in the scientific literature, with different purposes and goals, most of them share the following challenges in wear numerical simulations: 1. Solution of the vehicle-track coupled dynamics, with the current standard being the simulation based on multibody (MB) codes. 2. Calculation of the forces and creepages acting at the wheel-rail interface, which can be performed either globally or locally, i.e. considering or neglecting the distribution of stresses and relative speeds on the contact patch. 3. Calculation of the worn material according to a suitable wear law, correlating the contact patch parameters with the wear rates. This task can be accomplished by computing the worn material either locally in each cell of the contact patch or globally. 4. Optional update of the original profiles.
Since different solutions are witnessed to address the above challenges, the present paper aims to give a comprehensive review of the numerical models developed by several scholars and even railway companies to simulate the process of wheel and/or rail wear along the transversal profile, considering uniform shapes in the wheel circumferential direction and rail longitudinal direction. Clearly, wear of the wheel and rail surfaces is a slow process, which occurs on an extremely large scale in terms of vehicle mileage and load burden by the track. Hence wear models commonly assume that the transversal profiles can be considered as unchanged during a wear step of the simulation and then they introduce a simplifying hypothesis to magnify the amount of worn material and to update the profile. 405 Simulation of wheel and rail profile wear: a review of numerical models Figure 1 shows four different strategies that can be implemented in wear models to update the profiles. Please note that in Fig. 1 L and T stand for the current simulated track length and time, respectively, dT and dL refer to the increase in time and travelled length, while L f and T f correspond to the final track length and simulated time. The same symbols are used in all similar figures throughout the paper. The technique in Fig. 1a   followed by an update of the worn profile. This strategy clearly reduces the simulation time, but it can lead to spikes in the distribution of the wear depth along the transversal direction. An upgrade of this method can be the solution in Fig. 1b, in which after a single dynamic simulation, the profile is updated with an iteration scheme introducing a smoothing operation, until the desired mileage is reached. Please note that the approach sketched in Fig. 1b differs from the one shown in Fig. 1a in the number of calls to the wear module. In fact, the strategy depicted in Fig. 1a performs a single computation of the worn profiles, while the strategy sketched in Fig. 1b relies on an iteration of wear computations, and at each iteration a new coupling of the wheel and rail profiles is carried out, considering the profiles obtained at the previous step. The most common solution however is the one sketched in Fig. 1c, where an iteration of dynamic simulations with updated profiles is implemented. Finally, the most accurate solution, which however can require a huge computational time, especially if proper parallelization and co-simulation techniques are not introduced, is the online update of the wheel profile during the dynamic simulation, depicted in Fig. 1d.
Therefore, in the present review paper, the numerical models for the calculation of the wear of wheel and rail profiles in transversal direction are classified according to the strategy implemented for the worn profile update, according to the four strategies presented in Fig. 1. Furthermore, the key points and peculiarities of each surveyed model are brought out, describing the aim of the computation (wheel or rail profile or both) and possibly highlighting the cosimulation approach adopted. In fact, modern computer and networks allow for a speed up of the calculation thanks to co-simulation, thus performing different calculations in different software (and even hardware) environments. According to the authors' research group, three co-simulation approaches can be implemented in railway wear models, as shown in Fig. 2 [14]. The approach in Fig. 2a does not require a co-simulation, and both the dynamic simulation and the wear module are implemented in a MB software. The second strategy, shown in Fig. 2b, instead includes a MB model for the dynamic part and a separate module, implemented in a different software, which performs the calculation of the contact forces, the calculation of the worn material and the profile update. Finally, an intermediate solution, sketched in Fig. 2c is based on an external module which performs the post-processing calculation of the wear depth and the profile update, starting from the outputs of the MB contact module.
The present paper is organized as follows. First, the paper briefly focuses on the tribology of the wear mechanism, by highlighting the main phenomena involved. Then, attention is shifted to the wear laws used to relate wear to the contact conditions. Such laws typically express wear as a function of either the normal load or the work dissipated at the wheel-rail interface. Subsequently, the emphasis is placed on the review of the wear models in the literature. Finally, different models and strategies are compared, and also the future trends and challenges in the topic are underlined.

Tribology of wear mechanism
When sliding occurs at the interface of two bodies in contact, a material loss always occurs at a rate depending on the contact pressure, the sliding velocity, the adhesion level, the material properties, the temperature at the contact interface, etc. However, different mechanisms can be involved in the wear process. In 1990 Ashby and Lim [15] reviewed the main wear mechanisms for metals, according to the values of contact pressure, sliding velocity and temperature. The authors distinguished "hot" mechanisms, due to an increase in bulk and local temperatures of the contacting bodies and leading to melting and severe oxidative phenomena, and "cold" ones, less severe and dominated by plastic phenomena.
The wear of wheel and rails has been a big concern since the dawn of railway operations and much work was first carried out by researchers to better understand the wear mechanism from a metallurgical point of view. A great effort was provided by the British railway centre (BRR), which promoted several investigations starting from the early 50 s, with the works by Dearden, who investigated the main factors involved in rail steel wear and pointed out the economical and performance issues related to wear, also providing typical values of wheel and rail steel wear rates and suggesting simple methods to reduce the total wear [16,17]. The experimental campaign performed by Beagley using an Amsler wear machine led to the identification of two wear regimes, referred to as "mild" and "severe" [18]. The former was characterized by bright surfaces of the rollers and debris consisting of oxides, while the latter featured a matt grey surface of the rollers with debris consisting of metallic iron. Beagley suggested that the transition between the two regimes took place when the maximum shear stress due to rolling-sliding contact exceeded the shakedown limit. The same regimes were identified by Bolton et al. [19], who observed that in the mild regime, wear was a function of the normal load only as long as the friction coefficient was saturated, while wear in the severe regime proved to be proportional to the term T ⁄A, where A is the contact area and Tγ is the product of tangential force and creepage. A deeper investigation from the same research team [20] led to the identification of a third wear regime at high creepage values, characterized by ploughing scores on the disc surfaces and by a high wear rate. The authors referred to the three wear regimes as modes I, II and III, respectively, with modes I and II roughly corresponding to the mild and sever regimes previously defined. This work confirmed that the wear rate in modes I and II had a good correlation with the following expressions: where Δm is the mass lost, d is the distance rolled, N is the normal contact load, and k I , k II,1 and k II,2 are model constants fitted on the experimental data.
More recently, Lewis and Dwyer-Joyce [21], from Sheffield University, identified again three wear regimes, named as "mild", "severe" and "catastrophic", respectively, using a twin disc rig. The authors deeply focused on the transition mechanisms among the three regimes and showed that the transition between mild and severe modes occurred when the contact patch was in full slip conditions, while the second transition between the severe and catastrophic modes took place when the global surface temperature was in the range 200-250 °C, whereby also the yield strength of the wheel steel significantly decreased. Although Lewis and Dwyer-Joyce identified three regimes, their severe mode did not correspond exactly to the type II wear defined by Bolton and Clayton. In fact, in the BRR type II wear, the wear rate showed a linear dependency on the term T ⁄A, as stated by Eq. (2), while the severe mode described by Sheffield researchers featured a constant wear rate as a function of the term T ⁄A. Finally, Zhakarov et al. [22] noticed that by expressing the wear rate as a function of the parameter p z γ, where p z is the normal pressure, a fourth wear regime could be noticed,  2 Co-simulation approaches in wear modelling: a no co-simulation; b external code for wheel-rail contact and wear estimation; c external code for wear computation named heavy. In this wear mode, the wear rate decreases if p z γ increases and darker debris particles are generated with respect to the particles produced in the severe mode. The peculiarity of the experimental work [22] is that, conversely to what the BRR team did, Zhakarov et al. investigated conditions of pure lateral creepage rather than longitudinal creepage. The results showed that similar wear rates were observed to those predicted when simulating longitudinal creepage; however, since higher wear rates are usually observed in curved track sections, whereby lateral and spin creepages dominate over longitudinal creepage, the Russian research team stated that experimental campaigns should reproduce such conditions. Please note that in the above-mentioned works, the transition among the different wear regimes was considered dependent on both the contact pressure and the creepage value, while in tribological studies performed by Swedish researchers [23], who related the worn material to the product of normal load and sliding distance, an abrupt transition towards catastrophic worn rates was predicted when the contact pressure was above a threshold limit, approximately equal to 80% of the softer material hardness. Therefore, it should be cleared whether this transition was not observed by British researchers due to the limited values of contact pressure tested, which never reached the threshold value predicted by the Swedish researchers.
As noticeable from the previous works, the tribological studies of the wheel and rail wear mechanisms started in the early 70 s and were mainly performed by means of laboratory devices, especially pin-on-discs and twin-disc machines. However, although laboratory tests allow a better reproduction of test conditions and are less expensive with respect to onfield tests, they cannot completely reproduce the real tribology of the contact, which includes third body materials as well as natural and artificial lubrication. Olofsson [24] recently pointed out the effect of contaminants at the wheel-rail interface, which also have a strong influence on wear processes and material removal rate. Wang et al. [25] carried out an experimental campaign to identify wear regimes and transitions as a function of the T ∕A term, considering the presence of third body layers at the roller interface, such as water, oil, sand and friction enhancers. Not all simulators take this aspect into considerations; however in the next sections, focus will be given to some strategies developed by various authors to reduce the calculated wear considering lubrication effects.

Wear laws
As the wear of wheel and rail steels can occur according to several different mechanisms, mainly oxidation, plastic deformation and crack growth and propagation, the calculation of the material removal relying on an accurate description of the tribological processes would be extremely difficult to perform and would require an intensive computational effort in most cases.
As stated by Braghin et al. [26], the first attempts of calculating the wear of wheel and rails relied on basic equations relating the material removal rate to the angle of attack of the wheels. However, these models calculate no wear when the angle of attack is zero and this result is not supported by experimental data, which shows a little wear rate even when the vehicle is running on a tangent track with zero angle of attack.
Therefore, the main strategy followed nowadays in the railway field to predict the wear of wheel and rail materials consists in using simple expressions tuned on experimental data to relate wear to the main contact patch quantities. The main benefit of this approach is that once the contact patch parameters are obtained, the calculation of the worn material is straightforward with just one equation. However, since the wear laws need to be calibrated with experimental data, their range of application should be limited to the materials and contact patch conditions tested in the experimental campaign, and their extension can be misleading. Due to the high costs of on-track tests as well as their limitations in the reproduction of the experimental conditions, most data is obtained by means of laboratory devices, mainly pin-on-disc and twin-disc machines.
Nowadays, two main kinds of wear laws are used in the railway field, i.e. equations based on the normal load and expressions depending on the dissipated energy at the wheel-rail interface. Furthermore, most of these wear laws can be applied either globally or locally, i.e. either calculating the total amount of worn material from the global forces and creepages at the wheel-rail contact patch or considering the distribution of tangential pressures and slip velocity on the contact area. The global application of the wear law obviously saves much computational time, but if the worn profile needs to be updated during the simulation, a simplified hypothesis is required to spread the total calculated wear on the contact area at each time-step. On the other hand, the local strategy demands a higher computational effort, but it allows to directly obtain the wear distribution on the contact patch.

Wear laws based on normal load
Simple wear laws relating the worn material to the normal load were proposed since the dawn of tribological studies on railway wear. Bolton and Clayton suggested the wear law shown in Eq. (1), while Fries and Dávila [27] developed the exponential law as stated by where p z is the normal contact pressure, while k 1 and k 2 are two model parameters tuned on the experimental data from the BRR.
Nevertheless, the main wear law currently used by several scholars based on the normal load and sliding distance is Archard's wear law, which allows a straightforward calculation of the volume of worn material V worn : Simulation of wheel and rail profile wear: a review of numerical models where s is the sliding distance, H is the hardness of the softer contact material expressed in Pa, and k Arch is a fitting coefficient of the model. Equation (4) can be applied locally by dividing the left-hand and right-hand sides by the contact area and considering the distribution of normal pressure and sliding distance on the contact patch, as stated by where x and y are the longitudinal and transversal coordinates of the points lying on the contact patch, and z worn is the wear depth.
The wear law described by Eq. (5) is widely used in the railway field and has been adopted by the researchers of the Swedish Royal Institute of Technology (KTH), who developed a wear map to calculate the wear coefficient as a function of slip velocity v slip and normal contact pressure p z . The KTH wear map, shown in Fig. 3a, was obtained from experimental data collected by means of pin-on-disc and twin disc machines as well as from the results published by British researchers. The map identifies four wear zones, and each zone is characterized by a wear coefficient k i , with i = I, II, III and IV. Please note in Fig. 3a that the transition towards catastrophic wear (zone IV) occurs when the contact pressure is above 0.8H, with H being the hardness of the softer material. The map was built including experimental results for several rail and wheel steels [23,28]; therefore it can be considered as material independent. The wear law of the KTH researcher has hence a wide range of applicability, but it is less accurate with respect to other models which are based on the calibration of wear data of a specific wheel-rail steels pair. Therefore, for specific applications, it may be the case to find the proper coefficients by means of dedicated experimental campaigns [29,30]. Although the Archard wear law is mostly used in the form of Eq. (5) with the wear coefficient calculated from the KTH map, recently some researchers pointed out that different approaches for the calculation of the wear coefficients can improve the accuracy of the calculation. Ramalho [31] obtained proper weighing coefficients for the calculation of the wear coefficient, as a function of the maximum contact pressure, creepage and peripheral speed, while Cremona et al. [32] suggested that a universal kriging technique can extend the range of applicability of the KTH map.

Wear laws based on the dissipated energy
The wear laws depending on the frictional work at the wheel-rail contact pair arose from the hard work carried out at the BRR and described in the previous section. McEwen and Harvey [33] proposed a formulation that calculates the surface removed in each axial section of the wheel per each km rolled, ΔS w (mm 2 /km), and an equation that predicts the rail wear ΔS r (mm 2 /10 3 axle), in terms of surface removal per thousand of axle passages (see Eqs. (6) and (7)). In both formulations, wear is related to the Tγ index: where D is the wheel diameter in mm; F x , F y and M z are the longitudinal force, lateral force and spin moment on the contact area, while ξ, η and φ are the longitudinal, lateral and spin creepage, respectively. Figure 3b shows the wear laws for wheel and rail steels proposed by BRR researchers. Equations (6) and (7) were obtained by fitting experimental data obtained through disc-on-disc tests on Class D wheel tyre and BS11 rail steels. The application of Eqs. (6) and (7) should thus be limited to the materials and contact conditions tested during the experimental campaign. The main drawback of the BRR method is that it can only be applied globally, since the Tγ index is calculated from the global forces and creepages acting on the contact patch. A different law based on the specific energy dissipated at the contact patch was developed from work carried out at the Sheffield University (USFD) by Lewis et al. [34], who related the wheel wear rate calculated as the mass loss per distance rolled and per unit of contact area ẆU SFD (μg/(m·mm −2 )), to the index T ⁄A. The wheel wear law obtained from twin-disc tests performed on a R8T wheel steel disc paired with a UIC60 900A rail steel disc is described by Eq. (9) and plotted in Fig. 3c.
The USFD wear law can be applied not only globally but also locally if the index T ⁄A is computed in each cell of the contact patch and then the wear rate is calculated only in the slip region. The wear coefficients in Eq. (9) were fitted from twin-disc experimental data considering the global values of forces and creepages, but Tassini et al. [35] showed that in case the USFD wear law is applied following a local approach, wear coefficients should be slightly modified in order to guarantee a better agreement with the original results measured by USFD researchers. Nevertheless, if the wear rate is calculated with Eq. (9) applied locally, the wear depth distribution on the contact area can be determined according to Eq. (10), where Δx is the longitudinal discretization step of the contact patch and ρ is the wheel steel density.
Although the BRR and USFD wear laws calculate two different wear rates from two slightly different wear indexes, they have common points, i.e. they both identify three separate wear regimes, namely mild, severe and catastrophic and they both predict a constant wear rate in the severe wear regime.
The frictional work was used as a wear index also by researchers from Budapest University [36], who developed two models to calculate a stochastic wheel debris mass flow rate density ṁd (kg/(s·m −2 )), considered proportional to the energy flow density Ė d (W/m 2 ): The wear coefficient k Zob in Eq. (11) was a piecewise constant variable calculated as a function of the local energy flow density and showed a step transition between two regimes, mild and severe, see Fig. 3d, in which the energy flow density threshold value is indicated as Ė sw . The two models were developed for a local determination of the wear rate on the contact patch; however the first model can be referred to as an "accurate" model, since it required the determination of the distribution of the tangential stresses and of the adhesion and slip areas to calculate the energy flow density according to Eq. (12). On the other hand, the second model can be considered as an "approximate" model, since it calculated the global value of the energy flow density and spread it on the contact patch according to the distribution of the normal pressure, as stated by Eq. (13).
where p x and p y represent the longitudinal and lateral stresses; v x and v y represent the tangential slip velocities, while the spin creepage is neglected; V is the vehicle speed and p z,med is the mean value of the normal pressure acting on the contact area. Please note that Eq. (12) is such that the specific energy rate is equal to zero in the adhesion zone, where no slip occurs between wheel and rail, while Eq. (13) calculates a wear rate in each point of the contact patch, since it does not distinguish between adhesion and slip areas to save computational times. Therefore, the wear coefficient takes different values in the two models.
Finally, a wear law based on the dissipated friction is the Krause-Poll law [37], which relates the volume of worn material to the frictional work W fric at the wheel-rail contact patch by means of proportionality constants C m and C s which differ for two wear regimes, namely mild and severe (see Eq. (14)). The frictional work W fric is calculated by multiplying the frictional power dissipated at the contact patch P fric by the integration time step Δt, as stated by Eq. (15). The frictional power is simply obtained as the product of the Tγ index, calculated with Eq. (8), and the vehicle speed V, as stated by Eq. (16). The Krause-Poll law is similar to Zobory's model, since it distinguishes two separate wear regimes, namely mild and severe, by adapting proper wear coefficients. In both cases the wear coefficient is calculated as a function of the specific dissipated energy density, i.e. the contact power lost per unit of contact area. The transition between the two regimes takes place at a fixed value of the specific dissipated power, just like in Zobory's model.
A similar law, expressed in the more general form of Eq. (17), is also available in the commercial software Universal Mechanism under the name of Specht's model, where ms is a proportionality constant and lim is a threshold value.
The Krause-Poll law is implemented in the wear module of the MBS (multibody simulation) software Simpack, and it is applied globally. The Simpack wear module calculates the total worn volume and allocates 50% to the wheel and 50% to the rail, allowing the user to select either the Krause-Poll or Archard's equations. Since the software calculates a global value of the worn volume, the wear depth is spread along the width of the contact patch following an a priori distribution trend and uniformly in circumferential direction.
Although no clear evidence was found in the literature of a local application of the Krause-Poll wear law, the authors of the present paper recognize that if both left-hand and righthand sides of Eq. (14) are divided by the contact area, and the contact power lost as well as the friction work are calculated locally, a distribution of contact wear depth can be calculated according to Eqs. (18)- (20), in which ΔA is the area of an element in the global contact patch discretization grid. Please note that in Eq. (20) the spin contribution is neglected.

Final considerations on wear laws
As shown in the sections above, several laws exist in the railway field to predict the wear of wheel and rails. These equations are based on different parameters, typically related to normal load or frictional contact work, and they can compute different wear indexes, namely the total amount of worn material, in the form of worn mass or volume, or a wear rate, according to the experimental tests carried out to determine the wear coefficients. Table 1 provides a summary of the main wear laws surveyed in the previous subsections, highlighting the main features of each approach.
Although in the previous subsections the Archard wear law was presented as a law relating the worn material to the normal load rather than to the dissipated contact energy, it is interesting to notice that in the commercial MB software Universal Mechanism, the wear module is provided with a modified Archard law in which the volume of worn material V worn is computed as a linear function of the frictional work W fric , see Eq. (21), through a proportionality constant k UM . Obviously, such formulation is based on replacing N in Eq. (4) with T∕ (μ is the friction coefficient) and recognizing that the product Ts is indeed the frictional work. Therefore, the relationship between the coefficients k UM in Eq. (21) and k Arch in Eq. (4) can be described by Eq. (22). Nevertheless, the similarity between Eqs. (4) and (21) and (22) only holds when the contact patch is in full sliding conditions and the tangential load can indeed be written as the product of the normal load by the friction coefficient. When the contact patch is in partial slip, the relationship Archard KTH [38] Wear volume is proportional to normal load and sliding distance Wear coefficients are usually computed according to the KTH wear map, which includes for four regions, corresponding to three regimes Both local and global applications are possible Available in most commercial MB software packages (4) and (5) Wear laws based on the dissipated energy BRR severe wear [19] Wear mass is an explicit function of both tangential force and creepage, through the wear parameter T ∕A Wear law valid for severe regime only (2) BRR [33] Worn surface is an explicit function of both tangential force and creepage, through the wear number Tγ Only global application is possible Different coefficients for wheel and rail steels Three wear regimes are identified Wear in the severe regime is constant (6)- (8) USFD [34] Wear rate is a function of the wear parameter T ∕A , namely the wear number divided by the contact area Both global and local applications are possible Three wear regimes are identified Wear rate in the severe regime is constant (9) Zobory [36] Debris mass flow rate is proportional to the specific dissipated power Two wear regimes are considered The proportionality constant features an abrupt step as a function of the specific dissipated power Both accurate and simplified versions exist. The simplified version does not require to compute the distribution of the sliding speed on the contact patch (11)- (13) Krause-Poll [37] Worn volume proportional to the frictional work Two regimes are identified Transition between the two regimes depends on the value of the specific dissipated power Available in commercial MB software packages Possible extension to local application in Eq. (22) is not valid anymore, so the k UM coefficient in Eq. (21) must be properly tuned.
The VI-Rail wear module uses a modified Archard wear law, too, in which the worn mass Δm is related to the frictional work at the contact interface according to Eq. (23). Actually, the wear law implemented in VI-Rail is more similar to the Krause-Poll wear law, see Eqs. (14)- (17), as two regimes are identified and the switch between them depends on the value of the specific dissipated power. The default value of the switch specific power is 4 W/mm 2 , identical to the value adopted by the Simpack wear module, while the default values for the mild and heavy wear coefficients are k mild = 4.5 × 10 -7 g/J and k heavy = 1.25 × 10 -6 g/J. Although the Simpack wear modules computes the frictional work as the time integration of the dissipated power, based on the current time step size, the VI-Rail module calculates the dissipated work according to Eq. (24), as the product of the wear number Tγ and the loaded path length L p . As previously stated, the correspondence between the Archard wear law and the models based on the dissipated energy is valid only when the contact patch is in full slip. Therefore, the authors of the present review believe that naming laws based on the frictional work as "modified Archard laws" is not formally correct, also considering that in the original Archard's works [39,40] the worn volume was related the normal load and sliding distance rather than to the frictional work: In a recent work, Peng et al. [41] showed that although the BRR, KTH, USFD and Zobory wear law calculate different wear indexes, they all can be manipulated analytically in order to obtain similar expressions which relate the wear depth to a common wear parameter, using proper wear coefficients for global and local applications. Nevertheless, as the wear coefficients for all wear models were obtained with different experimental campaigns, the computed wear could vary significantly if different wear laws were used.
In 1984 Fries and Dávila [27], from Virginia Polytechnic Institute and State University, proposed a simple wear model to predict the worn material of the wheel tread on a freight vehicle running on straight track. The authors applied four wear laws in the same case study, namely a wear law based on the contact patch work, a wear law based on load and slip, similar to the Archard's model, an exponential law based on the maximum contact pressure, see Eq. (3), and the simple equation relating the worn volume to the normal load only, suggested by Bolton and Clayton, see Eq. (1). The four wear models showed a good agreement in the calculation of the wheel worn profile.
More recently, De Arizon et al. [42] compared the results obtained by two implementations of the Archard law, from Jendel and Emblom, respectively, by Zobory model and by BRR law. A comparison was first performed analytically for full slip conditions defining a common wear parameter, Ż w in m/s, which can be referred to as a wear velocity, and then numerically for rolling-sliding conditions and for a specific case study of a tramway vehicle facing curves of different radius. The authors showed that the two implementations of Archard's model provide results close to each other and to the results computed using Zobory's model, while the BRR law tends to calculate higher wear values.
As part of the AWARE (Reliable Prediction of the Wear of Railway Wheels) project, Pombo et al. [43] developed a wear model and tried to benchmark the influence of the selected wear law and of its global or local application, in terms of both prediction of worn material and computational times. The authors compared the results computed by the BRR wear law, the Archard model in the implementation of KTH and the USFD wear law in a specific case study of a railway vehicle running on the Cuneo-Ventimiglia Italian line. The authors found that in case of global application of the three laws, the BRR model computes higher wear on the flange, and this is in line with the results from De Arizon. The comparison performed in this work only dealt with the numerical results predicted by the three different wear laws, while a benchmark against experimental data and a tuning of the wear coefficients was not provided. Nevertheless, in a previous work [44], the authors had highlighted that the KTH coefficients should be reduced to give outputs in better agreement with the other models. From a computational point of view, the USFD model was the fastest, while BRR and KTH laws had similar performances and were slightly slower with respect to the USFD model. Finally, the authors showed that the application of the USFD model with a local approach led to results close to those obtained with a global approach, but the local approach was twice slower than the global one. Therefore, the authors concluded that the global application of the wear laws can be a good compromise between computational efficiency and numerical accuracy.
One last point can be mentioned focusing on the wear laws described above. In fact, these laws are typically obtained from tests performed on twin-disc machines; however Feldmaier et al. [45] showed that the wear rates predicted with a scaled two-roller machine in the severe regime corresponding to flange contact are lower than the wear rates experimentally calculated from tests on a scaled wheel machined to a real wheel profile. On the other hand, a better agreement exists in the wear rates corresponding to the mild regime, which occurs in tread contact. A possible explanation to this behaviour is that in the real wheel-rail contact, the contact pressure and slip ratio are both higher than on the two-roller machine. Therefore, great attention must be paid when a wear law obtained from a twin-disc rig is applied to the real wheel-rail contact scenario.

Wear models
The present section is the core of the paper, dealing with the presentation of the numerical models for the computation of wheel and/or rail profile wear in transversal direction. The models are surveyed based on the classification proposed in the introduction section, i.e. by grouping the numerical tools according to the strategy used for the profile update, see Fig. 1.
As previously stated, the vehicle dynamics and the wheel and rail profile wear phenomenon mutually interact, so that typically most wear computing tools include a dynamic module for the computation of the vehicle dynamic behaviour on a specific mission. This task is nowadays addressed by means of multibody simulations, usually implemented in commercial software packages. Therefore, most of the models that are surveyed in the following subsections rely on multibody simulations, which nowadays can easily be performed with no big concerns about the computational time. However, in the past, computational performance was a major concern, and simplifications were needed to guarantee acceptable computing times. Hence, some models developed with the aim to reduce the computational demands of the dynamic simulations are briefly shown in the following lines.
A simple but effective wear model was proposed by Fries and Dávila [27], who developed a computer programme able to consider the evolution of the wheel profile during the dynamic simulation. The input of the model was the PSD ( power spectrum density) of the rail alignment, and the dynamic problem of a 9 DOF (degree of freedom) quasilinear freight truck was first solved in the frequency domain and then transformed into a time series. Then, for each point in the time series, the contact patch location and the contact patch main quantities were calculated, so that wear could be computed using a proper wear law. Finally, the worn profile could be smoothed and updated for the following dynamic analysis. In an upgrade of this solution [46], the authors introduced a nonlinear model to be used in case of detection of the vehicle hunting behaviour. Figure 4a sketches the algorithm of the programme implemented by Fries Fig. 4 Simplified wear models: a wear model proposed by Fries and Dàvila (adapted from [27,46]); b model relying on quasi-static simulations suggested by Azpezetxea et al. [49] model is an example of combination of dynamic simulation and wear estimation, with a profile update strategy allowing to consider the mutual influence of the two phenomena, the main drawback of this tool was that it only accounted for the calculation of wear on tangent track. Despite stating that their model could be used for the prediction of both the wheel and rail worn profiles, the authors only showed results for the wheel wear.
A different solution was suggested by Kalker [47], who introduced a stochastic motion of the wheelset in order not to require a dynamic simulation, thus drastically reducing the computational time. The strategy for wheel wear calculation was based on the concept of a unit wear step, corresponding to 13,300 wheel revolutions, i.e. approximately 40 km for a wheel diameter of 1 m. In each wear step the wheel profile was assumed to be unchanged and to avoid launching longlasting dynamic simulations, Kalker introduced a probability distribution related to the probability of the wheelset lateral shift to fall into a specific range of lateral positions. The local distribution of friction energy on the contact patches was calculated by means of the FASTSIM routine, and the mass lost was related to the frictional work with a single proportionality constant, as the model only considered tread wear.
Another method to reduce the computational demands connected to dynamic simulations was proposed in 2005 by Magel et al. [48], who proposed a tool named PUMMEL. The authors computed the wear of a newly designed wheel profile on an 800 km run performing quasi-static calculations instead of dynamic simulations and without including a strategy for profile smoothing and update.
An approach based on quasi-static calculation was also recently proposed by Azpezetxea et al. [49]. In fact, although nowadays multibody dynamic simulations can be effectively performed on common computers, the total computational time can still be huge if many simulations must be performed. The strategy developed by Azpezetxea et al. [49], sketched in Fig. 4b, also included an algorithm for a deterministic identification of the so-called characteristics points of a selected route, based on the values of curve radius and unbalanced acceleration. Once these points have been determined, quasi-static analyses can be performed instead of dynamic simulations, so that the total computational time required is extremely reduced. As noticeable from Fig. 4b, this strategy is clearly an upgrade with respect to the approach proposed by Magel et al. [48], as at the end of each wear step the wheel profile is updated and a new set of quasi-static simulations is performed introducing the computed worn profiles until the desired mileage is reached and the computation is stopped.
Finally, Han and Zhang [50] proposed a different method based on an extensive experimental collection of wheel worn profiles, as they expressed the wear depth as a polynomial function of the wheel mileage and lateral coordinate. However, such a method can only be applied to the specific vehicle and track considered in the analysis, and it also requires a huge amount of experimental data for a proper tuning of the polynomial function.

Models performing a single dynamic simulation and a single wear calculation
The strategy of performing a single dynamic simulation followed by just one wear calculation, see Fig. 1a, can be misleading, as the wheel and rail profiles are kept constant during the whole dynamic simulation. This can eventually turn into the possible formation of wear depth spikes, which can be difficult to manage and smooth. Nevertheless, this approach was recently used by Wang et al. [51], who calculated the wear of the wheels of a high speed train exploiting the Simpack wear module, which does not allow to run new simulations with the worn profiles unless an external code for launching the simulations is developed. The main advantage of this approach is that the user does not need to develop a post-processing routine for wear calculation, as the commercial software provides a user-friendly routine for the definition of wear multipliers, wear coefficients and for the selection of the wear law. However, the main drawback of this solution is that the mutual interaction between wheel profiles and vehicle dynamics is neglected.
This approach is thus mainly introduced when short simulations are performed to assess the influence of specific properties on the computed wear. Biao Li et al. [52], from Singapore, studied the influence of vehicle speed, curve radius and superelevation on the rail wear. The proposed model was built using the commercial software UM, introducing a MB model of a passenger metro vehicle running on a 600 m-long track and computing the rail wear by means of the UM wear module using the modified expression of the Archard law, see Eq. (21). Researchers from China academy of railway science [53] recently proposed a wear model implemented in Matlab, including a MB model of a heavy-haul vehicle and based on a local application of the Archard wear law, to assess the influence of rail cant on the rail profile wear in both straight and curved sections.
Scholars from Central Queensland developed an interesting method to evaluate the wear intensity at the wheel-rail interface according to the value of the wear number Tγ, based on the interaction between a longitudinal train dynamics (LTD) module and a GENSYS locomotive-track MB model, also including the locomotive traction control. In a preliminary work, the authors compared the wear caused by two different locomotives featuring different values of the maximum traction torque running on a short curved track [54]. In this work, the dynamic simulation was performed in post-processing after running the LTD code, see Fig. 5a. An upgrade was then proposed allowing to co-simulate the LTD and the dynamics of two locomotives, with the new model also relying on parallelization to reduce the computational times, implementing the strategy sketched in Fig. 5b [55]. The upgraded model focused on the locomotive wheel wear and estimated the damages due to RCF. A similar strategy, based on the recorded values of the Tγ wear number to predict wear intensity and damages, was also used by Bracciali and Megna [56,57], who launched their simulations in the VI-Rail commercial software package. Nevertheless, although these models allow to estimate the wear intensity for the rails and wheels, they are currently not provided with a module able to compute the worn depth.
The models described above did not rely on co-simulation between the dynamic module and the wear calculation model; however some recent models were developed exploiting a co-simulation between the MB vehicle-track model and the wear computation module. Zhang et al. [58] proposed a tool based on the co-simulation between a locomotive-track Simpack model and a Simulink code for the calculation of the wheel wear depth, see Fig. 5c. The Simulink wear module implemented a global application of Archard law. The co-simulation interface allowed to send the global contact patch parameters from the Simpack model to the Simulink wear model online in each time step of the simulation, thus computing the wear depth while also running the dynamic simulation, with an increase in the computational speed with respect to the traditional solution of using the Simpack built-in wear module. The authors performed different simulations on short curved tracks and then also carried out a simulation on a track representative of the typical Chinese lines to compute the wear corresponding to a total mileage of 10,000 km. The authors recognized that the model should be upgraded in order to perform new simulations with updated profiles.

Models performing a single dynamic simulation and an iteration of profile updates
The strategy discussed in the previous section can be upgraded by performing an iteration of profile smoothing and update operations until the desired mileage or number of vehicle passages is reached, as shown in Fig. 1b. The hypothesis behind this approach is that the contact point position and the values of the contact forces and global creepages are little affected by changes in the wheel and rail transversal profiles. At each iteration of the wear depth calculation, an optional wear multiplier can be introduced to magnify the amount of wear depth along the profile transversal direction, so as to simulate a longer run with respect to the simulated one. Clearly, the introduction of a wear multiplier allows to reduce the total number of simulations, but it tends to gather wear always in the same profile positions, thus possibly generating wear depth spikes. Therefore, a proper value of the wear multiplier must be selected, and the wear depth distribution as well as the worn profiles need to be smoothed to avoid numerical issues. This strategy was the one initially used in the wear modules implemented in commercial MBS; however nowadays commercial wear modules are far more advanced.
One of the authors of the present review paper implemented such strategy with the aim to drastically reduce the computational times [59], developing the algorithm sketched in Fig. 6a MB software package Simpack, and the results underwent a statistical analysis to detect a small number of classes N c , under the hypothesis that the points that fall into a specific class give a similar contribution to the total wear. Several statistical procedures were investigated, and the one that proved to be the most accurate was based on a division into classes depending on the wheelset lateral shift, with all classes having the same number of elements (therefore, the wheelset lateral shift was not discretized with a constant step). After obtaining the statistical representation of the outputs computed by the MB software, the wheel profile was updated accounting for the statistical "weight" of each class. The calculation of the wheel worn profiles was carried out using an iterative strategy, thus performing an iteration of profile updates (in Fig. 6a, the total number of iterations is equal to N). At each iteration, the statistical values obtained from the original MB dynamic simulation were left unchanged; however a post-processing contact module was launched again introducing the worn wheel profile computed at the previous iteration, thus considering that the changes of the profile in transversal direction have a limited effect on the contact point position and on the global forces and creepages. The research group from China Academy of railway science still relied on the approach of launching a single dynamic simulation followed by an iteration of profile updates with the aim to accurately simulate the wheel-rail contact rather than to simply reduce the computational times, through the implementation of the algorithm depicted in Fig. 6b. The research team performed a co-simulation between MB dynamic simulations run in ADAMS Rail and a detailed local contact model implemented in the FE (finite element) commercial software ABAQUS [60,61]. The dynamic simulation was run just once in the MB environment, and then an iteration of wear steps took place in the post processing stage. At each iteration of the wear loop, the outputs of the global contact model, computed in the MB simulation, were loaded into a detailed 3D FE model for the steady-state solution of the wheel-rail rolling-contact based on the arbitrary Lagrangian-Eulerian approach. The FE model considered nonlinear material behaviour for wheel and rail steels and also took into account the dependency of

3
friction coefficient on the sliding speed. The wear depth was computed in each node of the FE model applying the KTH Archard wear law, and then the wear depths were spread along the wheel circumferential direction, to obtain a distribution of wear depth along the transversal profile. The profile was updated when the maximum wear depth exceeded 0.1 mm, also applying a cubic spline interpolation to smooth the profile and avoid numerical short-wave concavities. A new 2D mesh of the wheel profile was generated by means of a Laplacian smoothing technique with area weights, and then the 3D mesh was obtained by applying an axis symmetry condition to the 2D profile mesh.
Although the strategy sketched in Fig. 1b does not fully consider the mutual dependency of profile wear and vehicle-track dynamics, it is obvious that the researchers from China academy of railway science opted for this solution as an iteration of dynamic simulations and detailed FE analyses for wear calculation would require an extremely long computational time. Some years have passed since the development of the above model, and probably nowadays the computational speed might be increased by means of co-simulation and parallel computing techniques allowing to perform the wear calculation while running the MB simulation. Therefore, the approach could be upgraded in order to run new dynamic simulations with the worn profile computed at the previous iteration.
In fact, Bosso et al. [62] recently compared the results obtained both when considering and neglecting the influence of the vehicle-track coupled dynamics on the wear evaluation, and they showed that the worn profiles computed when leaving unchanged the outputs of a single dynamic simulation strongly differ from the ones computed when the dynamic simulation is launched at each iteration with worn profiles. Therefore, thanks to the computational power of modern computers, this approach is almost abandoned in favour of the strategy described in the next subsection, which is based on an iteration of dynamic simulations.

Models based on an iteration of dynamic simulations
Nowadays, the typical approach used for the determination of wheel and rail transversal profile wear relies on an iteration of dynamic simulations, each one corresponding to a wear step, in which the wheel and rail profiles are unchanged, see Fig. 1c.
The worn profile generated from the dynamic simulation is computed and used as the initial profile of the next iteration. The process is repeated until a desired mileage or number of vehicle passages is reached, or until the profile characteristic dimensions exceed the prescribed limit and reprofiling is necessary. This approach can rely on the application of a wear multiplier in the wear evaluation module, too, so that the computational time and the total number of iterations can be strongly reduced. Several models witnessed in the literature fit into this category; however, they can still be distinguished into models in which both the dynamic simulation and the wear calculation are performed in the same environment and models that rely on co-simulation, i.e. the dynamic model and the wear module run on two separate software platforms. Before each journey, a simplified contact code was run in the pre-processor to pre-calculate the contact patch sizes and positions according to the Hertzian theory. At the end of each transient analysis, the material loss was calculated with a dedicated routine in which the BRR equations were implemented, see Eqs. (6)-(8), starting from the outputs of the MBS, namely the position of the contact point, the contact patch size and the global forces and creepages. Wear was computed globally, due to the inefficiencies of computers in the early 90 s, and the worn material was spread parabolically within the contact patch. Two separate wear multipliers were used for the straight and track sections, to reduce the computational times, and their values were such that each loop described the wear due to a total travelling distance of 1100 km.
A step forward was taken by the research team from Budapest University, led by Professor Zobory, who proposed a model for the estimation of wheel and rail profile wear at the same time, based on the local computation of the removed material in each cell of the contact patch [36,64]. The model calculated distinct worn profiles for the left and right rail for each track section, as track sections were classified according to the curve radius. The Hungarian research group recognized that while the wear process in the wheel evolves as a function of the travelled distance, the rail wear process evolves as a function of the total number of wheelset passages. Therefore, to compute the worn profiles of rail and wheel simultaneously, the amount of worn material was magnified so that a single computation of wheel wear corresponded to a running distance of approximately 1000 km, while the computation of rail wear was magnified to represent approximately 40,000 wheelset passages on each track section. Hence, after each call to the multibody model performing the vehicle-track coupled dynamics simulation, a single worn profile was calculated for each wheel, with the worn material magnified to represent a travelled distance of 1000 km, while as many worn rail profiles as the number of track sections in the simulated track were obtained, with the rail worn material magnified to represent 40,000 wheelset passages. The computed worn profiles were then used as inputs for the following dynamic simulation and the wear computation was iterated until a desired mileage and number of wheelset passages was simulated. Clearly, the profile update operation required the application of a smoothing strategy, which in the authors' opinion is the most interesting aspect of the model proposed by Zobory's research group, as the smoothing procedure was applied in two stages only in case of a detection of an incompatibility between the worn wheel and rail profiles. The first stage performed a "physical" smoothing, i.e. a smoothing of the profiles performed to guarantee the compatibility of the wheel and rail profiles. The second stage was a "mathematical" smoothing based on a spline interpolation of the profile points altered  Fig. 7 a Wear model proposed by BRR researchers (redrawn from [63]); b model for wheel and rail profile wear evaluation adopted by Zobory's group; c model by Zobory's group for computation of metro wheel profile wear in the previous "physical" smoothing step. The smoothing procedure was carried out iteratively, and the process was stopped when the wheel and rail profiles after the mathematical smoothing were still physically compatible. Since at each iteration a single wheel worn profile was computed while the number of calculated worn rail profiles was equal to the number of track section, the smoothing procedure only modified the wheel worn profile, while keeping the rail worn profile unchanged. However, the original works by Zobory and Szabo lack details about this computation. In the authors' opinion, the strategy followed to smooth the worn profiles could be similar to the one depicted in Fig. 7b. According to the authors' hypothesis, as the final wheel worn profile should be compatible with all the worn rail profiles computed in the different track sections, the wheel smoothing process has to be carried out until compatibility is ensured on all rail profiles computed in each track section. Please note that in Fig. 7b N ts indicates the number of track sections, which corresponds to the number of computed rail profiles.
The strategy used by Zobory's group to calculate the wheel and rail profile wear within the same computation framework was based on the assumption of a "deterministic" operation on the considered railway line, i.e. under the assumption that a single vehicle runs on the line. Nevertheless, the Hungarian scholars mentioned that in case the wear of wheel and rail on a whole network had to be computed, considering many different vehicles with different wheel profiles running on many lines, a stochastic approach had to be developed, in order to take into consideration the probability that a certain vehicle with a certain wheel profile runs on a specific track section.
In a further work the research group from Budapest University proposed a method to compute the traction/braking forces and the speed profile of a metro vehicle using a LTD code, and to feed this data into the MB simulation [65]. In this work, the authors were mainly focused on the determination of the worn profiles of the wheels mounted on metro vehicles; therefore, a slightly different strategy was used to account for the possible contact between worn wheels and worn rails, see Fig. 7c. First, a simulation was carried out considering new rails and new wheels, and wheel profiles at different wear stages were computed. Then, a simulation was carried out with wheel profiles at different wear stages to compute worn rail profiles along the track. Finally, the evolution of the wheel profile wear could be evaluated using the worn rail profiles obtained from the previous preliminary simulation.
Since the early 2000s much work has been carried out at the Swedish Royal Institute of Technology, USFD and PoliMi to develop numerical models for the estimation of wear at the wheel-rail contact, mainly focusing on the calculation of the wheel worn profile. The models from these research institutes all shared the same main general architecture sketched in Fig. 1c, relying on numerical simulations performed in an MB environment and on a post-processing calculation of the worn material and of the worn profiles. Both modules were implemented in the same computational environment thanks to the development of user-written wear routines, as these models did not rely on co-simulation strategies. To compute the wear depth distribution locally in each element of the contact patch for each time step, all these models introduced a post-processing local contact model which performed the calculation of the slip velocity and tangential pressure distribution. Such local contact model may differ with respect to the tangential contact model used in the MB code for the computation of the tangential forces. In fact, when the wear depth is calculated locally in each element of the contact patch grid, it is essential to obtain an accurate estimate of the slip velocity and traction pressure distributions, while in the MB simulation a fast method for the computation of the tangential forces is essential to speed up the calculation.
Focusing on the work carried out at KTH, Jendel [38,66] implemented a dynamic simulation module and a wear estimation module in the commercial MB software GEN-SYS, providing the wear module with a local contact model based on the FASTSIM algorithm. The wear depth was calculated locally in each element of the contact patch according to the Archard wear law, reducing the wear coefficients obtained from the KTH map in straight sections and low rails to account for the effects of natural lubrication and in high rails to consider the effects of artificial lubrication. Reduction factors equal to 10-15 were suggested for high rails, while values in the range of 2-5 were recommended for low rails. Each wear step corresponded to a maximum wear depth of 0.1 mm and before the profile update, a cubic spline interpolation algorithm was applied to smooth first the wear distribution and then the worn profile.
The tool was upgraded by Enblom and Berg who introduced in the MB model the forces due to disc braking and modified the local contact model considering the contribution of the elastic slip in the calculation of the sliding velocity in the FASTSIM local contact model [67]. In this work the authors modified the update criterion, so that each wear step corresponded to a maximum wear depth of 0.1 mm or to a total running distance of 1500 km. The model proposed by Enblom is still employed as witnessed in the work by Li et al. [68], who used Simpack to run the dynamic simulations and implemented Jendel's strategy for the reduction of the KTH map wear coefficients to calculate the wear of high-speed train wheels running on Chinese lines.
The KTH numerical tool was also provided with a postprocessing routine for the estimation of the probability of crack initiation and propagation [69], also implementing this version of the model in Simpack rather than in GENSYS in a 421 Simulation of wheel and rail profile wear: a review of numerical models collaboration between KTH and Bombardier [70]. Initially, the tool could only estimate a probability of crack initiation, but in following works the RCF model was upgraded so that it was able to compute the crack length [71]. Similar models able to compute both wear and RCF were used by Hosseini et al. for the evaluation of RCF and wear evolution on ironore locomotive wheels [72,73].
The KTH wear model without the additional RCF routine was also used by Enblom and Berg to calculate the wear evolution of rail profiles [74]. Obviously, the wear of rail profiles depends on the total number of axle passages on a specific track section and different types of vehicles with different worn profiles and running in different operating conditions can run on the same track section. The simulation input dataset was therefore prepared according to a statistical procedure which implemented a hierarchical strategy for the definition of the operating conditions probability. At the highest level was the definition of the probability of a certain vehicle to run on the considered track section, while the second level accounted for the operation probability, with each operation defined in terms of vehicle speed and braking deceleration. The lowest level finally classified the occurrence probability of specific contact conditions in terms of friction coefficient, wear coefficients, contamination and wheel profile. The hierarchical strategy for the definition of the occurrence probability is sketched in Fig. 8a. The simulation strategy roughly followed the same approach used for the determination of the worn wheel profile. However, different worn rail profiles were computed for different track sections with different curve radius, but for each section the variation of the rail shape in longitudinal direction was neglected. More in detail, the authors considered four non-lubricated curves and two lubricated curves with radius in the range 303-802 m. Another difference with respect to the wheel wear model concerned the strategy for the profile update, which was performed so that each wear step corresponded to either a maximum wear depth of 0.1 mm or to a maximum number of axle passages, set equal to 1000 for non-lubricated curves and equal to 3000 for lubricated curves.
While work was being carried out at KTH to develop reliable models for wheel and rail uniform wear calculations, much effort was produced by researchers from USFD and Politecnico di Milano (PoliMi) to implement similar models. Scholars from PoliMi proposed a wheel wear prediction model based on an in-house built MB code and then implemented a post-processing routine provided with a local contact module which computed the distribution of traction pressure and sliding velocities according to the Kik-Piotroski or CONTACT algorithms, with the Kik-Piotroski method clearly proving to be a faster solution [75]. To speed up the dynamic simulation, the tangential contact forces in the MB code were computed according to the heuristic model by Shen et al. [76], which differed from the algorithm used in the post-processing local contact model (CONTACT/Kik-Piotrowski), as the heuristic model computes the tangential forces directly from the longitudinal, lateral and spin creepages, without requiring the determination of the slip velocity distribution on the contact patch. In the post-processing routine, the wear depth was calculated as proportional to the specific power dissipated at each element of the contact patch, see Eq. (25), in which k A is a proportionality constant. The whole 3D wheel geometry was discretized into regular squares, so that at each time step of the simulation, the wear depth could be assigned to the correct circumferential position on the wheel surface. Therefore, this method could theoretically be used to compute the wear in transversal and longitudinal direction at the same time. Nevertheless, for the sake of computational efficiency, the authors only focused on the wear of the wheel transversal profile. The profile was then updated and smoothed and used as initial profile in the following dynamic simulation. Fig. 8 a Hierarchical occurrence probability classification proposed by KTH researchers to define simulation input set for rail wear prediction (redrawn from [74]) and b wheel discretization adopted by USFD group (redrawn from [77,78]) Contextually, the research group from USFD implemented a similar model in the commercial MB software ADAMS Rail [77,78]. In this case, the local contact algorithm was similar to the FASTSIM method, and the wheel circumference was discretized into longitudinal strips, see Fig. 8b, so that at each time step, the position where the wear depth, calculated using the USFD wear model, was to be applied could be easily identified. However, since the authors were interested in the calculation of the wheel transversal profile, the wear depths of each longitudinal strip were summed and uniformly distributed along the circumferential direction in order to get a wear depth distribution along the lateral wheel profile coordinate.
A collaboration between PoliMi and USFD led to the development of a new model, named ProfCon, exploiting the same general idea of the tools described above [79]. The dynamic simulation was performed using the PoliMi in-house MB code, solving the normal contact problem by means of a multi-Hertzian approach and the tangential contact according to the heuristic method by Shen et al. The local contact model implemented in the wear module was instead based on a simplified version of FATSIM neglecting the spin creepage, and the wear depth was computed locally using the USFD wear law. Unlike the PoliMi and USFD models, the ProfCon tool did not rely on a discretization of the wheel surface along the running direction, as at each time step the wear depth was uniformly spread along the circumferential direction. The cumulative wear depth distribution was smoothed by applying a moving average and then the worn profile was interpolated with a cubic spline before starting the new dynamic iteration. Each wear step in the model corresponded to a maximum wear depth of 0.1 mm.
All the models described in the lines above in this subsection rely on a post-processing evaluation of the worn material and of the worn profiles. Nevertheless, other models exist in the literature that compute the worn material at each time step of the dynamic simulation, but that do not perform an online update of the worn profiles. The wear depth is therefore cumulated at each time step, and the profile update is performed at discrete intervals when a proper profile update criterion is satisfied.
Many models characterized by this feature were developed at Southwest Jiaotong University (SWJTU), solving the vehicle-track coupled dynamics using an in-house MB code. This choice was preferred as it allowed to implement a modified version of the CONTACT algorithm in the local wheel-rail rolling contact module. The general structure of the tool proposed by SWJTU researchers is shown in Fig. 9a, in which the reader can notice the upgraded CON-TACT algorithm, which receives as input the wear depth calculated in previous steps. The wear depth is considered when calculating the normal distance between the wheel and rail profiles, which is essential to determine the normal contact force. At first, the research group proposed a tool for the calculation of the rail wear in curved track sections, modelling a half vehicle running at constant speed [80,81]. The rails were treated as Euler beams considering their vertical, lateral and torsional flexibility by means of modal superposition techniques. The MB code solved the normal problem by means of nonlinear Hertzian springs, while the tangential forces were computed using the global model by Shen et al. The outputs of the MB code were then fed to a wear module comprising a local contact module based on the modified version of CONTACT, able to account for the wear depth calculated in previous steps. The wear depth was estimated locally through a dedicated routine, relating the wear depth in each element of the contact patch to the frictional work density, using the expression suggested by Igeland and Ilias [82]. In these preliminary works, the rail profile was not updated, so these models should belong to the category sketched in Fig. 1a and surveyed in Sect. 4.1. Nevertheless, they are presented here to show in the same subsection the models developed at SWJTU, which share the same main features and the peculiarity of computing the wear depth at each time step. Please note that in the SWJTU models, the rail worn profile was calculated in a specific location of a curved track section.
Later on, the SWJTU research group developed two separate models for wheel [83] and rail [84] wear prediction including dynamic models of a whole vehicle rather than a half vehicle. The two models shared the same main features, namely a dynamic simulation carried out using the in-house dynamic module and a similar strategy for the profile update, based on the total travelled distance or maximum wear depth for the wheel wear model and the total amount of wheel passages or maximum wear depth for the rail wear tool.
In the wheel wear model proposed in [83], the vehicle-track coupled dynamics module took into consideration the track flexibility and the rails were modelled as hinged Euler beams, using a modal superposition technique to consider the rail vertical, lateral and twisting flexibility. In contrast with the previous rail wear model, the wear depth was calculated using the strategy adopted by the KTH research team instead of the Igeland and Ilias's expression, and the authors still assumed that the contact patch was constant for a whole wheel revolution, so that the wear depth could be spread along the circumferential direction. To avoid numerical noise and spikes, a smoothing procedure based on cubic spline interpolation and Friedman's supersmoother [85] was applied to both the wear depth distribution and the worn profile.
As already mentioned, the rail wear model presented in [84] shared the same main features with the wheel wear module; however in the rail wear tool the rails were modelled as Timoshenko beams, and a discrete moving sleeper support was introduced. To account for the effect of lubrication, the authors modified the KTH Archard coefficients as suggested by Jendel, and also used two separate values of friction coefficient on the rail head and on the rail gauge. In this work, the authors performed simulations with wheel profiles at different wear stages to account for the possibility of different wheel worn profiles to run over the same track section.
Recently, Shi et al. [86] modified the SWJTU wheel wear tool using a global contact module developed by UniFi [87] for the detection of the wheel-rail contact points based on the trace curve method, which considers the wheelset yaw. The local wheel-rail contact module relied on the FAST-SIM algorithm, and the wheel wear depth was computed locally using the USFD wear law. The procedure for profile smoothing and update was the one suggested by scholars from UniFi, which will be better described in the following subsection.
The main advantage of the models developed at SWJTU is that the implementation of the wear model in an in-house code allows to easily modify the CONTACT algorithm so that the wear depth calculated at previous steps can be included in the solution of the wheel-rail rolling contact problem; however it is clear that using CONTACT can Another model able to account for the influence of the rail flexibility on the calculation of the position of the wheel-rail contact points and forces was proposed by Aceituno et al. [88], who implemented in an in-house MB code the finite element floating frame of reference method to account for the rail flexibility. In this model, sketched in Fig. 9b, the wheel wear depth was calculated at each time step; however the profile was updated only when a limit wear criterion was reached. To save computational time, as the modelling of the rail flexibility requires a huge effort, the wear was computed globally and then spread as a volume loss rate on the contact patch assuming the same trend as for the contact pressure: where V worn is the rate of worn volume; D is the wheel diameter; ΔS w is the area removed by wear, calculated according to the BRR law; p z is the normal pressure; V is the running speed; A e is the area of an element in the contact patch grid.
The works shown in the previous lines dealt with wear tools that either were implemented in in-house codes or that only relied on commercial MB codes for the vehicle-track dynamic simulations but that needed user routines for the calculation of the profile wear. Nowadays, commercial MB software are provided with much more powerful wear tools, which do not require the users to implement their own codes and routines for wear calculation, profile smoothing and update. Wang and Gao [89] recently showed an example of wheel wear evaluation based on the wear module available in UM. The UM wear module can update the wheel profiles to be used in the consecutive wear steps and it is also able to run several simulations in parallel and to apply weighing factors to account for different vehicle configurations, track layouts and speeds. The UM wear module was also exploited by Shadfar and Molatefi [90], to assess the influence of different parameters on the wheel wear evolution, as well as by Pires et al. [91] in the context of the design of new wheel profiles. On the other hand, Smetanka et al. [92] opted for the commercial software Simpack to run the wheel wear simulations; however the Simpack wear module is not able to automatically update the wheel profile by means of a simple graphic user interface once a wear step is over, so a code able to launch a cascade of Simpack simulations with updated profiles is required. Such code can be implemented using the QtScript language, which is integrated in the Simpack environment.
The strategy of using a commercial MB code provided with a wear evaluation tool has the main advantage of easy implementation, with no need for the users to write their own local contact and wear routines. On the other hand, the users can only set a limited amount of parameters for wear profile smoothing and update as well as for wear depth magnification, thus not having a full control on the computation and on the stability of the calculation.

Numerical tools relying on co-simulation techniques
As previously mentioned, using in-house wear evaluation tools allows to have a total control on the computation, with the chance also to include special features in the dynamic simulation and in the wheel-rail global and local contact modules. On the other hand, implementing a model in a commercial MB package guarantees an easier implementation and ensures the user to use reliable and validated routine. Therefore, a good compromise solution is the cosimulation between a commercial MB software, solving the vehicle dynamic simulation, and an in-house routine developed for the evaluation of the wheel wear and possibly of the contact forces, as shown in Fig. 2b and 2c. This solution, which is promoted by modern computer architectures, allows to speed up the computation, especially in case parallel computing techniques are implemented. The typical strategy followed by the models presented in the following lines is still the one depicted in Fig. 1c, with the update of the worn profiles performed at discrete steps at the end of each dynamic simulation. Once again, in case wear is computed locally in each cell of the contact patch, the wear module implemented in the external routine must also include a local contact model for the evaluation of the slip velocity and tangential pressure distributions on the contact patch at each time step.
A simple application of this approach was proposed by Zhu et al. [93], who slightly modified the KTH strategy performing the dynamic simulations in Simpack and then computed the wheel wear with a dedicated Matlab routine, which also determined the wheel worn profile shape.
An interesting example of co-simulation approach was proposed in the context of the AWARE transfer of knowledge project, as a collaboration among ALSTOM Ferroviaria, USFD, Technical University of Lisbon and University of Zilina [43,44,[94][95][96]. The tool, sketched in Fig. 10, included a dynamic module for the simulation of the vehicle dynamics, implemented in the commercial software VAM-PIRE, and a Matlab kernel dedicated to the computation of wear and to the profile smoothing and update for the next wear step. To account for the rail lubrication in tight curves, the authors defined different friction coefficient values at the wheel tread-rail head and wheel flange-rail gauge interfaces.
The tool was conceived with a modular approach, as the users could select either a local or global application as well as the desired wear law, chosen among the BRR, USFD and KTH laws. Clearly, in case wear was computed locally, the wear computational tool included a local contact algorithm (FASTSIM) to determine the distribution of local stresses and creepages on the contact patch. Furthermore, the user could also select two different profile smoothing strategies, namely a localized approach and a profile intersection technique. The first approach is based on the cumulation of all wear depths for the same wheel-rail relative displacement values, so an intensive smoothing procedure is subsequently needed to spread the total material removal on a larger region of the profile. On the other hand, with the intersection profile approach, the wheel and rail profiles are geometrically intersected for each considered value of the lateral shift in order to detect the profile transversal coordinates affected by the material removal, and the total wear is spread geometrically according to the profile intersection area shape.
The main drawback of the model presented above is that it cannot implement a parallel computing technique, as the Matlab routine is called in the post-processing stage after completing the dynamic simulation in VAMPIRE.
The combination of a VAMPIRE dynamic simulation with a MATLAB routine was also exploited by Bevan et al. [97] in the frame of the development of a model for the calculation of the wheel damage considering both wear and RCF. The Matlab routine was implemented to manage the pre-processing stage, for the determination of a statistical track representative of the real mission of the wheelset, and the post-processing stage, thus calculating the wheel wear, according to the KTH strategy, and the RCF damage index.
The research team from Università di Firenze (UniFi) proposed a model similar to the one developed in the context of the AWARE project; however in this case the contact model used in the dynamic simulation was implemented in an inhouse code implementing a rigid contact method allowing to compute any number of wheel-rail contact points by means of a semi-analytical method. The global contact model was implemented in a C/C + + external routine communicating online with SIMPACK, and it was demanded to detect the contact points and to compute the normal and tangential forces. The normal problem was solved according to the Hertzian theory, while the tangential problem was solved by means of the fast heuristic method by Shen et al. Initially, the model only aimed at computing the wear of the wheel profile, assuming a constant rail profile. The wear module was implemented in MATLAB, and it included a local contact module, based on the FASTSIM algorithm, and a wear evaluation module, based on the local application of the USFD wear law. Figure 11a sketches the tool developed by UniFi researchers for the computation of the wheel worn profile.
In preliminary works, at the end of each wear step, a single worn wheel profile was calculated as the average worn profile of all wheels, and it was smoothed with a moving average filter and then fed back into the dynamic module. The authors explored two scaling strategies for wear depth magnification, which are essential to reduce the computational speed, namely an approach based on a constant mileage for each wear step [98,99], and a second strategy relying on an adaptive wear step [100][101][102], so that each step corresponded to a prescribed maximum wear depth. To further speedup the computation, the research group from UniFi opted for a statistical representation of the reference railway track, introducing proper weighing coefficients to account for the layout of the real track.
The model was then also extended to compute both the wheel and rail profile wear, using a smart algorithm to consider the fact that wheel profile at any state of the wear evolution process can run on any kind of worn rail profile. Once again, in some works the authors computed a single average wheel profile for each wear step using a constant mileage step [103,104], while in [105] an adaptive step for both the wheel and rail profiles was adopted. On the other hand, in [106] an average profile was calculated for each wheelset, under the assumption of a full symmetry of the track and of the consequent wear of the right and left wheel of each wheelset. Anyhow, it is clear that the wear module could be easily modified so that a worn profile could be computed for all wheels. Furthermore, it is interesting to notice that in [105], the authors benchmarked their tool against the results computed by the Simpack wear module, which applies the Krause-Poll wear law globally, launching the dynamic simulation and wear calculation in the commercial software, and then using an external routine for the wheel profile update and smoothing operations. This strategy of using the external wear routine to only smooth and update the wheel profile was recently investigated by two authors of the present paper [14], to rely on the fast computation ensured by Simpack while also improving the numerical stability of the worn profile computation.
In recent years, some upgrades were proposed to improve the capability of the UniFi tool. In [107], the research team from UniFi upgraded the contact model used in the Simpack dynamic simulation, replacing the global contact model with a local model implementing the non-Hertzian Kik-Piotrowski method to solve the normal problem and a modified FASTSIM method, applied to non-elliptic contact areas, to deal with the tangential problem. On the other hand, in [108], the original wheel and rail combined wear prediction tool was also provided with a RCF module.
A great point of merit of the works by scholars from UniFi is the detailed description of the strategies used for the profile smoothing and for the calculation of both rail and wheel wear in the same tool, although the two wear processes take place in very different time scales.
Since the wheel wear process is quite faster with respect to the rail wear, any possible wheel worn profile can run on any worn rail profile. The authors thus proposed an interesting solution to deal with this issue, and the algorithm adopted is briefly sketched in Fig. 11b. First, the authors decided to split the total vehicle travelled distance, km tot , in a number n sw of wheel steps, so that each wheel step corresponds to a travelled distance equal to km tot /n sw . Similarly,   Fig. 11 UNiFi wear model: a general structure of wheel wear estimation tool; b algorithm for update of wheel and rail profiles proposed by UniFi researchers the total number of vehicles travelling on the track N tot was divided into n sr rail steps. To better understand the strategy for the calculation of wheel and rail worn profiles at each stage of the wheel and rail wear processes, a simplified hypothesis is made in the following lines, which is that at each wear step, a single average profile for all wheels is obtained, as the UniFi reserchers actually did at first. At the beginning of the simulation, the wheel and the rail have each their nominal profile. The first step is the dynamic simulation of the vehicle track interaction for the n sw wheel steps, which is performed in a "wheel wear loop" comprising n sw iterations. This means that at the ith iteration of such loop, the wheel profile is calculated and updated, while the rail profile is assumed to be constant for the entire simulated trip. Hence, a number of dynamical simulations equal to n sw is performed, and also n sw calls to the wear module are made.
At the end of each iteration of the wheel wear loop, the worn wheel profile is saved to a memory stack. Once the loop is completed, the rail profile can be updated according to the following approach. For each stage of the wear evolution process, a rail profile is obtained, thus obtaining n sw worn rail profiles, which are calculated in a "rail wear loop". Each rail profile is saved to a memory stack and once the n sw worn rail profiles are computed, they are averaged and a single worn profile is obtained. The rail wear loop is then completed. Please note that at this stage of rail profile calculation, only n sw calls to the wear routine are required, with no need for extra dynamic simulations.
Then, a new set of n sw dynamic simulations can be performed using a constant rail profile, which now is the previously determined worn rail profile at the end of the first rail step. Once again, n sw wheel worn profiles are obtained, starting from the nominal unworn wheel profile. Therefore, a second averaged rail profile can be easily determined with n sw calls to the Matlab wear estimation routine. The second rail step is thus completed. The rail steps are repeated n sr times, so that at the end the final results is an evolution of the wheel profile running on each rail worn profile. From a computational point of view, the model requires a number of dynamic simulation equal to n sw × n sr , and a total number of calls to the wear estimation module of 2 × n sr × n sw .
For what concerns the approach for the profile smoothing operation, the UniFi model is based on the following 7 steps: 1. Spread of the wear rates computed at each contact point on each wheel-rail pair for each simulation in the rail longitudinal direction and in the wheel circumferential direction. 2. Time integration of the wear rates to obtain the wear depths. 3. Sum of the contribution of all contact points to the overall wear depths.
4. Weighted average to consider the real layout of the track, and possible averages to obtain a single wheel worn profile or a wheel profile for each wheelset. 5. Wear scaling, which can be performed according to a constant or adaptive strategy, i.e. with a constant travelled distance or with a travelled distance limited to a threshold corresponding to the maximum wear depth allowed. 6. Smoothing of the wheel and rail profiles by means of a moving average. 7. Update of the profiles for the next simulation.
Please note that conversely to other wear models for rail worn profile calculation presented above, such as the models from Zobory's research group, from KTH and from SWJTU, the UniFi model computed a single worn rail profile applied to the whole length of the simulated track, instead of different worn profiles for separate track sections. The resulting worn profile can hence be considered as an "average" rail profile, as the authors defined a "statistical" track and introduced weighing factors to account for the occurrence probability of the simulated track sections on the real line.
The smoothing and update strategy suggested by researchers from UniFi was adopted by Tao et al. [109] in an implementation of a wear model based on the co-simulation between a dynamic simulation performed in Simpack and a wear module written in Matlab. The wear module included a local contact module, in which the normal problem was solved using the Hertzian theory and the tangential problem was solved using the Fastrip algorithm [110]. The computation of the wheel-rail forces in the Simpack MB simulation was instead carried out using the Hertzian theory to solve the normal problem and the FASTSIM algorithm to solve the tangential problem. Wear was calculated by applying the USFD wear law locally in each element of the contact patch grid and the wear depth distribution along the profile transversal coordinate was first processed with a moving average and then smoothed with a cubic spline, which was also applied to the worn profile. The updating criterion selected by the authors was an adaptive criterion on the maximum wear depth related to each wear step. The same strategy based on Simpack dynamic simulations followed by a wear module based on the FaStrip algorithm and on the USFD wear law was recently implemented by Ye et al. [111] in the context of optimizing the suspension of a freight wagon to minimize wear, introducing a Kriging surrogate model (KSM) written in a Matlab routine. The KSM was used to establish a relationship relating the wheel wear to vehicle speed and vertical stiffness of the primary suspension of the Y25 bogie, to reduce the number of simulations to be performed and obtain a quicker convergence to the optimum solution.
A model similar to the one proposed by Tao et al. was also developed by Ding et al. [112], who however used a multi-Hertzian approach to solve the normal problem in the post-processing local contact module, to deal with nonelliptic contact patches, and used the wear law proposed by Zobory to compute the worn material. Furthermore, in this work the profile update criterion corresponded to a maximum wear depth of 0.1 mm and the smoothing strategy was based on a wavelet packet filter, which treated the distribution of wear depth as a noisy signal.
The profile smoothing operations suggested by UniFi scholars were also adopted by Indian researchers from Indian Institute of Technology (IIT), who implemented a wear prediction tool based on the co-simulation between ADAMS VI Rail and a Matlab wear computation routine. This research team proposed both post-processing local contact models relying on a semi-Hertzian approach [113,114], and a simplified model in which the local contact model was based on the FASTSIM algorithm [115]. The worn material was calculated by means of a local application of the Archard law, and the profile update corresponded to a predefined travelled distance. The profile update was preceded by a Hermite polynomial interpolation on the wear depth distribution and a Gaussian filter on the worn profile.
All models presented above in this subsection require the determination of the worn profile which is then fed back to the dynamic simulation module to carry out the next iteration of the wear computation. To tackle this issue and to avoid the need of the profile update operation, Megohe et al. [116] proposed a wear computation tool in which measured worn profiles were subsequently fed to the dynamic simulation block. In fact, the goal of the authors was the estimation of the evolution of the rail worn area in three track sections with different curve radius and the comparison against experimental measurements. The tool relied on a co-simulation between a VI-Rail vehicle-track model and a Matlab routine for rail wear computation. The Matlab routine was launched in the post-processing phase, after the dynamic simulation, and it included a local contact model for the determination of the adhesion and slip areas and for the calculation of the tangential stress distribution by means of the FASTSIM algorithm. The wear depth in each cell of the contact patch was computed according to Archard law.

Models based on an online communication between dynamic simulation and wear module
All the numerical tools described in the previous subsections rely on simplifying hypothesis to deal with the mutual interaction between the vehicle dynamics and the profile wear evolution. An exact solution would clearly require the update of the worn profiles at each time step of the dynamic simulation. Obviously, due to the limited computing power, the models developed in the past could not afford such computation and simplified strategies were introduced, basically under the hypothesis that the wheel and rail profiles could be considered as unchanged for a whole dynamic simulation. This hypothesis was justified by the fact that the wear phenomenon occurs on a large time scale, and the material removal has a little effect on the vehicle dynamics as long as the simulated distance is short enough. Nevertheless, modern computers feature a larger computational power and allow for the implementation of parallel computing and co-simulation techniques; therefore it is not trivial to develop models based on an online interaction between the vehicle dynamic and the wheel and rail profile wear. Recently, researchers from SWJTU upgraded existing wear models to implement an online calculation of the wheel wear, updating the wheel profile after each revolution [117]. Figure 12 shows the algorithm implemented in the online wear estimation tool developed by SWJTU researchers. In Fig. 12, the variable θ represents the angular displacement of the wheel and it is used to control the profile update. The coupled dynamics of a locomotive and railway track were simulated in an in-house MB code. The rails were treated as Timoshenko beams, considering their flexibility and their vertical, lateral and torsional vibrations by means of a modal perform an online estimation of the wheel wear, the contact model used in the MB simulation and in the wear module had to be the same. The wear depth was computed locally in each element of the contact patch using the USFD wear law and then spread along the longitudinal direction, according to the approach proposed by UniFi scholars. Despite computing the wheel wear online and updating the wheel profile at each revolution, the authors still expressed the need for a smoothing operation to avoid numerical noise and short wavelength contributions. The smoothing operation first required a moving average applied to the wear depth distribution, and then a spline was applied to both the wear distribution and to the worn wheel profile. The online model was later upgraded adding a module for the locomotive traction control and solving the tangential problem with a model combining the FATSIM algorithm and Polach's model, able to consider the effects of variable friction coefficient and large creepage values [118]. Furthermore, the contact point detection method was developed using a multi-point contact method based on the tracing line method, which considers the effect of yaw.
Currently, the diffusion of online models is still limited to niche applications, as such models require a huge effort for a proper implementation. Furthermore, it is quite tricky to develop online models performing the dynamic simulation in a commercial MB software, as the synchronization between the commercial MB code and the wear module may be extremely tough to obtain. Therefore, models relying on co-simulation techniques between a MB code for the dynamic simulation and an external routine for wear calculation, with a discrete update of the worn profiles, are still the most widespread solution.

Discussion
The previous section has surveyed in detail the main numerical literature models for the calculation of wheel and wear uniform wear. Table 2 classifies the reviewed models according to the strategy adopted to update the worn profiles, also highlighting the benefits and shortcomings of each approach.
As noticeable from Table 2 and from the previous section, the most used strategy for wheel and rail profile evolution modelling relies on a dynamic simulation followed by a wear module including a local contact module. This approach is followed both when launching the dynamic simulation on a commercial MB software and when running the simulation on in-house codes. Often, the calculation of the wheel-rail contact forces in the dynamic simulation relies on simplified algorithms, while the local contact model implemented in the wear module is more complex. Kaiser et al. [119] pointed out that as the two contact models used in the tool can be inconsistent as they rely on different traction-creepage relationships, usually correction factors need to be introduced in the detailed model.
For instance, Enblom et al. [120] investigated the influence of the semi-Hertzian contact modelling on wear computation using the STRIPES contact algorithm and to obtain the same values of normal forces as computed by the Hertzian solution implemented in GENSYS, the interpenetration value in the local contact module was modified. At the same time, it is obvious that the algorithm selected for wear calculation strongly affects the final outputs and the computational speed. Tao et al. [121] compared different normal contact algorithms (Hertz, STRIPES, ANALYN, Kik-Piotrowski) and found that using a combination of Hertzian solution of the normal problem and FASTSIM method for the tangential problem is still the best compromise between accuracy and speed. Nevertheless, in the authors' opinion, a conformal contact model should be introduced when dealing with worn profiles as due to wear, the wheel and rail profiles tend to become more and more conformal. In such condition, the size of the contact patch increases, and the maximum contact pressure is expected to decrease, thus leading to a decrease in the wear rate.
Despite the huge computational power of modern computers, the reduction of computational times is still a big challenge in wear modelling. In fact, as the wear process evolves on a large time scale, the total simulated time is still high. The most common strategy adopted to deal with this aspect is the definition of wear multipliers which magnify the wear depth, so that a single wear step corresponds to either a prescribed fixed travelled distance/amount of traffic on the line (constant wear step) or to a maximum wear depth along the profile (adaptive step). The latter strategy is clearly adopted to limit the maximum wear depth in a specific point of the profile, to avoid numerical instabilities and wear depth peaks. Nevertheless, magnifying wear assigns the material removal to the same point on the profile, but once the material is removed from a little portion of the profile, that portion is not in contact in the next iteration. To partially deal with this issue, a smoothing procedure is commonly applied to spread the wear depth distribution on a wider section of the profile. However, this strategy requires to properly fix the profile update criterion and tune the smoothing operation to produce a reliable estimation of the worn profiles. Although great attention must be paid to the computational efficiency of the simulation, a direct quantitative comparison of the numerical models described in the previous section cannot be provided. In fact, the numerical tools shown above were Single dynamic simulation and single wear calculation (Fig. 1a) Fast computation of worn profile Suitable to quickly assess the influence of different parameters on wear Simple implementation in most wear tools available in MB software packages Influence between vehicle dynamics and profile shape is neglected If wear for a long route/high tonnage value is determined and wear multipliers are adopted, possible spikes in the worn profiles are produced Numerical accuracy is low Wheel wear models: [48,51] Rail wear models: [52,53] Wheel wear models: Zhang et al. [58] Rail wear models: SWJTU [80,81] Single dynamic simulation and iteration of profile updates (Fig. 1b) Simulated time is reduced The iteration of profile updates mitigates the production of local spikes in the worn profiles Influence between vehicle dynamics and profile shape is neglected Wear multipliers are required due to the low simulated time A proper value of wear multiplier must be identified -Wheel wear models: PoliTo [59], China academy of railway science [60,61] Iteration of dynamic simulations with updated profiles (Fig. 1c Rail wear models: KTH [74], SWJTU [84] Wheel and rail wear models: Budapest University [36,64] Wheel wear models: AWARE [43,44,[94][95][96], UniFi [98][99][100][101][102]107], SWJTU [109,112], IIT [113][114][115]125], Others [93,97,111,116] Wheel and rail wear models: UniFi [103][104][105][106]108] Online computation (Fig. 1d Fig. 1a to 1d require an increasing computational demand, while providing more accurate results. At the same time, the definition of wear multipliers that significantly magnify the predicted wear depth slows down the simulation but better reproduces the real mutual interaction between wear and dynamics. The core of all the wear models presented above is clearly the calculation of the worn material or wear rate as a function of the contact patch parameters through suitable wear laws, that rely on wear coefficients tuned on experimental data. In the authors' belief, there is currently a huge need of new experimental tests to obtain a proper set of wear coefficients. In fact, nowadays literature coefficients are applied outside of the range in which they were determined and even to wheel and rail steel materials that differ from the ones used in the experimental tests. In fact, the BRR and USFD wear laws were obtained for specific pairs of wheel and rail steel materials, so it is not formally correct to leave the coefficients unchanged when dealing with different materials. As a support to this statement, the work from Lewis et al. [122] recently showed that the R7T wheel steel has a different wear behaviour in the catastrophic regime with respect to the R8T wheel steel when paired with the same 900A rail steel. Moreover, the USFD wear law computes the wear rate relative to the wheel only, thus using the same law for the rail steel might be incorrect as the wheel and the rail steels can have different hardness properties. On the other hand, according to Jendel [38] the KTH wear map for the calculation of the Archard wear coefficient was built starting from experimental data collected at KTH and literature data from British researchers, thus obtaining a map that can be considered as material independent. Therefore, the same coefficients can be used for both the rail and the wheel steel, which means that the total worn volume arising at the contact interface can be split equally between the wheel and the rail, as done by the wear tool of the commercial MB software Simpack. Nevertheless, the wheel and rail steels may feature different values of hardness and hence they may wear out at different rates. To deal with this limit, maps that were obtained for the calculation of the wear coefficient considering only the wear of the wheel [21] or rail [28] only could be adopted, but these maps are valid only when the considered wheel steel is paired with the considered rail steel and vice-versa. Furthermore, most of the data used to calibrate the wear coefficients were collected when simulating pure longitudinal creepage conditions, nevertheless higher were rates are usually observed in curves, where the longitudinal creepage is negligible with respect to lateral and spin creepages. Therefore, considering the previous observations, the authors of the present paper strongly recommend that effort in the railway field is driven to obtain new wear coefficients and new wear laws obtained from experimental tests in which also spin and lateral creepage conditions are simulated. The new wear models should possibly be able to split the worn material between the wheel and rail steels proportionally to their hardness.
The literature review has highlighted that the main recent innovations in wear modelling have dealt with more detailed contact algorithms introduced in the dynamic simulator and in the wear module or with a few examples of online update of the worn profile. However, apart from the models proposing an online calculation of the worn profiles, the recent innovations have left the typical approach developed in the 90s unchanged. Thanks to the increase in computer power and to the availability of reliable MB and FE codes, the authors believe that tools intended to predict the wheel damage could also account for the effects of tread braking, which still is the common braking system adopted in freight wagons, but it can also be installed on some metro vehicles and locomotives. In fact, results published in the literature by Walia et al. [123] showed that even for powered wheels, the highest contribution to the wheel tread wear can be due to tread braking, so it is obvious that wheel wear due to tread braking should be considered. The task is clearly not easy as such numerical tool should be provided with a thermal module, because during tread braking the wheel surface heats up to high temperatures which drastically change the material properties. Research in this field is on-going [124,125], but still a complete model is not witnessed in the literature.
Furthermore, all models surveyed in the present paper use simplified strategies to spread the wear depth along the longitudinal direction, nevertheless contact conditions may vary during a wheel revolution due to wheel and rail defects and track irregularities. Therefore, an upgrade in wheel and rail damage modelling would require modelling wheel and rail wear in transversal and longitudinal direction altogether. Although in the past some models were proposed with strategies to mesh the whole wheel surface [75,77], they still only computed wear in transversal direction. For wheel wear prediction tools, meshing the whole wheel surface and then mapping the wear depth to the specific locations at each time steep with no need to spread the wear depth along the circumferential direction is essential to deal with this issue; however changing coordinates from the local discretization of the wheel-rail contact patch to the global mesh may not be an easy task and suitable strategies should be developed.

Conclusions
The present paper surveyed the numerical tools witnessed in the literature for the calculation of wheel and/or rail profile wear in transversal direction. All models need to face the same main challenges, namely the simulation of the interaction between vehicle dynamics and profile wear and the evolution of the profiles during their life. The models were classified according to the strategy adopted for the update of the worn profiles. The main conclusions that can be drawn from the review carried out by the authors are given in the following lines.
1. Nowadays, the most common solution to account for the mutual influence between wear and vehicle-track dynamics relies on performing a dynamic simulation considering constant profiles and then computing in post-processing the worn profiles which are fed to the dynamic simulation module for the next simulation. Previous modelling strategies relied on a single dynamic simulation which was then followed by a single calculation of the worn profiles or by a cascade of profile updates which however did not feed the worn profiles into a new dynamic simulation. The solution of running a single simulation and calculating a single worn profile is nowadays adopted only when general indications about wear process are investigated to keep the computation short. The solution based on a cascade of profile updates with no new calls to the dynamic module is instead limited to niche applications, for instance to model the wheel-rail contact in detail with FE codes. 2. To speed up the computation, the wear depth at the end of the simulation is magnified introducing wear multipliers, so that the wear calculation corresponds to either a prescribed travelled distance/amount of axle passages or to a maximum wear depth along the profile. However, the wear multipliers and the profile update strategy should be chosen carefully to obtain reliable results. Moreover, in this case smoothing of the wear depth distribution and of the worn profiles is typically introduced to avoid numerical instabilities. 3. Running the dynamic simulations in commercial MB codes eases the model implementation and speeds up the computation; however several models were proposed that implement the dynamic module in in-house codes, so that the contact module could be provided with special features, able to account for the wear depth computed in previous steps or to model the wheel or rail flexibilities. 4. The wear modules available in commercial MB software packages are typically used only when fast computational times with no need for extremely accurate results are a primary concern. In fact, these modules usually rely on simplified hypotheses and the user has a limited control on the wear depth calculation and on the profile smoothing operations. Therefore, since the user can only tune a limited amount of parameters, the stability of the computation can become tricky to manage. 5. To increase the computational performance without impacting on the accuracy of the outputs, a common solution that is possible with modern computer architectures is the co-simulation between a dynamic model, launched on a commercial software platform, and a wear module implemented in an external routine, that runs in post-processing or in parallel with the simulation. 6. Online update of the worn profiles is currently limited to applications based on in-house code for the dynamic simulation, mainly because the synchronization between the wear update strategy and the MB dynamic simulation may be tricky.
In the authors' opinion, the following main challenges should be faced in the near future to take a step forward in wear modelling in the railway field.
1. Most models rely on wear coefficients tuned on old experimental campaigns performed on a limited number of wheel and rail steel materials, and these coefficients are also applied outside of the range in which they were obtained and for different materials with respect to the ones considered in the tests. Furthermore, most of the wear coefficients currently adopted were tuned on tests performed when only considering pure longitudinal creepage conditions, and these coefficients are extended also to situations in which a combination of longitudinal, lateral and spin creepages is present. Therefore, new experimental campaigns should be carried out to obtain a wider set of wear coefficients. 2. At the same time, solutions used to deal with lubrication at the wheel flange-rail gauge interface are based on either reducing the friction coefficient or reducing the wear coefficients. However, the application of a lubricant at the wheel-rail interface substantially changes the contact interface from a tribological point of view. Therefore, it is essential to carry out experimental tests to address the phenomenon and adopt contact models able to describe the presence of a lubricant spread at the contact patch. 3. Due to wear, the wheel and rail profiles tend to become more and more conformal, and the wear rate consequently decreases. Therefore, even when running the dynamic simulation in commercial MB codes, contact models able to consider the conformity condition should be introduced, in order not to overestimate the damage of the profiles.
4. Thanks to the computational power of modern computers and to the availability of reliable commercial MB and FE codes, it is the authors' belief that wear models could be upgraded to consider wear in transversal and longitudinal directions at the same time, while more difficulties may arise in accounting for the wear due to the wheel-shoe contact in tread braked wheels, as the cosimulation between the MB simulation and a thermal code should be implemented.