Wear and subsurface stress evolution in tractive rolling contact

Contents lists available at ScienceDirect International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci Wear and subsurface stress evolution in tractive rolling contact Javier M. Juliá , Luis Rodríguez-Tembleque ∗ Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, Sevilla 41092, Spain A R T I C L E I N F O Keywords: Contact mechanics Tribology Wear Rolling contact Subsurface stress Orthotropic friction A B S T R A C T Wear phenomenon is inherent to tractive rolling contact problems e.g., in rolling bearings or in rail–wheel interaction. It takes place in the sliding regions of the rolling contact area, and it accumulates as particles from the rolling bodies cross the contact region. Wear modifies the solids’ surface, the contact tractions, the subsurface stresses, and the tangential forces transmitted between the rolling solids, which are fundamental contact variables in sectors such as rolling fatigue or vehicle system dynamics. Thus, ignoring wear in tractive rolling contact analysis could lead to underestimations of bearing fatigue lives or inaccurate lateral guiding forces in multibody vehicle models. This work presents a robust SAM-based formulation on rolling contact to study how wear, contact tractions, resultant rolling contact-forces, and surface and subsurface stresses evolve with the number of revolutions. For the first time, subsurface stress distributions are computed as a function of revolutions under orthotropic friction and wear conditions, highlighting the influence of tribological axes orientation and wear evolution on stress and force reactions. After validating the proposed formulation, several numerical examples are presented to show how considering orthotropic friction and wear laws impacts stress distributions and resultant rolling contact forces. These findings could provide important insights into the role of wear in tractive rolling contact for applications in engineering and industrial design. 1. Introduction Tractive rolling contact is defined as the case of rolling motion where, due to friction, an overall tangential force is transmitted between solids through the contact interface. In tractive rolling, the contact region is divided into stick and slip zones, which are determined as a function of the frictional contact law and the elastic deformation of the solids [1]. Tractive rolling contact has a tremendous interest for many industries and engineering disciplines, e.g., in rolling bearing design for machinery industry – or energy generation industry – [2–6], in rail–wheel contact or in vehicle system dynamics for transportation industry [7–9]. Due to its tremendous interest for the industry development, the research on tractive rolling contact has been developed for more than 100 years – since the pioneer works of Reynolds [10], Hertz [11] and Carter [12] –. In 1876, Reynolds [10] described and measured the creepage phenomenon between a rubber cylinder and a metal plate, proving that the contact region of a rolling contact is divided into stick and slip zones. Later, in 1882, Hertz [11] solved analytically the contact problem between two – smooth surface – elastic solids. However, the treatment of the tractive rolling contact starts with Carter [12] in 1926. Assuming plane strain conditions, he solved the 2D steady-state rolling contact of an elastic cylinder which transmits a tractive force ∗Corresponding author. E-mail address: [email protected] (L. Rodríguez-Tembleque). to the plane on which it is rolling. It allows him to obtain the relation between creepage and creepage forces, applied on locomotive wheels, where tangential forces are transmitted from the wheel to the rail. Later, this two-dimensional theory was extended to the 3D case by Johnson [13,14], considering the longitudinal and lateral creepages in two rolling balls, and by Vermeulen and Johnson [15], extending this solution to arbitrary smooth half-space solids. Bentall and Johnson [16] – and later Nowell and Hills [17] – studied numerically the 2D tractive rolling contact problem between elastically dissimilar rollers. The pioneer numerical methodologies which were able to deal with the 3D elastic tractive rolling contact problem – imposing creepage and spin – were proposed by Kalker [8,18–22]. More recently – and inspired in Kalker’s works –, important contributions have been presented in the tractive rolling contact context. Using the boundary elements method (BEM), González and Abascal [23, 24] developed several techniques for solving 2D steady-state and transient rolling contact problems, and Rodríguez-Tembleque and Abascal [25,26] proposed a methodology to 3D steady-state tractive rollingcontact, which was valid for structured and un-structured meshes in the contact region. In addition, Santeramo et al. developed a 2D BEM methodology to study steady-state viscoelastic circular contact problems [27]. Using Semi-Analytical Models (SAM), Guler et al. [28– https://doi.org/10.1016/j.ijmecsci.2025.110195 Received 30 September 2024; Received in revised form 4 March 2025; Accepted 25 March 2025 International Journal of Mechanical Sciences 294 (2025) 110195 Available online 12 April 2025 0020-7403/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC license ( http://creativecommons.org/licenses/bync/4.0/ ). J.M. Juliá and L. Rodríguez-Tembleque 31] solved 2D graded coatings and functionally graded orthotropic solids, and Xi et al. [32] analyzed 2D steady-state rolling contacts with dry friction. Three-dimensional SAM-based solutions for steadystate tractive rolling contacts were presented by Wang et al. [33] and Xi et al. [34]. Later, they were extended to multilayered bodies by Meshcheryakova and Goryacheva [35] and Xi et al. [36,37], and to viscoplastic multi-layered half-space by Manyo [38] and Wallace et al. [39,40]. SAM – in some cases combined with the Finite Element Method (FEM) – have also allowed to analyze thermo-elastic rolling contact problems [41], conformal rolling contact in wheel– rail or in spherical roller bearing applications [42,43], and non-steady state wheel–rail contact problems [44–47]. Some of these solutions were validated with the Finite Element Method (FEM). Using the FEM, Brhane and Mekonone studied the Von Mises stress evolution in a 2D twin-discs problem [48] and Sharma and Sadeghi [49] developed a 3D model to investigate the effects of fretting-wear on rolling contact fatigue. Furthermore, Zhao and Li [50] presented a 3D transient finite element model for wheel–rail rolling contact, which was extended in [51] to frictional wheel–rail rolling contact in elasto-plasticity, and Yang et al. [52] proposed an explicit finite element methodology for rolling contact problems. All those works allow engineers to have models to predict rolling element fatigue lives (due to the cyclic loading conditions, rolling contact fatigue damage can occur in rolling bearings or wheel–rail contact [2,53–56]), or to estimate the rolling contact-forces in multibody vehicle dynamic models [57–59]. However, due to the presence of slip zones in the contact area, wear phenomena should be also considered in the tractive rolling contact models, since it modifies rolling surfaces and, consequently, the surface (and subsurface) stresses, the slip zones and the resultant rolling contact-forces. For this reason, in the last two decades, a significant number of investigators have attempted to develop theoretical and semi-analytic models to predict wear in rolling contact problems. Olofsson et al. [60] proposed a model for predicting wear on lubricated spherical roller thrust bearings. Jendel et al. [61] and Enblom and Berg [62] developed models to predict wear on train wheel profiles. Using a semi-winkler based model, Telliskivi [63] and Telliskivi and Olofsson [64] simulated wear in rolling–sliding contact on disc-on-disc and in wheel–rail contact. Hegadekatte et al. [65] presented a simplified theoretical model to estimate the maximum wear depth evolution in twin-disc tribometers, and Rodríguez-Tembleque et al. [66,67] presented BEM-based formulations to compute wear under rolling contact conditions Recently, Liu et al. [68] presented a new numerical calculation of railway wheel wear based on the Universal Kriging technique, considering Archard’s wear law and the uncertainties of the wear coefficient. In this context, this work presents a robust SAM-based formulation to study how surface wear, contact tractions, resultant rolling contactforces, and subsurface stresses evolve (with the number of revolutions) in tractive rolling contact problems. We first analyze the tractive rolling contact between two identical spheres. Since the proposed formulation considers orthotropic friction and wear laws, this problem is analyzed under both, isotropic and orthotropic friction laws. Later, we study the tractive rolling contact in twin-discs system under orthotropic friction and wear laws. It allows us to compute the evolution of the worn disc profiles, the tangential contact force components, and the contact pressures or the subsurface von Mises stress distribution in the solids (value and location of the maximum subsurface stresses). Furthermore, we present a novel semi-analytical modeling approach that integrates wear evolution and subsurface stress redistribution in tractive rolling contact. Unlike previous studies that treat wear and stress analysis separately, the presented formulation explicitly associates wear progression with von Mises stress distributions over multiple rolling cycles. In addition, we introduce an innovative discretization scheme for the calculation of wear depth and incorporate orthotropic friction and wear laws within a unified rolling contact framework. These contributions allow us to model a more accurate representation of wear-induced changes in contact mechanics and lay the bases for further research. In particular, this methodology can be extended to include fatigue-life estimation [69–72], thermal effects [73,74], or anisotropic inclusions [75–77], which have been recently studied and are of significant importance in rolling contact analysis. Thus, to the best authors knowledge, this paper presents, for the first time, the calculation of subsurface stress distributions as a function of the number of revolutions, under orthotropic tractive rolling contact conditions. In engineering applications, there are many elements which have orthotropic surface properties due to their manufacturing process (e.g. machining, lamination...) or their nature (i.e. composite materials). The numerical studies presented in this work can highlight that, in tractive rolling contact problems, adequate friction and wear laws should be considered in order to compute precisely both surface and subsurface stresses, as well as the resultant contact forces, which are essential for fatigue-life calculation methods or for multi-body vehicle modeling, respectively. This manuscript is organized as follows. First, the kinematic equations for tractive rolling contact problems are presented in Section 2. Then, in Section 3, the elastic influence coefficient for the half-space approximation is presented. Section 4 presents the orthotropic frictional rolling contact problem description, and Section 5 presents the orthotropic wear law for rolling contact problems. The expressions to compute the subsurface stresses are showed in Section 6. The complete discrete formulation of the whole problem is described in detail in Section 7. It should be noted the novelties on the discrete wear depth calculation proposed in this work. Then, the proposed methodology is applied in Section 8 to study several tractive rolling contact problems under different orthotropic friction and wear conditions. These examples allow us to study how the surface contact tractions, the tangential contact force components or the subsurface von Mises stress distribution in the solids – their maximum values and their locations – evolve with the number of revolutions. Finally, some concluding remarks are presented in Section 9. 2. Rolling contact kinematics The kinematic description of two elastic bodies (𝛺𝛼, 𝛼=𝐴, 𝐵) under rolling contact conditions needs to define a coordinate system 𝑂𝑥𝑦𝑧 which moves with the contact area 𝛤𝑐. The 𝑂𝑥𝑦𝑧 system is located at the center of 𝛤𝑐. The 𝑧-axis is orthogonal to the contact plane and points to solid 𝛺𝐴 (see Fig. 1). The 𝑥-axis points to the direction of rolling and the 𝑦-axis completes the right-handed system. The normal gap between these two solids’ surfaces – in the contact zone position 𝐱 – at a certain time instance 𝜏 can be defined as 𝑔𝑛(𝐱, 𝜏) = 𝑔𝑛,𝑜(𝜏) + 𝑤(𝐱, 𝜏) + 𝑢𝑛(𝐱, 𝜏),(1) where 𝑔𝑛,𝑜(𝜏) = 𝑔𝑔−𝑔𝑜(𝜏) – being 𝑔𝑔 the geometric gap and 𝑔𝑜 the rigid body approachment –, 𝑤 is the wear depth and 𝑢𝑛 is the relative normal displacement: 𝑢𝑛(𝐱, 𝜏) = 𝑢(𝐴) 𝑧(𝐱, 𝜏) − 𝑢(𝐵) 𝑧(𝐱, 𝜏). The slip velocity of solid 𝛺𝐴 with respect to solid 𝛺𝐵 ( 𝐬(𝐱, 𝜏) = [𝑠𝑥(𝐱, 𝜏)𝑠𝑦(𝐱, 𝜏)]𝑇 ) can be written, according to [8,78], as 𝐬(𝐱, 𝜏) = 𝐜(𝐱, 𝜏) − 𝜕𝐮𝑡(𝐱, 𝜏) 𝜕𝐱𝐯(𝐱, 𝜏) + 𝜕𝐮𝑡(𝐱, 𝜏) 𝜕𝜏 ,(2) where 𝐜= 𝐱(𝐴)− 𝐱(𝐵) is the rigid body tangential slip velocity, 𝐯= −(  𝐱(𝐴)+ 𝐱(𝐵))∕2 is the rolling velocity and 𝐮𝑡=𝐮(𝐴) 𝑡−𝐮(𝐵) 𝑡 is the tangential displacement difference. Rolling contact occurs when ‖𝐜‖≪‖𝐯‖=𝑉, being 𝑉 the magnitude of the rolling velocity. Similarly to [8,78], we have assumed that 𝑥-axis points in the rolling direction and, therefore, 𝐯= [𝑉0]𝑇. Moreover, the rigid body tangential slip velocity is expressed as 𝐜= [𝑉(𝜉𝑥−𝜙𝑦)𝑉(𝜉𝑦+𝜙𝑥)]𝑇, where 𝜉𝑥, 𝜉𝑦 and 𝜙 denote the longitudinal, lateral and spin creepages, respectively. Consequently, Eq. (2) can be written as 𝑠𝑥(𝐱, 𝜏) = 𝑉(𝜉𝑥−𝜙𝑦 −𝜕𝑢𝑥(𝐱, 𝜏) 𝜕𝑥 +1 𝑉 𝜕𝑢𝑥(𝐱, 𝜏) 𝜕𝜏 ), 𝑠𝑦(𝐱, 𝜏) = 𝑉(𝜉𝑦+𝜙𝑥 −𝜕𝐮𝑦(𝐱, 𝜏) 𝜕𝑥 +1 𝑉 𝜕𝑢𝑥(𝐱, 𝜏) 𝜕𝜏 ), (3) International Journal of Mechanical Sciences 294 (2025) 110195 2 J.M. Juliá and L. Rodríguez-Tembleque Fig. 1. A elastic body (i.e. a sphere 𝛺𝐴), subjected to a normal load 𝑃 and a tangential load 𝑄, rolls over an elastic halfspace (i.e. 𝛺𝐵), being 𝛤𝐶 the resulting contact area. and for steady-state rolling contact problems, where 𝜕𝐮𝑡∕𝜕𝜏 = 0, as 𝑠𝑥(𝐱) = 𝑉(𝜉𝑥−𝜙𝑦 −𝜕𝑢𝑥(𝐱) 𝜕𝑥 ), 𝑠𝑦(𝐱) = 𝑉(𝜉𝑦+𝜙𝑥 −𝜕𝐮𝑦(𝐱) 𝜕𝑥 ). (4) 3. Half-space approximation To compute the relative displacements 𝐮(𝐱, 𝜏) we use the halfspace approximation proposed by Kalker [22] and employed later by several authors [79–81]. This approach includes several assumptions, among which is the quasi-static description of the problem, which allows the effects of inertia to be neglected. Therefore, under those assumptions explicitly described in previous work [82], the relation between the deformation 𝐮(𝐱, 𝜏) at point 𝐱 and the contact traction 𝐩(𝐱, 𝜏) = [𝐩𝑡(𝐱, 𝜏)𝑇𝑝𝑛(𝐱, 𝜏)]𝑇 on point 𝐱′∈𝛤𝑐 can be expressed as 𝐮(𝐱, 𝜏) = ∫ ∫𝛤𝑐(𝜏) 𝐀(𝐱,𝐱′)𝐩(𝐱′, 𝜏)𝑑𝑥′𝑑𝑦′.(5) In Eq. (5), 𝜏 is a pseudo-time instance and 𝐀(𝐱,𝐱′) is the kernel function matrix which quantifies on of the surface displacement components at 𝐱 induced by one of the traction components of unit magnitude acting at 𝐱′. The expression of the different components of the mentioned kernel function matrix is presented on Appendix A. 4. Tractive rolling contact restrictions 4.1. Normal contact restrictions From the normal gap 𝑔𝑛(𝐱, 𝜏) and the normal contact pressure 𝑝𝑛= 𝑝𝑧(𝐱, 𝜏) we can express the Signorini’s unilateral contact conditions for every point in the contact zone (i.e., 𝐱∈𝛤𝑐) as 𝑔𝑛≥0, 𝑝𝑛≥0, 𝑔𝑛⋅𝑝𝑛= 0 (6) Similarly to Ref. [83], it is possible to summarize the Signorini’s unilateral contact conditions as 𝑝𝑛−PR+(𝑝∗ 𝑛) = 0,(7) where PR+(∙) = 𝑚𝑎𝑥(0,∙) is the normal projection function and 𝑝∗ 𝑛 is the augmented normal traction, which is related to the normal contact pressure and the normal gap through a normal penalty parameter (𝑟𝑛∈ R+). 𝑝∗ 𝑛=𝑝𝑛+𝑟𝑛𝑔𝑛.(8) Fig. 2. Sliding rule in orthotropic tribological axes 𝑒1 and 𝑒2. The 𝛽 angle (counterclock-wise) represents the orientation of the tribological axes [𝑒1, 𝑒2] with respect to the half-space coordinate axes [𝑥, 𝑦]. The 𝜃 angle indicates the sliding direction relative to the tribological axes [𝑒1, 𝑒2] and it is different for each point in the contact zone. 4.2. Tangential contact restrictions Analogously to the normal contact constraints, the fulfillment of an orthotropic friction law is guaranteed by the complementary relations: 𝑓(𝑝𝑛,𝐩𝑡)≤0, 𝜆 ≥0, 𝜆 𝑓(𝑝𝑛,𝐩𝑡) = 0,(9) where, according to [67,82,84–88], 𝐩𝑡 is the tangential contact traction vector (𝐩𝑡= [𝑝𝑥(𝐱, 𝜏)𝑝𝑦(𝐱, 𝜏)]𝑇), 𝑓(𝐩𝑡, 𝑝𝑛) is the orthotropic friction limit function (𝑓(𝐩𝑡, 𝑝𝑛) = ‖𝐩𝑡‖𝜇−𝑝𝑛) and 𝜆=‖𝐬‖∗ 𝜇 [85]. In the previous expressions, the elliptic norms ‖∙‖𝜇 and ‖∙‖∗ 𝜇 are respectively defined so that ‖𝐩𝑡‖𝜇=√(𝑝𝑒1∕𝜇1)2 +(𝑝𝑒2∕𝜇2)2 , ‖𝐬‖∗ 𝜇=√(𝜇1𝑠𝑒1)2+ (𝜇2𝑠𝑒2)2. (10) In the expressions above, 𝜇1 and 𝜇2 are the principal friction coefficients in the tribological main axes {𝑒1, 𝑒2}. The tangential contact tractions and the slip velocity components can respectively be defined in the tribological axes as [𝑝𝑒1 𝑝𝑒2]=[cos 𝛽sin 𝛽 − sin 𝛽cos 𝛽][𝑝𝑥 𝑝𝑦](11) and [𝑠𝑒1 𝑠𝑒2]=[cos 𝛽sin 𝛽 − sin 𝛽cos 𝛽][𝑠𝑥 𝑠𝑦],(12) where 𝛽 is the counter-clock-wise angle between the tribological axes (𝑒1, 𝑒2) and the contact coordinate system axes (𝑥, 𝑦) [82] (see Fig. 2). Therefore, the angle 𝛽 allows orienting the directions where the tribological properties (i.e. 𝜇1, 𝜇2, 𝑖1, 𝑖2) are minimum or maximum due to the orthotropic contact characteristics. Thus, any value of 𝛽∈ (0◦,90◦) implies that the equivalent tribological properties (i.e. equivalent friction intensity coefficient or equivalent intensity wear coefficient) are a combination of the properties existing on the main tribological axes. The frictional contact constraints presented in Eq. (9) can also be expressed as 𝐩𝑡−PE𝜌(𝐩∗ 𝑡) = 0,(13) where 𝐩∗ 𝑡=𝐩𝑡−𝑟𝑡M2𝐬 (being M=𝑑𝑖𝑎𝑔(𝜇1, 𝜇2) and 𝑟𝑡∈R+) is the augmented tangential traction and PE𝜌(∙) ∶ R2⟶R2 is the tangential projection function defined in [85]: PE𝜌(𝐩∗ 𝑡) = {𝐩∗ 𝑡if ‖𝐩∗ 𝑡‖𝜇< 𝜌, 𝜌𝐩∗ 𝑡∕‖𝐩∗ 𝑡‖𝜇if ‖𝐩∗ 𝑡‖𝜇≥𝜌. (14) International Journal of Mechanical Sciences 294 (2025) 110195 3 J.M. Juliá and L. Rodríguez-Tembleque with 𝜌=|PR+(𝑝∗ 𝑛)|. Therefore, Eqs. (7) and (13) constitute the orthotropic frictional contact restrictions for every point: 𝐱∈𝛤𝑐, whose possible states are: No contact: 𝑝𝑛= 0, 𝑔𝑛≥0,𝐩𝑡=𝟎, Contact-Stick: 𝑝𝑛≥0, 𝑔𝑛= 0,𝐬𝑡=𝟎, Contact-Slip: 𝑝𝑛≥0, 𝑔𝑛= 0,𝐩𝑡= −𝑝𝑛M2𝐬∕‖𝐬‖∗ 𝜇. (15) 5. Wear in tractive rolling contact The Holm-Archard’s wear law [89] is widely used to estimate wear in contact problems, and it can be expressed for an infinitesimally small apparent contact area in terms of the wear depth rate (Refs. [66,67, 90–96]), as:  𝑤=𝑖𝑤|𝑝𝑛|‖𝐬‖, being 𝑖𝑤 the wear coefficient. However, in those cases where orthotropic tribological properties are considered in this work, an orthotropic wear law [67,87,88] should be assumed. Therefore, the wear law should be rewritten as  𝑤=|𝑝𝑛|‖𝐬‖𝑖,(16) where ‖𝐬‖𝑖=√(𝑖1𝑠𝑒1)2+ (𝑖2𝑠𝑒2)2, being 𝑖1 and 𝑖2 the principal intensity wear coefficients. In rolling contact, the wear depth rate of a solid particle – for every rotation – could be expressed as  𝑤(𝐱, 𝜏) = 𝜕𝑤(𝐱, 𝜏) 𝜕𝜏 −𝑉𝜕𝑤(𝐱, 𝜏) 𝜕𝑥 ,(17) where it has been assumed that 𝑥-axis is aligned to the rolling direction and, therefore, 𝐯= [𝑉0]𝑇. Consequently, the Holm-Archard’s wear law (16) can be rewritten for rolling contact problems as 𝜕𝑤(𝐱, 𝜏) 𝜕𝜏 −𝑉𝜕𝑤(𝐱, 𝜏) 𝜕𝑥 =|𝑝𝑛|‖𝐬‖𝑖,(18) and for steady-state conditions, where 𝜕𝑤∕𝜕𝜏 = 0, as 𝜕𝑤(𝐱) 𝜕𝑥 = − |𝑝𝑛|‖𝐬‖𝑖 𝑉.(19) Eq. (19) allows us to calculate the distribution of wear depth in the contact zone during each rotation. 6. Subsurface stresses Finally, the subsurface stresses at point 𝐱∈𝛺(𝛼) (𝛼=𝐴, 𝐵) –caused by a surface contact pressure at 𝐱′∈𝛤𝑐– can be computed as 𝝈(𝐱, 𝜏) = ∫ ∫𝛤𝑐(𝜏) 𝐓(𝐱,𝐱′)𝐩(𝐱′, 𝜏)𝑑𝑥′𝑑𝑦′,(20) where the influence coefficients of the kernel function 𝐓(𝐱,𝐱′) depend on the relative positions of the two surface points 𝐱 and 𝐱′, therefore 𝐓(𝐱,𝐱′) = 𝐓(𝐱′−𝐱). For the sake of clarity, the expression above is rewritten as 𝜎𝑖𝑗 (𝐱, 𝜏) = ∫ ∫𝛤𝑐(𝜏) 𝑇𝑆𝑥 𝑖𝑗 (𝐱′−𝐱)𝑝𝑥(𝐱′, 𝜏) + 𝑇𝑆𝑦 𝑖𝑗 (𝐱′−𝐱)𝑝𝑦(𝐱′, 𝜏) +𝑇𝑁 𝑖𝑗 (𝐱′−𝐱)𝑝𝑛(𝐱′, 𝜏)𝑑𝑥′𝑑𝑦′.(21) The explicit expressions of Kernel functions 𝑇𝑁 𝑖𝑗 (𝑥, 𝑦, 𝑧), 𝑇𝑆𝑥 𝑖𝑗 (𝑥, 𝑦, 𝑧) and 𝑇𝑆𝑦 𝑖𝑗 (𝑥, 𝑦, 𝑧) have been already presented in Ref. [82]. 7. Discrete rolling contact formulation For steady-state rolling contact problems, the surface wear and the subsurface stress evolution can be computed by solving the nonlinear equations set (1), (4), (5), (7), (13), (19), (21). For that purpose that equations set should be discretized and solved for every solids rotation (𝑘). Fig. 3. Computational domain covering the potential contact zone with a 𝑁𝑒=𝑁𝑥×𝑁𝑦 rectangular mesh formed by elements of size 𝛥𝑥×𝛥𝑦. 7.1. Discrete contact problem approximation For the numerical solution of the contact problem, a rectangular potential contact zone is considered and discretized by a regular mesh with 𝑁𝑒=𝑁𝑥×𝑁𝑦 elements of size 𝛥𝑥×𝛥𝑦 (see Fig. 3, where the coordinates of the center of the element 𝐼 are denoted by 𝐱𝐼= [𝑥𝐼, 𝑦𝐼,0]𝑇. Consequently, tractions 𝐩 in the surface integral Eq. (5) are approximated by element-wise constant functions and, therefore, the integral Eq. (5) can be rewritten, for each contacting element I, as the summation of the products of influence coefficients 𝐴𝑖𝑗 𝐼𝐽 (𝑖, 𝑗 corresponding to 𝑥, 𝑦 and 𝑧) and nodal pressures 𝑝𝐽 𝑗 at element 𝐽: 𝑢(𝑘) 𝐼𝑥 = 𝑁𝑒 ∑ 𝐽=1 𝐴𝑥𝑥 𝐼𝐽 𝑝(𝑘) 𝐽𝑥 +𝐴𝑥𝑦 𝐼𝐽 𝑝(𝑘) 𝐽𝑦 +𝐴𝑥𝑧 𝐼𝐽 𝑝(𝑘) 𝐽 𝑧 , 𝑢(𝑘) 𝐼𝑦 = 𝑁𝑒 ∑ 𝐽=1 𝐴𝑦𝑥 𝐼𝐽 𝑝(𝑘) 𝐽𝑥 +𝐴𝑦𝑦 𝐼𝐽 𝑝(𝑘) 𝐽𝑦 +𝐴𝑦𝑧 𝐼𝐽 𝑝(𝑘) 𝐽 𝑧 , 𝑢(𝑘) 𝐼𝑧 = 𝑁𝑒 ∑ 𝐽=1 𝐴𝑧𝑥 𝐼𝐽 𝑝(𝑘) 𝐽𝑥 +𝐴𝑧𝑦 𝐼𝐽 𝑝(𝑘) 𝐽𝑦 +𝐴𝑧𝑧 𝐼𝐽 𝑝(𝑘) 𝐽 𝑧 . (22) In the expressions above, 𝑢(𝑘) 𝐼𝑖 is the deformation in 𝑖-direction of element 𝐼 for every rotation (𝑘). The influence coefficient 𝐴𝑖𝑗 𝐼𝐽 represents the influence on the 𝑖-direction deformation of element 𝐼, caused by a unit 𝑗-direction traction on another element 𝐽. It is computed by integrating Eq. (5) over a single element 𝐽, with respect to an observation point at the center of element 𝐼 (see [22,97,98]). The explicit expressions of the influence coefficients 𝐴𝑖𝑗 𝐼𝐽 were presented by the authors in [82]. Assembling the coefficients 𝐴𝑖𝑗 𝐼𝐽 in the matrix 𝐀𝑖𝑗 : (𝐀𝑖𝑗 )𝐼𝐽 =𝐴𝑖𝑗 𝐼𝐽 , the Eq. (22) can be expressed as 𝐮(𝑘)=𝐀 𝐩(𝑘), i.e., ⎡⎢⎢⎣ 𝐮𝑥 𝐮𝑦 𝐮𝑧⎤⎥⎥⎦ (𝑘) =⎡⎢⎢⎣ 𝐀𝑥𝑥 𝐀𝑥𝑦 𝐀𝑥𝑧 𝐀𝑦𝑥 𝐀𝑦𝑦 𝐀𝑦𝑧 𝐀𝑧𝑥 𝐀𝑧𝑦 𝐀𝑧𝑧⎤⎥⎥⎦⎡⎢⎢⎣ 𝐩𝑥 𝐩𝑦 𝐩𝑧⎤⎥⎥⎦ (𝑘) ,(23) where vector 𝐮(𝑘) and vector 𝐩(𝑘) contain the deformation components and the contact tractions of each element 𝐼 for every rotation (𝑘), respectively. Similarly to [82], Eq. (23) is conveniently rearranged as [𝐮𝑛 𝐮𝑡](𝑘) =[𝐀𝑛𝑛 𝐀𝑛𝑡 𝐀𝑡𝑛 𝐀𝑡𝑡 ][𝐩𝑛 𝐩𝑡](𝑘) ,(24) where 𝐮(𝑘) 𝑛=𝐮(𝑘) 𝑧, 𝐩(𝑘) 𝑛=𝐩(𝑘) 𝑧, 𝐮(𝑘) 𝑡 contains the tangential deformation components of each element 𝐼 (i.e., (𝐮(𝑘) 𝑡)𝐼= [ (𝐮(𝑘) 𝑥)𝐼(𝐮(𝑘) 𝑦)𝐼]𝑇) and 𝐩(𝑘) 𝑡 the contact traction components of each element 𝐼 (i.e., (𝐩(𝑘) 𝑡)𝐼= [ (𝐩(𝑘) 𝑥)𝐼(𝐩(𝑘) 𝑦)𝐼]𝑇). Eq. (24) constitutes the discrete representation of Eq. (5). International Journal of Mechanical Sciences 294 (2025) 110195 4 J.M. Juliá and L. Rodríguez-Tembleque Fig. 4. Notation of the mesh elements coincident with the rotation plane (i.e. 1,2...𝐼 − 1, 𝐼, 𝐼 +1...𝑁𝑥) where the potential contact zone, the dimensions of the computational domain 𝑁𝑒=𝑁𝑥×𝑁𝑦 and the size of the elements 𝛥𝑥×𝛥𝑦 are shown. 7.2. Discrete contact kinematic variables The normal gap (1) on each element 𝐼, at rotation 𝑘, can be expressed as (𝐠(𝑘) 𝑛)𝐼= (𝐠(𝑘) 𝑧,𝑜)𝐼+ (𝐰(𝑘))𝐼+ (𝐮(𝑘) 𝑛)𝐼,(25) where 𝐠(𝑘) 𝑧,𝑜 is the bodies’ rapprochement vector, and vectors 𝐰(𝑘) and 𝐮(𝑘) 𝑛 contain the wear depth and the normal deformation of each element, respectively. The tangential slip velocity (4) on each element 𝐼, at rotation 𝑘, can be written as (𝐬(𝑘) 𝑡)𝐼= (𝐜(𝑘))𝐼+𝐷𝑟(𝐮(𝑘) 𝑡)𝐼 ,(26) where (𝐜(𝑘))𝐼=𝑉[𝜉𝑥−𝜙 𝑦𝐼 𝜉𝑦+𝜙 𝑥𝐼],(27) 𝐷𝑟(𝐮(𝑘) 𝑡)𝐼=−𝑉 𝛥𝑥[(𝐮(𝑘) 𝑥)𝐼+1 − (𝐮(𝑘) 𝑥)𝐼 (𝐮(𝑘) 𝑦)𝐼+1 − (𝐮(𝑘) 𝑦)𝐼],(28) and vector 𝐮(𝑘) 𝑡 contains the tangential deformation of each element. In the Eq. (28), 𝐷𝑟((𝐮(𝑘) 𝑡)𝐼) is the first-order upwind derivative approximation of the tangential deformation of each element 𝐼, i.e., the element 𝐼+1 being the one whose coordinates, in this case, are: 𝑥𝐼+1 = 𝑥𝐼+𝛥𝑥 and 𝑦𝐼+1 =𝑦𝐼 (see Fig. 4). 7.3. Discrete wear depth An innovative aspect of the proposed formulation is how the wear depth on each element 𝐼 is computed on every rotation (𝑘) from the discrete expression of Eq. (19), i.e., from 𝐷𝑟(𝐰(𝑘))𝐼=|(𝐩(𝑘) 𝑛)𝐼|‖(𝐬(𝑘) 𝑡)𝐼‖𝑖.(29) In Eq. (29), vector 𝐰(𝑘) contains the wear depth on each element 𝐼 and 𝐷𝑟(𝐰(𝑘) 𝑡)𝐼=−𝑉 𝛥𝑥((𝐰(𝑘))𝐼+1 − (𝐰(𝑘))𝐼)(30) is the first-order upwind derivative approximation of the wear depth of each element 𝐼, i.e., the element 𝐼+ 1 being the one whose coordinates are: 𝑥𝐼+1 =𝑥𝐼+𝛥𝑥 and 𝑦𝐼+1 =𝑦𝐼. The system of equations given in Eqs. (29) and (30) can be solved independently for each set of elements that share the same 𝑦-coordinate (i.e., 𝑦𝐼+1 =𝑦𝐼) (𝐰(𝑘))𝐼= (𝐰(𝑘))𝐼+1 + (𝐩(𝑘) 𝑛)𝐼‖(𝐬(𝑘) 𝑡)𝐼‖𝑖 𝛥𝑥 𝑉.(31) Since 𝑥𝐼+1 =𝑥𝐼+𝛥𝑥, the Eq. (31) can be repeated recursively, then we obtain the following equation (𝐰(𝑘))𝐼= (𝐰(𝑘))𝑁𝑥+ 𝑁𝑥−1 ∑ 𝐽=𝐼 (𝐩(𝑘) 𝑛)𝐽‖(𝐬(𝑘) 𝑡)𝐽‖𝑖 𝛥𝑥 𝑉,(32) where, according to Fig. 4, 𝐼 index refers to the elements 1,…,(𝑁𝑥−1) and the element 𝑁𝑥 refers to the element where the solid particles come to the potential contact zone. The assumption that we are studying a complete steady state rolling process allows us to know the wear depth on element 𝑁𝑥 from the previous rotation (𝑘− 1), i.e., (𝐰(𝑘))𝑁𝑥= (𝐰(𝑘−1))1. Therefore, Eq. (32) can be rewritten as (𝐰(𝑘))𝐼= (𝐰(𝑘−1))𝑁𝑥+ 𝑁𝑥−1 ∑ 𝐽=𝐼 (𝐩(𝑘) 𝑛)𝐽‖(𝐬(𝑘) 𝑡)𝐽‖𝑖 𝛥𝑥 𝑉.(33) Finally, the wear depth on each solid 𝛺𝛼 (𝛼=𝐴, 𝐵) surface is computed from the total wear depth ((𝐰(𝑘))𝐼), similarly to previous works [66,91], as (𝐰(𝑘))(𝐴) 𝐼=𝐻(𝐵) 𝐻(𝐴)+𝐻(𝐵)(𝐰(𝑘))𝐼,(𝐰(𝑘))(𝐵) 𝐼=𝐻(𝐴) 𝐻(𝐴)+𝐻(𝐵)(𝐰(𝑘))𝐼,(34) where 𝐻(𝐴) and 𝐻(𝐵) are the surface hardness of each solid (𝐴) and (𝐵), respectively. According to the Holm-Archard’s wear law [89], the wear volume is inversely proportional to the surface hardness of the materials in contact. 7.4. Discrete subsurface stresses equations The subsurface stresses can be computed from discrete expression of Eq. (21) for a set of interior points (i.e., 𝐱𝐼∈𝛺(𝛼) (𝛼=𝐴, 𝐵) being 𝐼= 1...𝑁𝑖), as (𝜎(𝑘) 𝑖𝑗 )𝐼= 𝑁𝑒 ∑ 𝐽=1 (𝐵𝑆𝑥 𝑖𝑗 )𝐼𝐽 𝑝(𝑘) 𝐽𝑥 + (𝐵𝑆𝑦 𝑖𝑗 )𝐼𝐽 𝑝(𝑘) 𝐽𝑦 + (𝐵𝑁 𝑖𝑗 )𝐼𝐽 𝑝(𝑘) 𝐽𝑧 (35) where (𝜎(𝑘) 𝑖𝑗 )𝐼 is the value for the stress at point 𝐱𝐼 on instant 𝜏𝑘: 𝜎𝑖𝑗 (𝐱𝐼, 𝜏𝑘), due to the known surface contact tractions 𝐩(𝑘). Moreover, the influence coefficients (𝐵𝑆𝑥 𝑖𝑗 )𝐼𝐽 , (𝐵𝑆𝑦 𝑖𝑗 )𝐼𝐽 and (𝐵𝑁 𝑖𝑗 )𝐼𝐽 are computed as (𝐵𝑆𝑥 𝑖𝑗 )𝐼𝐽 =∫𝑥+ 𝑥−∫𝑦+ 𝑦− 𝑇𝑆𝑥 𝑖𝑗 (𝑥, 𝑦, 𝑧𝐼)𝑑𝑥′𝑑𝑦′, (𝐵𝑆𝑦 𝑖𝑗 )𝐼𝐽 =∫𝑥+ 𝑥−∫𝑦+ 𝑦− 𝑇𝑆𝑦 𝑖𝑗 (𝑥, 𝑦, 𝑧𝐼)𝑑𝑥′𝑑𝑦′, (𝐵𝑁 𝑖𝑗 )𝐼𝐽 =∫𝑥+ 𝑥−∫𝑦+ 𝑦− 𝑇𝑁 𝑖𝑗 (𝑥, 𝑦, 𝑧𝐼)𝑑𝑥′𝑑𝑦′. (36) The influence coefficient (𝐵𝑁 𝑖𝑗 )𝐼𝐽 represents the influence to the 𝜎𝑖𝑗 - stress on 𝐱𝐼, caused by a unit normal traction applied on element 𝐽. Similarly, coefficient (𝐵𝑆𝑥 𝑖𝑗 )𝐼𝐽 (or (𝐵𝑆𝑦 𝑖𝑗 )𝐼𝐽 ) represents the influence to the 𝜎𝑖𝑗 -stress caused by a unit tangential traction on 𝑥-direction (or on 𝑦-direction), applied on element 𝐽. The subsurface influence coefficients (36) were calculated according to [99], and explicitly presented in [82]. Finally, Eq. (35) can be presented in the following matrix form: 𝝈(𝑘) 𝑖𝑗 =𝐁𝑆𝑥 𝑖𝑗 𝐩(𝑘) 𝑥+𝐁𝑆𝑦 𝑖𝑗 𝐩(𝑘) 𝑦+𝐁𝑁 𝑖𝑗 𝐩(𝑘) 𝑧(37) where 𝝈(𝑘) 𝑖𝑗 is a 𝑁𝑖-dimensional vector containing the value of the stress component 𝜎𝑖𝑗 on every subsurface point 𝐱𝐼∈𝛺(𝛼) for the 𝑘 rotation. 7.5. Solution scheme The discrete equations set (22)–(36) allows us to compute the subsurface stresses evolution caused by orthotropic wear and contact conditions. That set of equations can be solved using an iterative Uzawa scheme The efficiency of this method has already been evaluated in mechanical contact problems in several works [85,100,101], where International Journal of Mechanical Sciences 294 (2025) 110195 5 J.M. Juliá and L. Rodríguez-Tembleque it is also compared with other existing methods (i.e., Newton algorithms), and using the finite element method or the boundary element method. Although its convergence is slow (i.e., it exhibits only linear convergence in general), Uzawa’s method has demonstrated to be very robust and that is why it has also been employed in rolling contact problems [66,67]. To compute the variables on every rotation (𝑘), (I) Initialize 𝐩(0) 𝑛=𝐩(𝑘−1) 𝑛, 𝐩(0) 𝑡=𝐩(𝑘−1) 𝑡, 𝐰(0) =𝐰(𝑘−1) and iterate using (𝑛= 0,1,2,3,…) index. (II) Compute the normal gap vector and the tangential displacements: [𝐠𝑛 𝐮𝑡](𝑛+1) =[𝐠𝑛,𝑜 𝟎](𝑘) +[𝐰 𝟎](𝑛) +[𝐀𝑛𝑛 𝐀𝑛𝑡 𝐀𝑡𝑛 𝐀𝑡𝑡 ][𝐩𝑛 𝐩𝑡](𝑛) , (III) Compute the slip velocity: (𝐬(𝑛+1) 𝑡)𝐼= (𝐜(𝑘))𝐼+𝐷𝑟(𝐮(𝑛+1) 𝑡)𝐼 (IV) Update contact tractions on every element I: (𝐩(𝑛+1) 𝑛)𝐼=PR+( (𝐩(𝑛) 𝑛)𝐼+𝑟𝑛(𝐠(𝑛+1) 𝑛)𝐼) (𝐩(𝑛+1) 𝑡)𝐼=PE𝜌( (𝐩(𝑛) 𝑡)𝐼−𝑟𝑡M2(𝐬(𝑛+1) 𝑡)𝐼), where 𝜌= (𝐩(𝑛+1) 𝑛)𝐼. (V) Update accumulated wear depth: (𝐰(𝑛+1))𝐼= (𝐰(𝑘−1))𝑁𝑥+ 𝐼 ∑ 𝐽=1 (𝐩(𝑛+1) 𝑛)𝐽‖(𝐬(𝑛+1) 𝑡)𝐽‖𝑖 𝛥𝑥 𝑉. (VI) Compute the error: 𝛹(𝑛+1) =‖𝐩(𝑛+1) 𝑛−𝐩(𝑛) 𝑛‖+‖𝐩(𝑛+1) 𝑡−𝐩(𝑛) 𝑡‖. (a) If 𝛹(𝑛+1) ≤𝜀, the solution for the rotation (k) is reached. If the applied boundary condition are the external normal load (𝑃(𝑘)), the resultant contact loads on the contact zone (𝛤𝑐) have to be also computed before reaching the solution for rotation (k): 𝑃(𝑛+1) =𝛥𝑥𝛥𝑦 𝑁𝑒 ∑ 𝐼=1 (𝐩(𝑛+1) 𝑛)𝐼,(38) 𝐐(𝑛+1) =𝛥𝑥𝛥𝑦 𝑁𝑒 ∑ 𝐼=1 (𝐩(𝑛+1) 𝑡)𝐼,(39) therefore, (a.1) If 𝑃(𝑛+1) > 𝑃 (𝑘)+𝜀𝑙𝑜𝑎𝑑 or 𝑃(𝑛+1) < 𝑃 (𝑘)−𝜀𝑙𝑜𝑎𝑑 , modify 𝐠(𝑘) 𝑜,𝑛 →𝐠(𝑘) 𝑜,𝑛 ∓𝛥𝐠(𝑛) 𝑜,𝑛 and return to (II). (a.2) Otherwise, the solution for instant (k) is reached, so 𝐩(𝑘) 𝑛=𝐩(𝑛+1) 𝑛, 𝐩(𝑘) 𝑡=𝐩(𝑛+1) 𝑡, 𝐠(𝑘) 𝑛=𝐠(𝑛+1) 𝑛, 𝐬(𝑘) 𝑡= 𝐬(𝑛+1) 𝑡 and 𝐰(𝑘)=𝐰(𝑛+1). (b) Otherwise, return to (II) evaluating: 𝐩(𝑛) 𝑛=𝐩(𝑛+1) 𝑛, 𝐩(𝑛) 𝑡= 𝐩(𝑛+1) 𝑡 and 𝐰(𝑛)=𝐰(𝑛+1), and iterate until the convergence is reached. (VII) Finally, compute the subsurface stresses for the rotation (𝑘), i.e., 𝝈(𝑘) 𝑖𝑗 , according to the Eq. (37). Once the solution for rotation (k) is reached, the solution for the next rotation is achieved by returning to (I). 8. Numerical analysis In this section, the proposed formulation is applied to analyze two tractive rolling contact problems. The first example is presented in Section 8.1 and solves the rolling contact between two spheres under different isotropic and orthotropic friction conditions. This example allows us to deeply study the influence of friction intensity, and the importance of 𝛽 – i.e., the angle between the tribological axis (𝑒1) and the rolling direction (𝑥-axis) –, in the surface and subsurface stresses, and in the resultant tangential contact forces. In the second example (Section 8.2), we analyze the role of wear in tractive rolling contact problems. For this purpose, a twin-discs rolling system is analyzed under orthotropic friction and wear laws. This example shows the evolution of the solids profiles, the surface and subsurface stresses, and the resultant tangential contact forces with the number of revolutions. Moreover, the influence of the 𝛽 angle is also analyzed. In both sections, the 𝛽 angles of 0◦, 45◦ and 90◦ have been selected as representative for describing the evolution of the main variables in rolling contact problems based on the orientation of the tribological axes. At 𝛽= 0◦ the tribological axes align with the reference axes (i.e. 𝑒1 parallel to 𝑥 and 𝑒2 parallel to 𝑦); at 𝛽= 45◦ the main tribological axes do not coincide with any reference axis, representing a case where the resulting properties are a combination of the properties existing in both directions (i.e., an equivalent property); and at 𝛽= 90◦ the situation is the opposite of 𝛽= 0◦, since the 𝑒1 axis is aligned with the 𝑦-direction, while the 𝑒2 axis is aligned with the 𝑥-axis (see Fig. 2). 8.1. Tractive rolling contact between two identical spheres 8.1.1. Rolling contact under isotropic friction First, we validate the proposed methodology by solving the tractive rolling contact problem presented by Manyo in [38], i.e., two identical spheres of radius 𝑅= 337.5 mm and the same elastic properties (i.e., shear modulus 𝐺(𝐴)=𝐺(𝐵)=𝐺= 1 MPa and the Poisson’s ratio 𝜈(𝐴)=𝜈(𝐵)= 0.28) which are pressed against each other and then subjected to a rolling movement relative to each other. An illustration of the problem is shown in the Fig. 5(a). The applied normal force 𝑃 is worth 0.4705 N and the tangential force 𝑄 (following 𝑥 axis) is 0.657𝜇𝑃 (𝜉𝑥= −0.0031, 𝜉𝑦= 0 and 𝜙= 0), where the friction coefficient is 𝜇= 0.4013 (i.e., 𝜇1=𝜇2=𝜇 in the friction law presented in Section 4.2). Only the steady state case is considered. According to the Hertz solution, the applied normal force (𝑃= 0.4705 N) produces a maximum normal pressure value 𝑝𝑜= 0.01834 MPa and a contact zone whose radius is 𝑎𝑜= 3.5 mm. Fig. 5(b) presents the normalized tangential contact traction 𝑝𝑥∕𝜇𝑝𝑜 distributions on the planes 𝑦= 0 and 𝑦= 0.75𝑎𝑜, showing a good correlation between the computed results and the solution presented by Manyo [38]. Therefore, it validates our model. In our simulations, the computation domain size considered for this problem has been 𝐿𝑥=𝐿𝑦= 4.08 mm, and the domain has been discretized with 91 × 91 mesh elements. These mesh parameters have been selected to ensure the mapping of the contact zone throughout the simulation with a number of elements that ensures the accurate calculation of results in an efficient manner. In this solution, the – normal and tangential – penalty parameters have been fixed to 𝑟𝑛= 8.0⋅10−3 and 𝑟𝑡= 5.0⋅10−7, and the convergence is achieved for 𝜀= 10−11. Fig. 5(c) shows the normalized contact pressure 𝑝𝑛∕𝑝𝑜 distribution on the surface contact zone and Fig. 5(d) presents normalize tangential contact traction 𝑝𝑡∕𝜇𝑝𝑜 distribution. We can thus distinguish the sliding zone at the back of the contact zone and that of adhesion at the front of the contact zone. In addition, Fig. 5(e–f) show the normalized von Mises 𝜎𝑉 𝑀 ∕𝑝𝑜 subsurface stress distributions on the planes 𝑦= 0 and 𝑦= 0.75𝑎𝑜, respectively. In these figures, it is possible to observe the existence of two von Mises stress maximums. One of the maximum values is located inside the solid. The other one is located on the solid surface, at the boundary between the adhesion and the sliding zones. This is in perfect accordance with the subsurface stress distributions observed in [35,38] for tractive rolling contact problems. International Journal of Mechanical Sciences 294 (2025) 110195 6 J.M. Juliá and L. Rodríguez-Tembleque Fig. 5. (a) Scheme of two elastically similar spheres (𝛺𝐴 and 𝛺𝐵) in rolling contact pressed together by a normal load 𝑃. (b) Comparison of normalized tangential contact traction 𝑝𝑥∕𝜇𝑝𝑜 distributions on the planes 𝑦= 0 and 𝑦= 0.75𝑎𝑜 between computed results and the results obtained by Manyo [38]. (c) Normalized contact pressure 𝑝𝑛∕𝑝𝑜 distribution and (d) normalized tangential contact traction 𝑝𝑡∕𝜇𝑝𝑜 distribution on the surface contact zone. (e–f) Normalized von Mises equivalent subsurface stress distributions 𝜎𝑉 𝑀 ∕𝑝𝑜 on the planes 𝑦= 0 and 𝑦= 0.75𝑎𝑜, respectively. The 𝑥 and 𝑦 axis are non-dimensional as they are expressed relative to the initial contact semi-width (𝑎𝑜). International Journal of Mechanical Sciences 294 (2025) 110195 7 J.M. Juliá and L. Rodríguez-Tembleque Fig. 6. (a) Normalized tangential resultant components 𝑄𝑥∕𝜇𝑃 and 𝑄𝑦∕𝜇𝑃 computed as a function of the creep parameter  𝜉, for tractive rolling contact under isotropic friction compared with the resultant components from Johnson [15] and Kalker [19]. (b) Normalized tangential tractive rolling contact resultants 𝑄𝑥∕𝜇1𝑃 and 𝑄𝑦∕𝜇1𝑃 –under orthotropic friction–, as a function of the creep parameter  𝜉, for the following values of the angle between the tribological axes and the coordinate system ones: 𝛽= {0◦,45◦,90◦}. Additionally, Fig. 6(a) presents the normalized tangential resultant components 𝑄𝑥∕𝜇𝑃 and 𝑄𝑦∕𝜇𝑃 as a function of the following creep parameter,  𝜉= −𝜉𝑥 16𝐺𝑎2 𝑜 3𝜇𝑃 (4 − 3𝜈).(40) The numerical results are presented by comparison with the theoretical solution proposed by Vermeulen and Johnson in [15], 𝑄𝑥=𝜇𝑃 (1 − (1 −  𝜉)3), 𝑄𝑦= 0,(41) and with the empirical theoretical solution proposed by Kalker in [19], 𝑄𝑥=𝜇𝑃 (3 2  𝜉arccos(  𝜉)+1−(1 +  𝜉2∕2)√1 −  𝜉2), 𝑄𝑦= 0.(42) In both the cases, during the tractive rolling contact regime, 𝑄𝑥 increases with  𝜉 until it reaches a constant value, i.e., when the sliding rolling contact regime is achieved. Zaazaa and Schwab [102] pointed that, during the tractive rolling contact regime, for all  𝜉 , the Vermeulen and Johnson theory predicts higher values of the tangential resultant than Kalker’s theory. In Fig. 6(a), our numerical results show a very high agreement between our numerical results – square markers – and Kalker’s theoretical solution – continuous red lines –. 8.1.2. Rolling contact under orthotropic friction After validating the accuracy of our numerical scheme for tractive rolling contact problems, the algorithm is extended to include the orthotropic tribological conditions presented in Section 4.2, being 𝜇1≠ 𝜇2, i.e., 𝜇1= 0.4 and 𝜇2= 0.2. In this case, we study the influence of 𝛽, i.e., the angle between the tribological axis (𝑒1,-axis) and the rolling direction (𝑥-axis), considering several values for 𝛽= {0◦,45◦,90◦}. As we have selected 𝜇1≥𝜇2, the value of the effective friction intensity coefficient in the rolling direction will decrease with the 𝛽 angle. The computed normalized tangential tractive rolling contact resultants 𝑄𝑥∕𝜇1𝑃 and 𝑄𝑦∕𝜇1𝑃 are presented in Fig. 6(b) as a function of the creep parameter  𝜉= −𝜉𝑥16𝐺𝑎2 𝑜∕3𝜇1𝑃(4 − 3𝜈), for 𝛽= {0◦,45◦,90◦}. In all the cases, 𝑄𝑥 increases with  𝜉 until it reaches a constant value, i.e., when the sliding rolling contact regime is achieved. It should be noted that, the lower effective friction intensity coefficient values we have, the sooner we reach the sliding rolling contact regime. Moreover, Fig. 6(b) reveals that the value of the 𝑄𝑥 component decreases with the effective friction intensity coefficient in the rolling direction. When the sliding regime is achieved, 𝑄𝑥(𝛽= 90◦) = 𝑄𝑥(𝛽= 0◦)∕2, since 𝜇2=𝜇1∕2 and 𝜇1 is the effective friction intensity coefficient when 𝛽= 0◦ (and 𝜇2, when 𝛽= 90◦). Thus, those are expected results. For 𝑄𝑦 component, expected results are also observed for 𝛽= 0◦ and 𝛽= 90◦, i.e., a zero value for 𝑄𝑦 component is obtained for 𝛽= 0◦ and 𝛽= 90◦. However, a nonzero value for 𝑄𝑦 is observed for 𝛽= 45◦. This is due to the fact that nonzero 𝑝𝑦 contact traction values are obtained for 0◦< 𝛽 < 90◦, when an associative sliding rule is considered for our orthotropic friction law [82]. This is clearly observed in Fig. 7, where the normalized tangential contact traction module ‖𝐩𝑡‖𝜇∕𝑝𝑜 distributions (upper figure), and the normalized tangential contact traction distribution components (𝑝𝑥∕𝑝𝑜 and 𝑝𝑦∕𝑝𝑜) on the planes: 𝑦= 0 and 𝑦= 0.75𝑎𝑜 (lower figure), are presented for 𝛽= {0◦,45◦,90◦}. To conclude with this section, Fig. 8 presents the normalized von Mises stress distributions (𝜎𝑉 𝑀 ∕𝑝𝑜) in xz plane (upper figure) and in yz plane (lower figure), for 𝛽= {0◦,45◦,90◦}. In the xz plane there are stress maximums both inside the solid and at the surface of the solid, but only the interior one is maintained as the equivalent friction intensity coefficient decreases (i.e., the 𝛽 angle increases). Moreover, in both planes – the xz plane and the yz plane –, we observe that the maximum value of the von Mises stress inside the solid remains approximately the same (i.e., 𝜎𝑉 𝑀 ∕𝑝𝑜≈ 0.6) as the 𝛽 angle increases, i.e., as the equivalent friction intensity coefficient decreases. However, at the solid surface, it clearly decreases as 𝛽 increases. The reason behind this evolution of the subsurface von Mises stresses is that the maximum stress located in the interior of the solid is associated with the frictionless contact of two solids itself [35,82], while the maximum stress values appearing at the surface are associated with friction. That is why as the equivalent friction intensity coefficient decreases (i.e. 𝛽 increases) the value of the stresses closer to the surface decrease while those closer to the interior of the solid remain approximately constant. Finally, it should be also noted that, in the yz plane, due to nonzero 𝑝𝑦 contact traction values are obtained for 0◦< 𝛽 < 90◦, the von Mises stress distribution is not symmetric with respect to the xz plane (see von Mises stress distribution in the yz plane for 𝛽= 45◦). 8.2. Tractive rolling contact in twin-discs system After validating the proposed formulation to compute wear in rolling contact problems (see Appendix B), this section studies the influence of wear solving the twin-discs system presented in Fig. 9(a). Thus, this case incorporates wear into tractive rolling contact problems, International Journal of Mechanical Sciences 294 (2025) 110195 8 J.M. Juliá and L. Rodríguez-Tembleque Fig. 7. Normalized tangential contact traction ‖𝐩𝑡‖𝜇∕𝑝𝑜 distributions on the xy plane (upper figures), and the normalized tangential contact traction distribution components (𝑝𝑥∕𝑝𝑜 and 𝑝𝑦∕𝑝𝑜) on the planes: 𝑦= 0 and 𝑦= 0.75𝑎𝑜 (lower figures), for the following values of the angle between the tribological axes and the coordinate system ones: (a) 𝛽= 0◦, (b) 𝛽= 45◦ and (c) 𝛽= 90◦. The 𝑥 and 𝑦 axis are non-dimensional as they are expressed relative to the initial contact semi-width (𝑎𝑜). Fig. 8. Normalized Von Mises stress distributions (𝜎𝑉 𝑀 ∕𝑝𝑜) in xz plane (upper figures) and in yz plane (lower figures), for the following values of the angle between the tribological axes and the coordinate system ones: (a) 𝛽= 0◦, (b) 𝛽= 45◦ and (c) 𝛽= 90◦. The 𝑥 and 𝑦 axis are non-dimensional as they are expressed relative to the initial contact semi-width (𝑎𝑜). allowing us to observe the evolution of contact pressures and subsurface stresses as the contact zone is modified. Additionally, considering orthotropic friction and wear laws allows us to model a more realistic representation of cases where solid surfaces (or materials) exhibit direction-dependent friction and wear behavior. The upper disc 𝛺𝐴 of radii 𝑅𝐴,𝑥 = 32.5 mm and 𝑅𝐴,𝑦 = 32.5 mm is the driving disc and it is subjected to a normal force 𝑃= 300𝑁 and to a rotation velocity 𝜔𝐴= 300 rev/min. The lower disc 𝛺𝐵 is the driven disc and its radii are 𝑅𝐵,𝑥 = 32.3 mm and 𝑅𝐵,𝑦 = ∞. Both discs have the same material properties, i.e., elasticity modulus 𝐸(𝐴)=𝐸(𝐵)= 208×103 MPa and the Poisson’s ratio 𝜈(𝐴)=𝜈(𝐵)= 0.30. In this example, the discs are subjected to tractive rolling contact conditions, being the 𝜉𝑥= −0.005, 𝜉𝑦= 0 and 𝜙= 0. Moreover, orthotropic tribological laws are considered for friction and wear. The following values have been taken for the principal axes friction coefficients: 𝜇1= 0.60 and 𝜇2=𝜇1∕2, and for the intensity wear coefficients: 𝑖1= 2.0⋅10−6 MPa−1 and 𝑖2=𝑖1∕2. Similarly to the example presented in Section 8.1.2, when the value of the 𝛽 angle increases, the equivalent friction intensity coefficient decreases. In the following simulations, the computation domain size parameters (𝐿𝑥 and 𝐿𝑦) have been set, respectively to 𝐿𝑥= 0.35 mm and 𝐿𝑦= 1.4 mm, and the domain has been discretized with 41 × 121 mesh elements. These mesh parameters have been selected to ensure the mapping of the contact zone throughout the simulation with a number of elements that ensures the accurate calculation of results in International Journal of Mechanical Sciences 294 (2025) 110195 9 J.M. Juliá and L. Rodríguez-Tembleque Fig. 15. Comparison of the evolution of the maximum wear depth, in millimeters, as a function of the number of rotations (𝑁) between the computed results and the GIWM method presented by Hegadekatte et al. [65,107]. Table 1 GIWM algorithm for twin-disc rolling system. (I) First, the variables are initialized for 𝑁= 0: (a) The semi-minor and semi-major axes lengths of the contact ellipse: 𝑎0=𝑎𝑜 and 𝑏0=𝑏𝑜, and the equivalent radii: 𝑅𝑒𝑞,𝑥 = (1∕𝑅𝐴,𝑥 + 1∕𝑅𝐵,𝑥)−1 and 𝑅𝑒𝑞,𝑦 = (1∕𝑅𝐴,𝑦 + 1∕𝑅𝐵,𝑦)−1. (b) The initial normal elastic deformation [108]: 𝑢0=𝑃∕(2𝐸∗√𝑎0𝑏0), where 𝐸∗= ( (1 − 𝜈2 𝐴)∕𝐸𝐴+ (1 − 𝜈2 𝐵)∕𝐸𝐵)−1. (c) Constant average pressure 𝑝0 is assumed over the contact area: 𝑝0=𝑃∕(𝜋𝑎0𝑏0). (d) The total wear depth is 𝑤𝑡𝑜𝑡𝑎𝑙,0=𝑤0+𝑢0, being 𝑤0= 0. (e) Finally, the kinematic variables are initialized as: 𝑉𝐴,0= 2𝜋𝑅𝐴,𝑥,0𝜔𝐴, 𝑉𝐵,0=𝑉𝐴,0(1 + 𝜉𝑥) and 𝑠0=|𝑉𝐴,0−𝑉𝐵,0|. (II) For revolution 𝑁+ 1, compute: (a) The wear depth: 𝑤𝑁+1 =𝑤𝑁+𝑖𝑤𝑝𝑁(2𝑎𝑁𝑠𝑁∕𝑉), where 𝑉= (𝑉𝐴,𝑁 +𝑉𝐵,𝑁 )∕2. (b) The semi-contact width: 𝑏𝑁+1 =√2𝑅𝑒𝑞,𝑦𝑤𝑡𝑜𝑡𝑎𝑙,𝑁 − (𝑤𝑡𝑜𝑡𝑎𝑙,𝑁 )2. (c) The curvature radius: 𝑅𝑒𝑞,𝑥,𝑁+1 = (1∕𝑅𝐴,𝑥,𝑁+1 + 1∕𝑅𝐵,𝑥 )−1, where 𝑅𝐴,𝑥,𝑁+1 =𝑅𝐴,𝑥,𝑁 − (𝑤𝑁+1 −𝑤𝑁). (d) The semi-contact width: 𝑎𝑁+1 =𝜅√4𝑃 𝑅𝑒𝑞,𝑥,𝑁+1∕(2𝜋𝑏𝑁+1 𝐸∗), where 𝜅=𝜋∕4 is the correction factor to adapt the equation to an elliptical area. (e) The normal pressure: 𝑝𝑁+1 = (𝜋∕4)√𝑃 𝐸∗∕(2𝜋𝑏𝑁+1𝑅𝑒𝑞,𝑥,𝑁+1). (f) The normal elastic deformation: 𝑢𝑁+1 =𝑃∕(2𝐸∗√𝑎𝑁+1𝑏𝑁+1 ). (g) The wear depth: 𝑤𝑡𝑜𝑡𝑎𝑙,𝑁+1 =𝑤𝑁+1 +𝑢𝑁+1 . (h) Finally, update the slip velocity: 𝑠𝑁+1 =|𝑉𝐴,𝑁+1 −𝑉𝐵,𝑁+1 |, where 𝑉𝐴,𝑁+1 = 2𝜋𝑅𝐴,𝑥,𝑁+1𝜔𝐴 and 𝑉𝐵,𝑁+1 =𝑉𝐴,𝑁+1(1 + 𝜉𝑥). (III) This process will finish when 𝑁+ 1 = 𝑁𝑚𝑎𝑥. If 𝑘+ 1 < 𝑁𝑚𝑎𝑥, 𝑁=𝑁+ 1 and go to (II). being 𝑙=√(𝑥′−𝑥)2+ (𝑦′−𝑦)2 and 𝐺, 𝜈 and 𝐾 are the material parameters defined as 1 𝐺=1 2(1 𝐺(𝐴)+1 𝐺(𝐵)),𝜈 𝐺=1 2(𝜈(𝐴) 𝐺(𝐴)+𝜈(𝐵) 𝐺(𝐵)), 𝐾=𝐺 4(1−2𝜈(𝐴) 𝐺(𝐴)−1−2𝜈(𝐵) 𝐺(𝐵)). (45) The terms of Eq. (43) can be also expressed as: 𝐀(𝐱,𝐱′) = 𝐀(𝐱′−𝐱), i.e., indicating that the influence coefficients depend on the relative positions of the two surface points 𝐱 and 𝐱′. Appendix B. Validation of the wear depth computing scheme This appendix presents the validation of the proposed wear computing scheme presented in Section 7.5. For this purpose, we solve the twin-disc tribometer system presented in [66], where upper disc 𝛺𝐴 of radii 𝑅𝐴,𝑥 = 32.5 mm and 𝑅𝐴,𝑦 = 125 mm is subjected to a normal force 𝑃= 300𝑁 and to a rotation velocity 𝜔𝐴= 300 rev/min, being the Creep = 0.5% (𝜉𝑥= −0.005). The radii of lower disc 𝛺𝐵 are 𝑅𝐵,𝑥 = 32.3 mm and 𝑅𝐵,𝑦 = ∞. In this problem, both discs have the same material properties, i.e., elasticity modulus 𝐸(𝐴)=𝐸(𝐵)= 208×103 MPa and the Poisson’s ratio 𝜈(𝐴)=𝜈(𝐵)= 0.30, and isotropic friction and wear laws are considered: 𝜇= 0.6 and 𝑖𝑤= 2.0⋅10−6 MPa−1. Hegadekatte et al. [65,107] proposed a theoretical scheme, i.e., global incremental wear model (GIWM) – see Table 1 –, to compute the evolution of the maximum wear depth (𝑤𝑁) as a function of the number of revolutions (𝑁). Therefore, in Fig. 15, we can compare the evolution of the maximum wear depth – with 𝑁 – that we have obtained with the solution obtained with the GIWM. We can see a good agreement between both results, what allows us to validate our model. Appendix C. Surface contact variables evolution This Appendix shows some additional results obtained in the twindisc system problem presented in Section 8.2. In Fig. 16, we can see the evolution of wear depth distributions in the contact zone with the number of revolutions, i.e., 𝑁= {50,250,1000,4000,15000}, for 𝛽= {0◦,45◦,90◦}. As we have mentioned in Section 8.2, 𝛽 angle defines the angle between the tribological 𝑒1 axis and the rolling direction (𝑥-axis). According to the friction and wear coefficients selected for the – orthotropic – friction (𝜇1= 0.60 and 𝜇2= 0.30,) and wear (𝑖1= 2.0⋅10−6 MPa−1 and 𝑖2= 10−6 MPa−1) laws, when 𝛽 increases its value, the effective friction and wear coefficients in the rolling direction decreases. For this reason, the wear depth values observed in Fig. 16 for 𝛽= 0◦ are higher than the values obtained for 𝛽= 45◦ or 𝛽= 90◦. In Fig. 17, we can see the evolution of dimensionless normal contact traction (𝑝𝑛∕𝑝𝑜) with the number of revolutions, i.e., 𝑁= {50,250,1000,4000,15000}, for 𝛽= {0◦,45◦,90◦}. We can see that the normal contact pressure module evolves very similar for 𝛽= {0◦,45◦,90◦}. However, the contact width – perpendicular to the rolling direction – for 𝛽= 0◦ is the greatest one, since, according to Fig. 16, the severer wear was obtained for 𝛽= 0◦. The semi contact width in the rolling direction (𝑎) decreases with 𝑁, meanwhile the semi contact width perpendicular to the rolling direction (𝑏) increases with 𝑁. Moreover, after 𝑁= 15000 revolutions, the maximum normal contact pressure location moves from 𝑦= 0 plane – i.e., the center of the contact zone – when 𝑁= 0, to 𝑦= ±𝑏 plane – i.e., the contact zone limits –. Finally, Fig. 18 presents the normalized tangential slip velocity module (‖𝐬𝑡‖𝑖∕𝑉) evolution with the number of revolutions, i.e., 𝑁= {50,250,1000,4000,15000}, for 𝛽= {0◦,45◦,90◦}. We can see that, in International Journal of Mechanical Sciences 294 (2025) 110195 16 J.M. Juliá and L. Rodríguez-Tembleque Fig. 16. Wear depth evolution, expressed in millimeters, on the 𝑥𝑦 plane with the number of revolutions, i.e., 𝑁= {50,250,1000,4000,15000}, for the following values of the angle between the tribological axes and the coordinate system ones: (a) 𝛽= 0◦, (b) 𝛽= 45◦ and (c) 𝛽= 90◦. The 𝑥 and 𝑦 axis are dimensionless as they are expressed relative to the initial contact semi-width (𝑎𝑜 and 𝑏𝑜), respectively. Fig. 17, the normal contact pressure module evolves very similar for 𝛽= {0◦,45◦,90◦}. However, in Fig. 18, the greatest values of the normalized tangential slip velocity module are obtained for 𝛽= 0◦ (i.e., for the highest values of the effective friction and wear intensity coefficients in the rolling direction). This evolution of the normalized tangential slip velocity module and the analysis of Eqs. (16) and (19) explain why the highest values of wear depth are obtained in each rotation for 𝛽= 0◦ (i.e., for the highest values of the effective friction and wear intensity coefficients in the rolling direction). Data availability Data will be made available on request. International Journal of Mechanical Sciences 294 (2025) 110195 17 J.M. Juliá and L. Rodríguez-Tembleque Fig. 17. Dimensionless normal contact traction (𝑝𝑛∕𝑝𝑜) evolution colorgreen on the 𝑥𝑦 plane with the number of revolutions, i.e., 𝑁= {50,250,1000,4000,15000}, for the following values of the angle between the tribological axes and the coordinate system ones: (a) 𝛽= 0◦, (b) 𝛽= 45◦ and (c) 𝛽= 90◦. All the results are non-dimensional: the normal contact traction is presented relative to the initial maximum contact pressure (𝑝𝑜) and the 𝑥 and 𝑦 axis are expressed relative to the initial contact semi-width (𝑎𝑜 and 𝑏𝑜), respectively. International Journal of Mechanical Sciences 294 (2025) 110195 18 J.M. Juliá and L. Rodríguez-Tembleque Fig. 18. Normalized tangential slip velocity module (‖𝐬𝑡‖∕𝑉) evolution on the 𝑥𝑦 plane with the number of revolutions, i.e., 𝑁= {50,250,1000,4000,15000}, for the following values of the angle between the tribological axes and the coordinate system ones: (a) 𝛽= 0◦, (b) 𝛽= 45◦ and (c) 𝛽= 90◦. The 𝑥 and 𝑦 axis are expressed relative to the initial contact semi-width (𝑎𝑜 and 𝑏𝑜), respectively. International Journal of Mechanical Sciences 294 (2025) 110195 19 J.M. Juliá and L. Rodríguez-Tembleque References [1] Johnson KL. Contact mechanics. Cambridge: Cambridge University Press; 1987. [2] Harris TA, Kotzalas MN. Advanced concepts of bearing technology. Boca Raton, Florida: Taylor & Francis Group, LLC; 2007. [3] Guo Y, Keller J. Rolling element bearing dynamics in wind turbines. NREL/PR, National Renewable Energy Laboratory; 2018. [4] Doll GL. Rolling bearing tribology: Tribology and failure modes of rolling element bearings. Elsevier Science; 2022. [5] Li W, Tan Y, Tao Y, Bai D. Wear analysis of rolling bearing contact surfaces based on finite element method. Phys Scr 2023;99:015916. [6] Kunzelmann B, Rycerz P, Xu Y, Arakere NK, Kadiric A. Prediction of rolling contact fatigue crack propagation in bearing steels using experimental crack growth data and linear elastic fracture mechanics. Int J Fatigue 2023;168:107449. [7] Garg VK, Dukkipati RV. Dynamics of railway vehicle systems. Canada: Academic Press; 1984. [8] Kalker J. Three-dimensional elastic bodies in rolling contact. The Netherlands: Kluwer Academic Publisher; 1990. [9] Wu B, Wang W, Pan J, Hu Y, Xu R, Ye D, Yan W, Zhang R. Study on corrugated wear on high-speed railways based on an improved finite element model of wheel-rail rolling contact. Tribol Int 2023;179:108199. [10] Reynolds O. On rolling friction. Phil Trans R Soc 1876;166:155–74. [11] Hertz H. Über die Berührung fester elastischer Körper. J Für Die Reine Und Angew Math 1882;92:156–71. [12] Carter FW. On the action of a locomotive driving wheel. Proc R Soc 1926;112:151–7. [13] Johnson KL. The effect of spin upon the rolling motion of an elastic sphere upon a plane. J Appl Mech 1985;25:332–8. [14] Johnson KL. The effect of a tangential force upon the rolling motion of an elastic sphere upon a plane. J Appl Mech 1985;25:339–46. [15] Vermeulen PJ, Johnson KL. Contact of nonspherical elastic bodies transmitting tangential forces. J Appl Mech 1964;31:338–40. [16] Bentall R, Johnson K. Slip in the rolling contact of two dissimilar elastic rollers. Int J Mech Sci 1967;9:389–404. [17] Nowell D, Hills D. Tractive rolling of dissimilar elastic cylinders. Int J Mech Sci 1988;30:427–39. [18] Kalker JJ. On the rolling contact of two elastic bodies in the presence of dry friction. Wear 1968;11:1–303. [19] Kalker JJ. The tangential force transmitted by two elastic bodies rolling over each other with pure creepage. Wear 1968;2:421–30. [20] Kalker JJ. A fast algorithm for the simplified theory of rolling contact. Veh Syst Dyn 1982;11:1–13. [21] Kalker JJ. Wheel-rail rolling contact theory. Wear 1991;144:243–61. [22] Kalker JJ, Jacobson B. Rolling contact phenomena. Springer Vienna; 2000. [23] González JA, Abascal R. An algorithm to solve coupled 2D rolling contact problems. Internat J Numer Methods Engrg 2000;49(9):1143–67. [24] González JA, Abascal R. Solving 2D transient rolling contact problems using the BEM and mathematical programming techniques. Internat J Numer Methods Engrg 2002;53(4):843–75. [25] Abascal R, Rodríguez-Tembleque L. Steady-state 3D rolling-contact using boundary elements. Commun Numer Math Eng 2007;33:905–20. [26] Rodríguez-Tembleque L, Abascal R. A 3D FEM-BEM rolling contact formulation for unstructured meshes. Int J Solids Struct 2010;47:330–53. [27] Santeramo M, Putignano C, Vorlaufer G, Krenn S, Carbone G. Viscoelastic steady-state rolling contacts: A generalized boundary element formulation for conformal and non-conformal geometries. J Mech Phys Solids 2023;171:105129. [28] Güler MA, Alinia Y, Adibnazari S. On the rolling contact problem of two elastic solids with graded coatings. Int J Mech Sci 2012;64:62–81. [29] Güler MA, Adibnazari S, Alinia Y. Tractive rolling contact mechanics of graded coatings. Int J Solids Struct 2012;49(6):929–45. [30] Güler MA, Alinia Y, Adibnazari S. On the contact mechanics of a rolling cylinder on a graded coating. part 2: Numerical results. Mech Mater 2013;66:134–59. [31] Alinia Y, Güler MA, Adibnazari S. On the contact mechanics of a rolling cylinder on a graded coating. part 1: Numerical results. Mech Mater 2014;68:207–16. [32] Xi Y, Almqvist A, Shi Y, Mao J, Larsson R. A complementarity problem–based solution procedure for 2D steady-state rolling contacts with dry friction. Tribol Trans 2016;59(6):1031–8. [33] Wang Z, Jin X, Keer LM, Wang Q. A numerical approach for analyzing three-dimensional steady-state rolling contact including creep using a fast semi-analytical method.. Tribol Trans 2012;55(4):446–57. [34] Xi Y, Almqvist A, Shi Y, Mao J, Larsson R. Linear complementarity framework for 3D steady-state rolling contact problems including creepages with isotropic and anisotropic friction for circular hertzian contact. Tribol Trans 2017;60(5):832–44. [35] Meshcheryakova AR, Goryacheva IG. Stress state of elastic bodies with an intermediate layer in rolling contact with slip. Phys Mesomech 2021;23:91–101. [36] Xi Y, Björling M, Almqvist A. A numerical model for solving three-dimensional rolling contact problems with elastic coating layers. Tribol Lett 2021;69(4):139. [37] Xi Y, Li B, Almqvist A. Semi-analytical model for 3D multilayered rolling contact problems with different creepage combinations. J Tribol 2024;146:044103. [38] Manyo EY. Modélisation avancée du contact pneu-chaussée pour l’étude des dégradations des chaussées en surface (Ph.D thesis) (Ph.D. thesis), Université de Limoges; 2019. [39] Wallace E, Chaise T, Nelias D. Rolling contact on a viscoelastic multi-layered half-space. Int J Solids Struct 2022;239–240:111388. [40] Wallace ER, Chaise T, Duval A, Nelias D. Transient tractive rolling contact between elastically dissimilar and multi-layered bodies. Int J Solids Struct 2023;265–266:112124. [41] Yu Y, Suh J. Numerical analysis of three-dimensional thermo-elastic rolling contact under steady-state conditions. Friction 2021;1–15. [42] Blanco-Lorenzo J, Santamaria J, Vadillo EG, Correa N. On the influence of conformity on wheel-rail rolling contact mechanics. Tribol Int 2016;103:647–67. [43] Blanco-Lorenzo J, amd J. Santamaria SL, Meehan PA, Vadillo EG. Frictional contact analysis in a spherical roller bearing. J Comput Des Eng 2022;10:139–59. [44] Alonso A, Giménez JG. Non-steady state modelling of wheel-rail contact problem for the dynamic simulation of railway vehicles. Veh Syst Dyn 2008;46:179–96. [45] Guiral A, Alonso A, Baeza L, Giménez JG. Non-steady state modelling of wheel–rail contact problem. Veh Syst Dyn 2013;51:91–108. [46] Baeza L, Giner-Navarro J, Thompson DJ, Monterde J. Eulerian models of the rotating flexible wheelset for high frequency railway dynamics. J Sound Vib 2019;449:300–14. [47] Baeza L, Bruni S, Giner-Navarro J, Liu B. A linear non-Hertzian unsteady tangential wheel-rail contact model. Tribol Int 2023;181:108345. [48] Brhane AG, Mekonone ST. Numeric simulation of steel twin disc system under rollingsliding contact. Tribol Mater 2023;2(4):181–8. [49] Sharma A, Sadeghi F. Effects of fretting wear on rolling contact fatigue. Tribol Int 2024;192:109204. [50] Zhao X, Li Z. The solution of frictional wheel-rail rolling contact with a 3D transient finite element model: Validation and error analysis. Wear 2011;271 (1–2):444–52. [51] Zhao X, Li Z. A three-dimensional finite element solution of frictional wheel-rail rolling contact in elasto-plasticity. Proc Inst Mech Eng J 2015;229(1):86–100. [52] Yang Z, Deng X, Li Z. Numerical modeling of dynamic frictional rolling contact with an explicit finite element method. Tribol Int 2019;129:214–31. [53] Sadeghi F, Jalalahmadi B, Slack TS, Raje N, Arakere NK. A review of rolling contact fatigue. J Tribol 2009;131(4):041403. [54] Khoramzad E, Nia SH, Casanueva C, Berg M. Impact of slip velocity-dependent friction coefficient on surface traction, wear, RCF and curve squeal noise prediction in wheel-rail contact. Int J Veh Mech Mobil 2025;1–19. [55] Wang HH, Wang WJ, Han ZY, Wang Y, Ding HH, Lewis R, Lin Q, Liu QY, Zhou ZR. Wear and rolling contact fatigue competition mechanism of different types of rail steels under various slip ratios. Wear 2023;522:204721. [56] Wang K, Lai J, Xu J, Liao T, Wang P, Chen R, Qian Y, Li L, Ma X. Multiscale analysis of wheel-rail rolling contact wear and damage mechanisms using molecular dynamics and explicit finite elements. Tribol Int 2023;185:108574. [57] Shabana A, Zaazaa KE, Escalona JL, Sany JR. Development of elastic force model for wheel/rail contact problems. J Sound Vib 2004;269:295–325. [58] García-Agúndez A, García-Vallejo D, Freire E. Linearization approaches for general multibody systems validated through stability analysis of a benchmark bicycle model. Nonlinear Dynam 2021;103:557–80. [59] García-Agúndez A, García-Vallejo D, Freire E. Analytical and numerical stability analysis of a toroidal wheel with nonholonomic constraints. Nonlinear Dynam 2024;112:2453–76. [60] Olofsson U, Andersson S, Björklund S. Simulation of mild wear in boundary lubricated spherical roller thrust bearing. Wear 2000;241:180–5. [61] Jendel T. Prediction of wheel profile wear-comparisons with field measurements. Wear 2002;253:89–99. [62] Enblom R, Berg M. Simulation of railway wheel profile development due to wear-influence of disc braking and contact environment. Wear 2005;258:1055–63. [63] Telliskivi T. Simulation of wear in a rolling-sliding contact by a semi-Winkler model and the Archard’s wear law. Wear 2004;256:817–31. [64] Telliskivi T, Olofsson U. Wheel-rail wear simulation. Wear 2004;257:1145–53. [65] Hegadekatte V, Kurzenhäuser S, Kraft OA. Predictive modeling scheme for wear in tribometers. Tribol Int 2008;41:1020–31. [66] Rodríguez-Tembleque L, Abascal R, Aliabadi MH. A boundary element formulation for wear modeling on 3D contact and rolling-contact problems. Int J Solids Struct 2010;47:2600–12. [67] Rodríguez-Tembleque L, Abascal R, Aliabadi MH. Anisotropic wear framework for 3D contact and rolling problems. Comput Methods Appl Mech Engrg 2012;241:1–19. [68] Liu B, Bruni S, Lewis R. Numerical calculation of wear in rolling contact based on the archard equation: Effect of contact parameters and consideration of uncertainties. Wear 2022;490–491:204188. International Journal of Mechanical Sciences 294 (2025) 110195 20 J.M. Juliá and L. Rodríguez-Tembleque [69] Afridi AH, Zhu H, Camacho ET, Deng G, Li H. Numerical modeling of rolling contact fatigue cracks in the railhead. Eng Fail Anal 2023;143:106838. [70] Nezhad MS, Larsson F, Kabo E, Ekberg A. Numerical prediction of railhead rolling contact fatigue crack growth. Wear 2023;530–531:205003. [71] Polancec T, Lesicar T, Tonkovic Z, Glodez S. Modelling of rolling-contact fatigue pitting phenomena by phase field method. Wear 2023;532–533:205068. [72] Zhou Y, Zhu C, Song C, Chen Y, Zhou M. The effect of friction on surface crack initiation in rolling contact fatigue considering damage accumulation. Fatigue Fract Eng Mater Struct 2023;46:3487–500. [73] Wen B, Tao G, Wen Z. Prediction of locomotive wheel wear evolution considering thermo-mechanical coupling. Wear 2025;205805. [74] Hilal M, Bouallala M, Ouafik Y. Variational and numerical approaches to contact problems with Coulomb friction in thermoviscoelasticity. Math Methods Appl Sci 2025. [75] Peela A, Mikitisin A, Steinweg F, Janitzky T, Schwedt A, Broeckmann C, Mayer J. Influence of subsurface non-metallic inclusions on rolling contact fatigue behavior in bearing steels - experimental and numerical investigations. Wear 2025;566–567:205767. [76] Ravi G, Waele W, Nikolic K, Petrov R, Hertelé S. Numerical modelling of rolling contact fatigue damage initiation from non-metallic inclusions in bearing steel. Tribol Int 2023;180:108290. [77] Fu P, Zhao J, Zhang X, Miao H, Wen Z, Kang G, Kan Q. Three-dimensional tractive rolling contact analysis of functionally graded coating-substrate systems with interfacial imperfection and frictional anisotropy. Compos Struct 2023;307:116671. [78] Kalker JJ. The computation of three-dimensional rolling contact with dry friction. Internat J Numer Methods Engrg 1979;14(9):1293–307. [79] Wayne Chen W, Jane Wang Q. A numerical model for the point contact of dissimilar materials considering tangential tractions. Mech Mater 2008;40(11):936–48. [80] Bazrafshan M, de Rooij M, Schipper D. On the role of adhesion and roughness in stick-slip transition at the contact of two bodies: A numerical study. Tribol Int 2018;121:381–8. [81] Zhao Y, Liu HC, Morales-Espejel GE, Venner CH. Response of elastic/viscoelastic layers on an elastic half-space in rolling contacts: Towards a new modeling approach for elastohydrodynamic lubrication. Tribol Int 2023;186:108545. [82] Juliá JM, Rodríguez-Tembleque L. Subsurface stress evolution under orthotropic wear and frictional contact conditions. Int J Mech Sci 2022;234:107695. [83] Alart P, Curnier A. A mixed formulation for frictional contact problems prone to Newton like solution methods. Comput Methods Appl Mech Engrg 1991;92(3):353–75. [84] Juliá JM, Rodríguez-Tembleque L. Wear and subsurface stress evolution in a half-space under cyclic flat-punch indentation. Lubricants 2023;11(6). [85] Rodríguez-Tembleque L, Abascal R. Fast FE-BEM algorithms for orthotropic frictional contact. Internat J Numer Methods Engrg 2013;94:687–707. [86] Michalowski R, Mróz Z. Associated and non-associated sliding rules in contact friction problems. Arch Mech Stosow 1978;30:259–76. [87] Mróz Z, Stupkiewicz S. An anisotropic friction and wear model. Int J Solids Struct 1994;31:1113–31. [88] Mróz Z, Kucharski S, Paczelt I. Anisotropic friction and wear rules with account for contact state evolution. Wear 2018;396:1–11. [89] Rabinowicz E. Friction and wear of materials. John Wiley & Sons; 1965. [90] McColl IR, J D, Leen SB. Finite element simulation and experimental validation of fretting wear. Wear 2004;256:1114–27. [91] Rodríguez-Tembleque L, Abascal R, Aliabadi MH. A boundary element formulation for 3D fretting-wear problems. Eng Anal Bound Elem 2011;35:935–43. [92] Paczelt I, Kucharski S, Mróz Z. The experimental and numerical analysis of quasi-steady wear processes for a sliding spherical indenter. Wear 2012;274–275:127–48. [93] Paczelt I, Mróz Z. Solution of wear problems for monotonic and periodic sliding with p-version of the finite element method. Comput Methods Appl Mech Engrg 2012;249–252:75–103. [94] Stupkiewicz S. An ALE formulation for implicit time integration of quasi-steady-state wear problems. Comput Methods Appl Mech Engrg 2013;260:130–42. [95] Cavalieri JF, Cardona A. Three-dimensional numerical solution for wear prediction using a mortar contact algorithm. Internat J Numer Methods Engrg 2013;96:467–86. [96] Lengiewicz J, Stupkiewicz S. Efficient model of evolution of wear in quasi-steady-state sliding contacts. Wear 2013;303:611–21. [97] Willner K. Fully coupled frictional contact using elastic halfspace theory. J Tribol 2008;130(3). [98] Pohrt R, Li Q. Complete boundary element formulation for normal and tangential contact problems. Phys Mesomech 2014;17:334–40. [99] Wang QY, Bathias C, Kawagoishi N, Chen Q. Effect of inclusion on subsurface crack initiation and gigacycle fatigue strength. Int J Fatigue 2002;24:1269–74. [100] Joli P, Feng Z-Q. Uzawa and newton algorithms to solve frictional contact problems within the bi-potential framework. Internat J Numer Methods Engrg 2008;73:317–30. [101] Ning P, Feng Z-Q, Quintero JAR, Zhou Y-J, Peng L. Uzawa algorithm to solve elastic and elastic–plastic fretting wear problems within the bipotential framework. Comput Mech 2018;62(6):1327–41. [102] Zaazaa KE, Schwab AL. Review of joost Kalker’s wheel-rail contact theories and their implementation in multibody codes. In: Proceedings the ASME international design engineering technical conferences and computers and information in engineering conference 2009, DETC2009, PART C. 2009, p. 1889–900. [103] Erena D, Vázquez J, Navarro C, Domínguez J. A fretting fatigue model based on self-steered cracks. Theor Appl Fract Mech 2022;117:103144. [104] Rangel D, Erena D, Vázquez J, Araújo JA. Prediction of initiation and total life in fretting fatigue considering kinked cracks. Theor Appl Fract Mech 2022;119:103345. [105] Fang C, Peng Y, Guan Y, Zhou W, Gao G, Meng X. A new numerical method for the tribo-dynamic analysis of cylindrical roller bearings. Nonlinear Dynam 2023;111:11275–95. [106] Rudnytskyj A, Vorlaufer G, Leimhofer J, Jech M, Gachot C. Estimating the real contact area in lubricated hot rolling of aluminium. Tribol Int 2023;180:108283. [107] Hegadekatte V, Huber N, Kraft OA. Modeling and simulation of wear in a pin on disc tribometer. Trib. Lett 2006;24:51–60. [108] Oliver WC, Pharr GM. An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments.. J Mater Res 1992;7:1564–83. International Journal of Mechanical Sciences 294 (2025) 110195 21