Numerical Modeling of Wear in a Thrust Roller Bearing under Mixed Elastohydrodynamic Lubrication

: Increasing efforts to reduce frictional losses and the associated use of low-viscosity lubricants lead to machine elements being operated under mixed lubrication. Consequently, wear effects are also gaining relevance. Appropriate numerical modeling and predicting wear in a reliable manner offers new possibilities for identifying harmful operating conditions or for designing running-in procedures. However, most previous investigations focused on simplified model contacts and the wear behavior of application-oriented contacts is relatively underexplored. Therefore, the contribution of this paper was to provide a numerical procedure for studying the wear evolution in the mixed elastohydrodynamically lubricated (EHL) roller/raceway contact by coupling a finite element method (FEM)-based 3D EHL model with surface topography changes following a local Archard-type wear model, a Greenwood–Williamson-based load-sharing approach and the Sugimura surface adaption model. This study applied the operating conditions of an 81212 thrust roller bearing, considering realistic geometry and locally varying velocities. The calculated wear profiles in the raceway featured asymmetries, which were in good agreement with the experimental results reported in the literature and could be correlated with the velocity and slip distribution. In addition, the effects of speed, load and oil viscosity were investigated by means of four load cases and two lubricants, demonstrating the broad range of applying the numerical approach. Conceptualization, A.W. and M.M.; methodology, A.W. and M.M.; software, A.W.; formal analysis, A.W. and M.M.; writing—original draft preparation, A.W. and M.M.; writing—review and editing, A.W., M.M., S.T. and S.W.; visualization, A.W. and M.M.; supervision, S.T. and S.W. All authors have read and agreed to the published version of the manuscript.


Introduction
According to Holmberg and Erdemir [1], tribological contacts account for about 23% of the world's total energy consumption, whereas overcoming friction contributes to about 20% and remanufacturing worn parts and spare equipment due to wear and wear-related failures amount to 3%. By utilizing new surface, material and lubrication technologies for friction reduction and wear protection in vehicles, machines and other equipment, energy losses could be reduced by up to 40% [1]. However, wear gets increasingly promoted by the expanding use of low-viscosity lubricants and the tendency of operating in the friction optima in the mixed lubrication regime [2]. Moreover, wear can be considered critical as it can lead to catastrophic failures and breakdowns, thus having a negative impact on productivity and costs. Therefore, appropriate numerical modeling over different scales and predicting wear in a reliable manner offers new possibilities for understanding and redesigning tribological systems in machine elements or engine components. Consequently, the comprehension and simulation of wear in tribo-contacts were the subject of several studies ranging from boundary to hydro-(HL) or elastohydrodynamically lubricated (EHL) contacts. Thereby, Podra et al. [3,4] numerically investigated the wear of a dry pin-on-disc and a conical spinning contact based upon the Archard [5] wear model. The calculation of the contact pressure was accomplished using the finite element method (FEM) and by means of the Winkler surface model [6]. Hegadekatte et al. [7,8] also examined the wear of a pin-on-disc contact based upon the FEM and the Archard model, however considered local wear by shifting the nodes of an adaptive mesh in the normal direction. Andersson et al. [9] investigated the wear of a dry sphere on flat contact. The contact pressure was calculated by a discrete convolution and fast Fourier transformation (DC-FFT) method utilizing the half-space theory and assuming linear elastic-perfectly plastic material behavior as described by Liu et al. [10]. Ashraf et al. [11] implemented a FEM-based wear simulation for a dry composite alloy contact using a free-mesh. Sfantos et al. [12] suggested a boundary element formulation for three-dimensional dry sliding wear based on Archard's model, which was applied to a pin-on-disc contact as well as to a hip arthroplasty wear problem. As aforementioned approaches did not take deterministic surface topography into account, Terwey et al. [13,14] implemented a contact and wear model based upon the half-space theory for boundary lubricated thrust roller bearings, considering surface roughness. Thereby, the wear coefficient of Archard's law was determined using continuum damage mechanics (CDM). A comparison with test rig experiments revealed that the severe wear regime could not be satisfactorily described, whereas the good agreement of numerical prediction and test results was achieved for mild wear.
In addition to previously mentioned methods for the numerical wear calculation of dry and boundary lubricated contacts, several approaches for the mixed hydro-or elastohydrodynamically lubricated regime can be found in the literature. Zhu et al. [15] suggested an approach for the numerical wear calculation in lubricated contacts based upon a deterministic mixed elastohydrodynamic lubrication model [16,17]. Thereby, the surface topography was directly incorporated into the film thickness equation and the wear volume was determined by means of Archard's wear model. Lorentz et al. [18] developed a deterministic mixed lubrication micro-scale model considering two rough rubbing bodies, an adhesion model, heat generation as well as a lubrication domain. Reichert et al. [19] extended the model by a wear calculation, whereas the wear coefficient for an Archard type wear model was determined based upon the Johnson-Cook damage law [20]. An stochastic approach for wear calculation in lubricated EHL line contacts based upon a load-sharing concept and a CDM-based ARCHARD type wear model [21], considering simplified thermo-elastohydrodynamic analysis [22], was introduced by Beheshti and Khonsari [23]. Hao et al. [24] proposed a computational fluid dynamics (CFD)-based wear calculation of a cylinder sliding over a ring in the mixed lubrication regime. Thereby, hydrodynamics was stochastically coupled with the asperity contact pressure, which was calculated according to an elastic-plastic contact model [25]. The cylinder geometry profile was updated stepwise following Archard's wear law. Zhang et al. [26] investigated the wear and roughness evolution in a mixed lubricated line contact based upon a stochastic approach. The finite difference method (FDM) was used for the EHL simulation and microhydrodynamics were considered by flow factors according to Patir and Cheng [27,28]. The asperity contact pressure was determined by means of the Kogut-Etsion model [29][30][31]. Moreover, the change in surface height probability density function (PDF) was calculated in accordance to Sugimura and Kimura [32][33][34]. Local wear was computed by Archard's wear law using the asperity contact pressure. Moreover, Zhang et al. [35] analyzed the profile modification of a cylindrical roller in the mixed lubrication regime while taking thermal effects into account.
In brief, a wide range of approaches for the numerical calculation of wear processes have been developed in recent decades. Some were limited to dry contacts while others considered mixed HL or EHL contacts. Moreover, deterministic [16,17,36] and stochastic [37] lubrication models could be distinguished. The latter featured advantages regarding the computational effort due to a reduced required discretization. However, the wear behavior of application-oriented contacts with realistic geometries and locally varying velocities as encountered in actual machine elements is relatively underexplored due to the numerical complexity. Thus, most investigations focused on simplified model contacts. Furthermore, only a few selected load cases and lubricants were usually considered. For these reasons, this paper aims at providing a numerical procedure to predict the wear evolution in the mixed EHL roller/raceway contact of a thrust roller bearing. Moreover, the influence of different load-cases and lubricants on the wear behavior was studied in detail.

Numerical Modeling
For the numerical wear calculation of a rolling element and raceway in a thrust roller bearing, a full-system finite element modeling (FEM) approach for EHL contacts was expanded by a mixed lubrication model as well as a local Archard type wear model to describe the surface topography evolution. The calculation model was part of the computation software TriboFEM. The investigated tribo-system as well as the main physical and numerical characteristics are thoroughly described and reasoned below.

Load Cases, Material and Lubricant Properties
In order to allow qualitative comparability with the experimental results from literature, the wear calculations were carried out by means of a thrust roller bearing 81212. Four load cases with the speeds, simulated test durations, axial loads and initial Hertzian pressures in the roller/raceway contact (assuming a pure line contact, i.e., without profiling) according to Table 1 were studied. Thus, each washer surface element experienced a constant number of 2,850,000 overrollings. Moreover, a realistic velocity distribution was taken into account, whereby u1 was considered as the relative velocity of the washer in the rolling x direction: The roller maintained a constant peripheral speed of: Perpendicular to the rolling direction, the washer had a relative velocity of: whereas the roller had no velocity component in the y direction: Underlying geometry and coordinates are illustrated in Figure 1. It should be noted that the Cartesian coordinate system depicted in Figure 1 was rotated with the center of the considered rolling element. Thus, the resulting sliding velocity could be written as A Young's modulus of 210,000 MPa and a Poisson's ratio of 0.3 were considered as the mechanical properties of the rollers and washers. Moreover, the surface height probability density function was assumed to be Gaussian with an equivalent standard deviation of σeq = 0.15 µm. Two mineral reference oil types with rather low viscosities were studied, for which the relevant lubricant properties are summed up in Table 2. Density and viscosity were considered to be pressure-dependent following the equations from Dowson and Higginson [38] and Roelands [39], respectively: Shear-thinning and thermal effects were neglected within the scope of this study and are subject of on-going research work.

Asperity Contact Model
An approach to calculate the real contact area and the load between an elastic rough surface and a rigid smooth plane based upon Hertzian-like deformations of asperities was first suggested by Greenwood and Williamson [40]. Since then, the model has been expanded, for example to take into account curved surfaces [41], two rough surfaces with misaligned asperities [42], elliptic paraboloid asperities [43] and anisotropic surfaces [44]. In addition, Chang et al. [25], Halling and Nuri [45], Zhao et al. [46] as well as Jamari and Schipper [47] developed analytical asperity contact models that allowed the consideration of elastic-plastic material properties. Moreover, Kogut and Etsion [29][30][31] as well as Jackson and Green [48,49] introduced FEM-based asperity contact models also taking into account elastic-plastic material behavior. Due to its broad applicability and simplicity, the Greenwood-Williamson model was used within the present study. It should be noted that all dimensionless values denoted by * were normalized by σ.
The asperity contact pressure in the dimensional form was given by Since the Greenwood-Williamson model used the mean of the summit heights as a reference and the mean of the surface heights was usually considered in EHL and by the applied Sugimura model, a transformation of the derived parameters was required. According to Nayak [50], Bush et al. [51] and McCool [44], the following relation between the standard deviation of asperity heights and surface heights applied: whereas the bandwidth parameter α was defined by The parameters m0 = σ 2 , m2 and m4 denote the zeroth, second and fourth spectral moments of a surface profile. Considering the density of summits: and the mean summit radius: led to the following relation between the standard deviation of the surface heights and the asperity heights: The distance between the mean height of the asperities and the mean height of the surface according to Bush et al. [51] was: Finally, the relation between the separation of a rough surface and a flat surface based on the asperity heights and surface heights reads: The relation between the PDF of the summit heights and the PDF of the surface heights followed: Taking into account two rough surfaces by the Greenwood-Williamson model is possible by calculating equivalent values for the standard deviation, the asperity radius and the asperity density according to Beheshti and Khonsari [52]:

EHL Modeling
In the context of higher loaded non-conformal contacts, CFD-based approaches, as suggested by Hartinger et al. [53,54], usually lead to high computational effort and instabilities at high loads. Consequently, numerical EHL modeling has frequently been done by applying the Reynolds equation [55] for hydrodynamics coupled with the solution of the elastic deformation. Therefore, iterative finite difference and half-space-based multilevel and multi-integration (MLMI) solvers, detailed in Venner and Lubrecht [56], are well established. Furthermore, Habchi et al. [57] introduced the full-system FEM-based approach to solve the system of fully coupled equations. This offered advantages in terms of computational effort and convergence behavior and most notably, enabled the use of commercially available software. Hence, the latter approach was used within the scope of this study.
Therefore, the Reynolds equation in a slightly modified form: velocity / Couette is solved in its weak form on the upper surface Ωc of the elastically deformed equivalent body ( Figure  2) to describe the lubricant's hydrodynamics. The first term describes the influence of the pressure gradient, while the second one accounts for the boundary velocities of the contacting bodies and the wedge-shape of the lubricant gap. Zero pressure boundary conditions were applied at the in-and outlet. A mass-conserving algorithm, as introduced by Marian et al. [58], dealt with cavitation effects. Therefore, the lubricant's density and viscosity were multiplied with the fractional film content. The latter was defined as the ratio of the lubricant layer to the gap height: with the penalty factor γ(ph), which is zero if ph < 0, otherwise γ(ph) = ξ, where ξ is a sufficiently high algebraic number. The Galerkin least squares (GLS) method (Hughes et al. [59]) and isotropic diffusion (Zienkiewicz et al. [60]) were utilized for the numerical stabilization.
Based upon the FEM, the elastic problem was calculated for one body Ω with the equivalent Young's modulus: by applying the linear elasticity equation: A free triangular mesh with a refinement in the contact center of the upper surface was applied and the body was provided with a zero-displacement boundary condition on the bottom and free boundary conditions on the sides. The dimensions of the cuboid were chosen large enough to exclude any effects on the calculation results [57]. Furthermore, a free triangular mesh with a refinement in the contact center of the upper surface was applied. Geometry and meshing are schematically illustrated in Figure 2.
describes the height of the separating fluid film in terms of the distance and the shape of the undeformed geometry as well as of the elastic deformation. The geometry of the rolling element was composed of a quadratic approximation of the radius in and the logarithmic profile perpendicular to the direction of motion: The load balance equation: ensures the equilibrium of forces, taking into account simultaneously occurring asperity contact (Section 2.1.) and hydrodynamic lubrication by a stochastic mixed lubrication [61].
In order to obtain a consistent model without regenerating and remeshing the computational domains at different load cases as well as to avoid instabilities due to memory variables with limited accuracy, all the solution variables were normalized to characteristic quantities of the Hertzian theory or, in the case of the lubricant, initial values: In addition, the computational domain was decreased in y direction by the factor scale in order to reduce the number of degrees of freedom and thus the computation time.

Wear Modeling
Widely used wear models were proposed by Kragelski [62] and Fleischer et al. [63]. The first developed a molecular mechanical fatigue-based wear model and the approach of the latter considered the frictional energy density. However, the most common model for calculating wear, originally intended for dry contacts, goes back to Archard [5,64] and Holm [65]. Although various other wear modeling approaches, based upon empirical relationships or limited to specific applications, exist in the literature, Archard's wear model was applied within the scope of this contribution for its generality and simplicity. Thereby, the total wear volume Vw was considered proportional to the normal force FN, the sliding distance s and the proportionality factor k, also known as the wear coefficient: To calculate the wear in a lubricated contact, the Archard wear model, which in its conventional form is valid for dry contacts, was modified assuming that wear could only occur on a solid asperity contact: As a result, the local wear depth was derived by , h x y k s p x y .
Since the wear coefficient depends on numerous factors, such as material properties, surface conditions and boundary layers, it was frequently determined experimentally and could vary between 10 −1 and 10 −15 mm 3 /Nm for metals [66,67]. As it can be assumed that even under the highest loads, some protecting boundary layer will be formed and that at least some lubricant will remain in the contact between the two opposing asperities, it seems appropriate in the context of the present study to apply a wear coefficient that is valid for the boundary lubrication regime in order to calculate the wear depth. In accordance to Czichos and Habig [68], a local wear coefficient of k = 6 × 10 -8 mm 3 /Nm was applied.

Surface Topography Model
The Sugimura model [32][33][34] was used to take into account the wear-induced change of statistical surface parameters within the asperity contact pressure calculation.
Assuming an arbitrary distribution of the initial surface heights, the upper proportion of the probability distribution S1 was removed, as schematically illustrated in Figure 3a, depending on the wear depth W and the mean wear-related height loss w : where: Thereby, ψ denotes the probability density function of the wear-induced height loss and could be derived from the wear particle size distribution: , where u, v, and w are the wear particle dimensions. Additionally, f(u,v,w) denotes the probability density function of the wear particle size. Assuming the proportionality of u and v to w results in: with the marginal density function of the wear particle density function: Following Kimura and Sugimura [34], an exponential distribution for the wear particle size, as depicted in Figure 3b, was assumed. Finally, the fraction of the removed surface PDF was then redistributed to the remaining density function accordingly (Figure 3a). The surface height probability density function at the time t + Δt was written by considering the shift of the mean height Δz0: where the shift of the mean height satisfies:

Overall Numerical Procedure
The global numerical solution scheme is presented in Figure 4. Based upon the initial parameters on bulk material, surface topography and lubricant rheology (see Section 2.1.), the asperity contact pressure in dependency of the mean surface distance was calculated according to the Greenwood-Williamson model as described in Section 2.2 using MATHWORKS MATLAB. Subsequently, the system of highly nonlinear EHL equations, see Section 2.3., was solved fully coupled utilizing COMSOL MULTIPHYSICS. For more fundamental aspects about FEM in EHL and further information about the implementation, the interested reader is referred to Habchi [69]. After convergence, the wear evolution was calculated based upon the resulting solid contact pressure using Archard's law as described in Section 2.4 and the profile was updated accordingly. Moreover, the PDF of the surface heights was adapted using the Sugimura model, see Section 2.5. These steps were also implemented in MATHWORKS MATLAB and repeated until the intended operating period of the thrust roller bearing was reached. It should be noted that each of the individual aforementioned sub-models (Section 2.2-2.5) used within the scope of this contribution were chosen based upon their similar degree of complexity and good implementability. However, they could in principle be replaced by any other or more sophisticated models.

Results and Discussion
The lubrication conditions of all four load cases and both oil types in terms of fluid film gap and total and solid contact pressure are depicted in Figure 5,6 and Figure 10-13. The value ranges of the result variables are set equal in all figures to ensure easy comparability. The scales were chosen in such a way that the relevant effects are properly visualized and actual extrema are partially not mapped. Furthermore, the wear profiles of the rolling elements and washers as well as the probability density function and the asperity contact pressure function are shown in Figure 7-9 for several timesteps. The calculated mass losses are then summarized in Table 3.
Due to the roller radius in and the logarithmic profile perpendicular to the direction of motion, the initial conditions featured typical characteristics for elliptical EHL contacts with a horseshoeshaped lubricant gap (Figure 5a-d), a total pressure distribution similar to Hertz (Figure 5e-h) as well as a constriction in the film height as well as a barely pronounced Petrusevich spike in the outlet region. Due to the smaller contact area, the pressures exceeded the Hertzian values determined for a perfect line contact (Table 1). In addition, the speed profile (Figure 1b) resulted in a slight asymmetry and tilting of the contact area. The clearly recognizable asperity contact pressure (Figure 5i-l) followed the lubricant height in shape and indicated an operation under mixed EHL conditions. The effects of the load case and oil on the lubricant gap and the pressure formation basically corresponded to theoretical expectations. Higher speed (load case 2) and higher viscosity (Figure 6) led to a larger and higher load (load case 3, 4) to a lower fluid film height. Accordingly, the asperity contact pressure behaved just oppositely. Moreover, the higher load led to a significantly larger elastically deformed contact area. It should be noted that although the contact width y was scaled in the same way for all cases, the contact length x was plotted normalized to the Hertzian width for reasons of reduced space requirements. Thus, the contact area actually expanded further in the x direction at higher loads as well.  The continuous progress of the simulated time led to material loss and changes of the surface topography of the washers and the rollers. In brief, the wear-depth distribution, as depicted in Figure  7, followed the product of asperity contact pressure and sliding velocity. The inner (negative y values) and outer halves (positive y values) of the raceway showed a slightly asymmetrical shape due to different velocities. Thereby, wear was higher in the area of the negative slip. This stands in good agreement with the experimental results from Terwey et al. [13,14]. Despite differences in absolute values due to the less severe conditions (lower loads and higher velocities) and the operation being in the mixed instead of the boundary lubrication regime, the qualitative profiles agree well. Furthermore, a more significant difference in the wear depth between the inner and outer halves of the raceway could be found on the washer. This was due to the fact that the calculated wear volume was distributed over a smaller circumference at the inner half of the raceway than at the outer area. Higher loads and lower sliding velocities resulted in a higher wear depth. This can be explained by lower lubricating film heights resulting in higher asperity contact pressures. For the same reason, the use of a lower viscosity oil (FVA1) compared to a higher viscosity oil (FVA 2) led to higher wear depths.
Proceeding wear also resulted in an adjustment of the probability density function of asperity heights, see Figure 8. By applying Sugimura's wear model, the probability density function of the surface heights was recalculated and adapted in each time step of the wear simulation assuming an exponentially distributed wear particle size. With a further continuing wear process, the PDF reached a stationary shape. This can be seen the most for the wear-intensive load case 3 with lower-viscosity oil FVA 1 (see Figure 8c). Here, the PDF of the surface heights after half of the simulation time only differed slightly from the PDF at the end of the wear simulation, indicating a run-in surface topography.  The PDF of the asperity heights were further applied to the calculation of the asperity contact pressures, which were a prerequisite for the application of the mixed friction model based upon the load sharing concept. As can be seen in Figure 9, the contact pressures were shifted to the bottom left with increasing simulation time, which was caused by the worn-out asperity peaks in terms of the modified probability density functions. Furthermore, in most cases the asperity contact pressure curves also reached a stationary state with slight differences to the final state being distinguishable for less harmful cases, see for example Figure 9f). Changes in the surface topography were also accompanied by transformations of the lubrication conditions, see Figure 10-13. In addition to the horseshoe-shaped constriction of the lubricant gap at the contact outlet region, there was a similarly pronounced minimum along the y = 0 axis. This was due to the low wear in this area, despite the high total and asperity contact pressure, since the conditions of the pure rolling (slide-to-roll ratio (SRR) = 0) prevailed. Furthermore, the asymmetry and slight tilting resulting from the speed and wear profiles could still be observed. Generally, the distributions of the lubricant gap and the pressure continued to correspond to the principles of the influences of speed, load and viscosity. In most cases, the fluid and pressure distribution only featured small differences between the half and end of the simulation time. While some load cases still showed significant asperity contact pressure after the simulated runtime, see Figure 12i,kFigure 12i; Figure   12k), others indicated a kind of completed running-in process with hardly any solid contact compared to the initial state, see for example Figure 13l).    Finally, the calculated wear masses for both bearing washers and all the rolling elements of the studied 81212 thrust roller bearing are summarized in Table 3. Again, it can be seen that the highest wear mass was calculated for load case 3 (highest load, lowest speed) in combination with the lower viscosity oil FVA 1. In contrast, load case 2 (lowest load, highest speed) in combination with the higher viscosity oil FVA 2 resulted in the lowest total wear.

Conclusions
Within the scope of this contribution, a numerical simulation method for studying the wear evolution of a thrust roller bearing under mixed elastohydrodynamic lubrication was introduced. Therefore, realistic geometries and operating conditions, especially a local velocity distribution, were considered. During the simulation time, the macroscopic lubrication conditions obtained by a FEMbased 3D EHL model were coupled with macroscopic and microscopic surface topography changes following a local Archard-type wear model, a Greenwood-Williamson-based load-sharing approach and the Sugimura surface adaption model. The simulations, carried out by means of four load cases and two different mineral oil types, indicated clear differences in the wear behavior, whereby stronger wear was favored by a lower velocity, higher load and lower oil viscosity. Thereby, asymmetrical wear profiles on the rollers and washers were revealed, whereas a negative slip on the outer half resulted in more a pronounced wear compared to the positive slip on the inner half of the raceway. This could be attributed to the velocity and slip distribution. Thus, the obtained results corresponded qualitatively well with the experimental results from the literature. In the future, it is intended to further enhance the EHL model by taking into account thermal and non-Newtonian effects. The use of locally modified probability density and asperity contact functions will also improve the prediction quality. However, the presented numerical simulation procedure already enables the determination of harmful load case and fluid property combinations. Moreover, the approach could be used to derive suitable running-in procedures.