Rail damage prediction is a complex task because it depends on numerous tribological parameters and the dynamic conditions produced by the vehicles operating at different speeds and configurations. Shakedown maps and Whole-Life-Rail-Model/T-Gamma have been used to predict rail damage, but they involve assumptions that may reduce their accuracy. This paper proposes a simulation modelling method to predict rail surface damage based on a locomotive digital twin, calibrated shakedown maps and friction measurements. The method improves the accuracy of rail damage predictions by including slip-dependent friction characteristics, co-simulation of locomotive traction mechatronic system and the mechanical properties of the wheel and rail materials measured through tensile tests. A set of operating conditions are simulated on a high-performance computing cluster, with stress results being post processed into calibrated shakedown heatmaps. The method clearly indicated the influences of the different operating conditions on rail damage for specific combinations of wheel-rail materials and vehicle-track configurations.