Geometrical Optimization of the EHL Roller Face/Rib Contact for Energy Efficiency in Tapered Roller Bearings

In the context of targeted improvements in energy efficiency, secondary rolling bearing contacts are gaining relevance. As such, the elastohydrodynamically lubricated (EHL) roller face/rib contact of tapered roller bearings significantly affects power losses. Consequently, this contribution aimed at numerical optimization of the pairing’s macro-geometric parameters. The latter were sampled by a statistical design of experiments (DoE) and the tribological behavior was predicted by means of EHL contact simulations. For each of the geometric pairings considered, a database was generated. Key target variables such as pressure, lubricant gap and friction were approximated by a meta-model of optimal prognosis (MOP) and optimization was carried out using an evolutionary algorithm (EA). It was shown that the tribological behavior was mainly determined by the basic geometric pairing and the radii while eccentricity was of subordinate role. Furthermore, there was a trade-off between high load carrying capacity and low frictional losses. Thereby, spherical or toroidal geometries on the roller end face featuring a large radius paired with a tapered rib geometry were found to be advantageous in terms of low friction. For larger lubricant film heights and load carrying capacity, spherical or toroidal roller on toroidal rib geometries with medium radii were favorable.


Introduction
Tapered roller bearings (TRBs) are being widely used for high load supporting applications in mechanical, automotive and transmission engineering where their main characteristics-high radial and uniaxial load carrying capacity, demountability and adjustable clearance-can be fully exploited. They are therefore frequently found, for example, in rail or passenger vehicle wheel hub assemblies, gas turbine engines or worm and bevel gear stages or differentials [1]. Some investigations were therefore concerned with maximizing the load carrying capacity, fatigue, stiffness, load distribution and dynamics. In particular, analytical and numerical approaches for the design of geometrical bearing features have been applied [2][3][4][5]. Thereby, machine learning and data mining methods as well as optimization algorithms were also employed [6][7][8]. Despite secondary contacts or the micro-geometry were analyzed in few studies [9][10][11], most efforts focused on the macro-geometry of primary contacts between rolling elements and raceways, for example by profiling/crowning [12].
In the context of diminishing resources, rising environmental awareness and stricter legal requirements, secondary contacts are gaining relevance. For example, the roller face/rib-contact significantly affects the friction and power losses, especially when the TRB is subjected to high axial loads. Here, the lubrication conditions have a decisive role, which is why some investigations also dealt with its numerical modeling. Basically, the roller end face/rib contact differs from classical EHL point-contacts in terms of geometry and kinematics, whereby the effects due to the geometric pairing as well as spinning friction in particular have to be considered [13,14] Zhang et al. [15] studied the end face/rib-contact using a finite difference (FD) method with Gauss-Seidel iterations to solve the Reynolds equation simultaneously with the equations for elasticity based upon the half space theory. It was shown that the hydrodynamic fluid film formation together with the superimposed elastic deformation have both to be taken into account (elastohydrodynamic lubrication, EHL) and that the hydrodynamic lubrication (HL) theory alone is not sufficient [16]. It was also demonstrated that the curvatures of the roller and rib significantly influence the contact conditions. The favorable ratios of the end face radius to the rib radius in the range 0.6 to 0.8 were later confirmed by Wang et al. [17] with simulation studies extended to include the effects induced by surface roughness using the flow factor model from Patir and Cheng [18] and the asperity contact model following Greenwood and Tripp [19]. Colin et al. [20] investigated the starved EHL rib/roller end contact and indicated effects of contact geometry, load, rotational speed and oil supply conditions on the performance in terms of fluid film height and traction. Similarly, Fujiwara et al. [21] found a large ratio of roller end to rib face radius of slightly below 0.85 to be favorable in terms of skewing and fluid film thickness.
In summary, some research efforts were made to identify roller face and rib radii ratios with low frictional losses and high load carrying capacity. However, this was limited to comparatively simple geometric pairings as well as few parametric studies. In principle, the design space and modern manufacturing possibilities may offer much more options, which can be assessed with the aid of contact simulation tools as well as numerical machine learning and multi-objective optimization algorithms [22]. Therefore, this contribution is aimed at introducing a numerical optimization process for macro-geometric parameters in the EHL roller face/rib-contact featuring complex kinematics and geometries.

Materials and Methods
The general optimization procedure utilized within the scope of this contribution and as illustrated in Figure 1 has already been successfully applied to derive tailored micro-texture geometries for EHL contacts by Marian et al. [23]. Based on the geometric parameters describing the roller end face and rib geometries, they were sampled by a statistical design of experiments (DoE) and the tribological behavior was predicted by means of EHL contact simulations. Special attention has been paid to modeling the complex geometry pairings and kinematics, which exceeds most cases of classic EHL point-or linecontacts. For each of the geometric pairings considered, a database was generated on the basis of which key target variables such as coefficient of friction (COF), maximum pressure and minimal lubricant film height were approximated by a meta-model and an optimization was carried out. Finally, the predicted optima were verified again by EHL simulations. This whole process was done separately for each geometric pairing to finally compare their respective optima and derive generalized design recommendations. The relevant aspects are described in detail in the following Sections 2.2-2.6. At first, however, the parameters and properties of the TRB, load case and lubricant considered in this study are defined in Section 2.1.

TRB Load Case, Kinematics and Lubrication
A standard TRB 30207, as illustrated in Figure 2a, with an inner ring rotational speed n I of 1000 min −1 under a pure axial load of 12 kN, to get an equal load at the rib for each roller, was studied. For the materials, typical values for rolling bearing steel (100Cr6) were considered with a Young's modulus E of 210,000 N/mm 2 and a Poisson's ratio ν of 0.3. ratios with low frictional losses and high load carrying capacity. However, this was limited to comparatively simple geometric pairings as well as few parametric studies. In principle, the design space and modern manufacturing possibilities may offer much more options, which can be assessed with the aid of contact simulation tools as well as numerical machine learning and multi-objective optimization algorithms [22]. Therefore, this contribution is aimed at introducing a numerical optimization process for macro-geometric parameters in the EHL roller face/rib-contact featuring complex kinematics and geometries.

Materials and Methods
The general optimization procedure utilized within the scope of this contribution and as illustrated in Figure 1 has already been successfully applied to derive tailored microtexture geometries for EHL contacts by Marian et al. [23]. Based on the geometric parameters describing the roller end face and rib geometries, they were sampled by a statistical design of experiments (DoE) and the tribological behavior was predicted by means of EHL contact simulations. Special attention has been paid to modeling the complex geometry pairings and kinematics, which exceeds most cases of classic EHL point-or line-contacts.  For each of the geometric pairings considered, a database was generated on the basis of which key target variables such as coefficient of friction (COF), maximum pressure and minimal lubricant film height were approximated by a meta-model and an optimization was carried out. Finally, the predicted optima were verified again by EHL simulations. This whole process was done separately for each geometric pairing to finally compare their respective optima and derive generalized design recommendations. The relevant aspects are described in detail in the following Sections 2.2-2.6. At first, however, the parameters and properties of the TRB, load case and lubricant considered in this study are defined in Section 2.1.

TRB Load Case, Kinematics and Lubrication
A standard TRB 30207, as illustrated in Figure 2a, with an inner ring rotational speed nI of 1000 min −1 under a pure axial load of 12 kN, to get an equal load at the rib for each roller, was studied. For the materials, typical values for rolling bearing steel (100Cr6) were considered with a Young's modulus E of 210,000 N/mm 2 and a Poisson's ratio ν of 0.3. The kinematics for the roller face/rib contact of the TRB (see Figure 2b) were described by a discretized velocity matrix and mapped to the surface matrix of the EHL simulation (see Section 2.4) as depicted in Figure 2c. Thereby, the velocities in x-and ydirection of both contacting bodies were calculated at each grid node. When assuming the outer ring to be stationary, the rotational speed of each roller element can be written as follows [24]: with the pitch diameter given by: the average diameter of the roller dW and the rib angle β. Here, the grid had lateral dimensions (Xstart, Xend, Ystart, Yend) to cover eight times the Hertzian width (see Section 2.4) and was meshed with a number of 1024 nodes each in the x-(Nx) and y-direction (Ny), which also corresponded to the EHL calculation domain. The discrete distance between the nodes was given by: Each node was expressed by a column indexed by i in the x-direction and a row indexed by m in the y-direction. Thus, the distance da of the contact point in the middle of the grid to the rotation axis of the rolling element or the ring was derived: The kinematics for the roller face/rib contact of the TRB (see Figure 2b) were described by a discretized velocity matrix and mapped to the surface matrix of the EHL simulation (see Section 2.4) as depicted in Figure 2c. Thereby, the velocities in xand y-direction of both contacting bodies were calculated at each grid node. When assuming the outer ring to be stationary, the rotational speed of each roller element can be written as follows [24]: with the pitch diameter given by: the average diameter of the roller d W and the rib angle β. Here, the grid had lateral dimensions (X start , X end , Y start , Y end ) to cover eight times the Hertzian width (see Section 2.4) and was meshed with a number of 1024 nodes each in the x-(N x ) and y-direction (N y ), which also corresponded to the EHL calculation domain. The discrete distance between the nodes was given by: Each node was expressed by a column indexed by i in the x-direction and a row indexed by m in the y-direction. Thus, the distance d a of the contact point in the middle of the grid to the rotation axis of the rolling element or the ring was derived: With the angle γ resulting from the number of columns and the distance of notes in x-direction, the distance d a to the rotation axis for each node was obtained: The velocity v tangential to the axis of rotation for each node was thus determined by [24]: v(i, m) = 2·d a ·π·n I/W , whereas the components of this vector in the xand y-direction were calculated by: v x = v· cos γ and v y = v· sin γ.
Thus, the velocity distribution of both, the inner ring and the rolling element, in each node could be determined. It should be noted that both have different algebraic signs and the former had larger values than the latter.
The contact was assumed to be lubricated with a sufficient amount (fully flooded) of FVA 3 reference mineral oil. The respective properties are summarized in Table 1. Pressure p dependent density ρ and viscosity η were considered following the models from Dowson and Higginson [25]: as well as Roelands [26]: respectively. Non-Newtonian rheology, i.e., shear-thinning of the fluid, was modeled by means of the Ree-Eyring approach [27]:

Roller Face/Rib Contact Geometry
While the geometry of the roller end face can be described by the face radius r R and the eccentricity e, the rib is characterized by the raceway/rib angle β as well as the rib radius r B , see Figure 2b.
The position of the contact point is another constraint. Therefore, the contact point lies in the middle of the rib height. The variations of the radii and the eccentricity yield the rib angle. Within this contribution, the principle pairing possibilities of sphere, cone and torus were evaluated by means of parameter variations and limitations as shown in Table 2. Based on these parameters, the geometric pairing was transferred to the EHL simulation (Section 2.4). Before, this required further processing, for which parts of the so-called PIMP-method as introduced by Wirsching et al. [28] was applied. At first, the relative positions of the two contacting bodies and their local coordinate system within the global one were determined. The contact zone position was calculated using the position vectors q (x, y, z) and the vectors n (x, y, z) normal to the bodies' surface in the local coordinate systems. The latter are in the same position and oriented in opposite directions. Subsequently, the geometry of each body was described analytically in the Cartesian coordinate system G. For the sphere, this was defined as: for the cone as: and for the torus as: The vectors q and n determined the position and orientation of the point of contact and a projection plane for each body. A ray tracing method projected the geometry on the projection planes in the Cartesian coordinate system G: The two projection planes were parallel and congruent, and the contact point was located in the center. By converting above equations, a function defining the distance between the two rigid bodies to the projection plane was derived. For the case of the sphere, this can be written as follows: where t = 0 is the point of contact and both geometries of the two bodies were described by the distance functions g 1 (x, y) and g 2 (x, y). With: a substitute geometry was finally generated to equivalently describe the contact of the two rigid bodies [28]. The corresponding expressions for the cone could be written as: − q x n x + q y n y − m 2 q z n z + q x n x + q y n y − m 2 q z n z Lubricants 2021, 9, 67 6 of 19 and the torus: x + q 2 y + q 2 z · q x n x + q y n y + q z n z + R 2 · −q x n x − q y n y + q z n z − r 2 · q x n x + q y n y + q z n z · t respectively.

Design of Experiments
In order to generate a sufficient database for the deduction of correlations between input and target variables with a minimum of computational effort, methods of statistical DoE can be applied. Here, a so-called latin hypercube sampling (LHS) from the group of uniformly distributed test field was used, since this is particularly suitable for generating a wide range of approximation or meta-models (see Section 2.5). Thereby, the trial points are partitioned in such a way that they fill the design space as uniformly as possible and provide information about almost every area with little computational effort. The LHS elements are created by subtracting a random number between zero and one from each element of a matrix with columns consisting of random permutations of the numbers {1, . . . , n f }, with n f being the number of input factors, and then dividing this value by the number of trial points n p . For this, the Statistics and Machine Learning toolbox of MATLAB (The Mathworks Inc., Natick, MA, USA) was utilized to obtain LHS designs optimized regarding the distances of the trial points by means of the Maximin criterion [29,30]. The number of trial points was chosen as 20 for n f = 1, 50 for n f = 2 and 100 for n f = 3, see Table 2. The concave geometry pairings had the additional constraint, that the absolute value of the rib radius r B had to be larger than the roller face radius r R . Therefore, the number of trial points was doubled and the points outside the constraint were deleted.

EHL Contact Simulation
Numerical EHL modeling is usually done by a coupled solution of the lubricant's hydrodynamics and the elastic deformation of the contacting bodies. For this purpose, the software tool TELOS 5.0 (Schaeffler Technologies AG & Co. KG, Herzogenaurach, Germany) was used for a quasi-static, isothermal, lubricated contact simulation with a non-Newtonian rheology model. For a better conditioning and stabilized solution of the system of equations, which are briefly introduced below, the coordinates x and y, the lubricant gap h, the hydrodynamic pressure p as well as viscosity η and density ρ were normalized with the Hertzian contact widths a', the Hertzian approach of both bodies c', the maximal Hertzian pressure p 0 as well as the ambient viscosity η 0 and density ρ 0 : Due to the geometries differing strongly from classical point contacts, the normalization was done using the values of a ball-on-plane pairing with a ball diameter of 100 mm.
The geometry was described by the lubricant gap equation taking into account the distance of the rigid bodies H 00 , the approximated undeformed geometry of the respective pairing G(X, Y) as well as the combined elastic deformation: The calculation of the elastic contact was based on a variational principle [31]. Here, the pressure distribution and the actual contact zone minimized the absolute complementary distortion energy. An iterative algorithm following Polonsky and Keer [32] was used to solve the contact problem described here. This was based on the conjugate gradient method and the fast Fourier transformation, allowing a stable and efficient solution. The quasi-stationary Reynolds differential equation was applied to account for the fluid flow and hydrodynamic pressure generation: Here, the sum velocities u 1,2 and v 1,2 in xand y-direction mapped the motion of the contacting bodies (see Section 2.1). Non-Newtonian effects were considered according to the rheology model of Ree-Eyring: Furthermore, the following boundary conditions were assumed: here dΩ denoted the edges of the domain (23) P(X, Y) = 0; ∀(X, Y) ∈ Ω, Ω denoted all points in the comupational space (24) The equilibrium between normal load F and the hydrodynamic pressure was ensured by the load balance equation: The whole system of equations was solved iteratively based upon the finite difference method as well as the elastic half-space theory while applying multigrid methods and fast Fourier transformation as well as relaxation for numerical damping and stabilization [33]. The discretization of the domain was controlled stepwise with a binary system from 64 to 2048 points and the size of the calculation domain was chosen to be rectangular with lateral dimensions of eight times the Hertzian contact width a'. The numerical solution scheme is shown in Figure 3 and for more numerical details, the interested reader is referred to [34][35][36]. After the converged calculation, the COF was calculated: Lubricants 2021, 9, x FOR PEER REVIEW Furthermore, the following boundary conditions were assumed: The equilibrium between normal load F and the hydrodynamic pressure w sured by the load balance equation: = ∬ ( , ) .
The whole system of equations was solved iteratively based upon the finite dif method as well as the elastic half-space theory while applying multigrid methods a Fourier transformation as well as relaxation for numerical damping and stabilizatio The discretization of the domain was controlled stepwise with a binary system fro 2048 points and the size of the calculation domain was chosen to be rectangular w eral dimensions of eight times the Hertzian contact width a'. The numerical s scheme is shown in Figure 3 and for more numerical details, the interested reade ferred to [34][35][36]. After the converged calculation, the COF was calculated:

Machine Learning
Based on the simulation results of the DoE introduced in Section 2.3, the tribo behavior can be described by approximation-or meta-models [38,39]. Thus, the e

Machine Learning
Based on the simulation results of the DoE introduced in Section 2.3, the tribological behavior can be described by approximation-or meta-models [38,39]. Thus, the effect of the different factor combinations in the form of geometric parameters on the target variables (pressure, fluid film height, COF) can be determined and an optimization can be performed in a computationally efficient way. The meta-model, a representation of the simulation results, may be based on different approaches. One of the most elementary forms is polynomial regression, where an unknown function is approximated by a polynomial solution of usually low degree [29,40]. An extension with local character based on local weighting functions is the method of moving least squares (MLS) [41]. Usually, linear or quadratic terms are utilized for the base function. Furthermore, the Kriging model, also known as Gaussian process regression, is an accurate surrogate model due to its flexibility in approximating different and complex response functions by interpolating the data points and providing a confidence interval of the prediction [42]. The assessment of the prediction quality of the respective approach represents a central aspect of meta-modeling. An automated optimization was enabled by a so-called Meta-model of Optimal Prognosis (MOP). The concept developed by Most and Will [40,43] builds upon a two-step elimination process of non-significant variables based on a coefficient of determination (COD) and importance (COI) as well as the automatic selection of the most suitable of aforementioned approximations based on a coefficient of prognosis (COP). The latter is advantageous due to its automatic scaling to values between 0 and 1. For example, a COP of 0.9 corresponds to a prediction quality of 90% for new data points. Here, the initial data set was split into training and test data with a ratio of 70/30 in such a way that the response ranges in the data sets displayed maximum conformity to the overall data set. The MOP was implemented in the software OptiSlang version 8.1.0 (Dynardo GmbH, Weimar, Germany).

Optimization
For the determination of global optima of the geometric parameters based upon the MOP, various algorithms can be utilized. Basically, there are gradient-based methods [44], response surface methods (RSM) [45] and methods inspired by nature, such as particle swarm optimization (PSO) [46] or evolutionary/genetic algorithms (EA/GA) [47]. The latter were used within the scope of this study since they are an efficient and flexible way to optimize unknown problems and were also implemented in the software OptiSlang. Thereby, a population of solution candidates for the optimization problem was simulated in a selection process with underlying evolutionary principles. At first, the initial individuals were evaluated with respect to the optimization goals, which defined subsequent mating selection and the generation of child individuals from each parent. The new generation was created by recombining the characteristics of the parent individuals as well as a mutation and integrated into the existing population under environmental selection, replacing parent individuals. After each cycle, the termination conditions in form of a maximum number of iterations and a quality stagnation for a defined number of generations were checked. The EA setting underlying this contribution is summarized in Table 3.

Results
In the following, sample results of the EHL simulations are shown in Section 3.1, followed by the results of meta-modeling in Section 3.3, optimization in Section 3.3 as well as verification in Section 3.4.

Pressure, Film Thickness and Friction
The calculated pressure and lubricant gap distributions of one representative parameter set near the center of the parameter field for each pairing (see Table 4) are illustrated in Figure 4a-l). The main direction of motion of the rolling element and inner ring and thus the main fluid flow was in the x-direction. Accordingly, the pressure in the contact inlet increased and reached a maximum in the contact center before dropping again at the contact outlet. No pronounced Petrusevich spike was observed. Correspondingly, there was an elastic flattening in the contact center as well as a minimum in the lubricant film height with a slight horseshoe-shape near the contact outlet. In the y-direction, the distributions were nearly axisymmetric with small form deviations due to the variation in flow velocities.   For all geometric pairings, a different degree of ellipticity of the elastically deformed contact area with the longer semi-axis in the x-direction was found. This was less distinct for the sphere/cone and torus/cone pairings, but strongly pronounced for the sphere/torus and torus/torus pairings. Moreover, there were considerable differences in the magnitudes of the pressure maxima and the lubricant film heights between the different geometry pairings, see Table 4. In particular, the sphere/cone and torus/cone pairings exhibited larger pressures due to smaller elastically deformed contact areas while the sphere/torus and torus/torus geometries resulted in higher lubricant gaps. Accordingly, differences in the COF were observed. It was noticeable that the sphere/cone and torus/cone pairings featured rather low COFs and the sphere/torus and torus/torus geometries exhibited higher levels.

Influence of the Roller Face/Rib Geometry
In the following, the effects of the various geometric parameters on friction as well as the lubricant gap as an indicator of potential wear and the maximum hydrodynamic pressure will be shown. The results from the EHL simulations from the LHS are provided in the Supplementary Material and the calculated response surfaces of the MOP are displayed in Figure 5a-r).
Generally, red color indicates higher and blue color lower values of studied variables. The maximum hydrodynamic pressure and the COF were most accurately approximated by anisotropic kriging except for the pairing sphere/cone, where MLS was used for the pressure and linear regression for the COF. MLS was also applied for minimum lubricant gap of the sphere/cone pairing. The lubricant gap of all other pairings was represented by linear regression. The corresponding COPs are listed in Table 5. Table 5. COPs of the MOP of each geometry pairing shown in Figure 5. For the sphere/cone geometry (Figure 5a-c), the maximum hydrodynamic pressure, the minimum lubricant gap and the COF only depended on the face radius r R as single geometrical variable. Apparently, a smaller radius resulted in higher pressure and COFs as well as smaller lubricant gaps and vice versa. The same behavior was basically also observed for the torus/cone (Figure 5d-f), sphere/torus concave (m-o) and torus/torus concave (p-r) pairings. Furthermore, the eccentricity e of roller radius had a rather small or no influence on calculated results (Figure 5d-f,j-l,p-r) since the response surfaces barely changed in the eccentricity direction. Moreover, for all geometries, lower maximum pressures and larger lubricant gaps were achieved for larger radii r B and r R . However, the lowest COFs were partly obtained within the studied parameter field, see Figure 5i,l,o,r). Lubricants 2021, 9, x FOR PEER REVIEW 11 of 18

Optimization of the Roller Face/Rib Geometry
Multi-objective optimization by means of the EA was performed based upon the MOP response surfaces in order to minimize the COFs or maximize the minimum lubricant gap. This led to a Pareto front of which the extrema in the respect to the optimization goals (lowest friction or highest film height) are summarized in Figure 6. The resulting maximum pressure, minimum film height and COF are illustrated in Figure 6a-c). It can be seen that both objectives led to contradicting results since the lowest friction (sphere/cone and torus/cone geometry in Figure 6c) also resulted in rather small lubricant gaps (sphere/cone and torus/cone geometry in Figure 6b) while higher lubricant gaps also led to higher COFs. The hydrodynamic pressure behaved similar to the minimal lubricant gap and thinner lubricant gaps resulted in higher hydrodynamic pressures. Generally, toroidal rib geometries generated lower hydrodynamic pressures. Accordingly, the optimized geometry parameters in Table 6 differed in dependency of the objective. When targeting a higher minimum film height, the optima were found near the extrema of the studied parameter field (largest possible radii) while the optima were found rather in the middle of the rib radius and the roller face radius field for the pairings torus/torus, sphere/torus concave and torus/torus concave when minimizing friction.
or no influence on calculated results (Figures 5d-f,j-l,p-r) since the response surfaces barely changed in the eccentricity direction. Moreover, for all geometries, lower maximum pressures and larger lubricant gaps were achieved for larger radii rB and rR. However, the lowest COFs were partly obtained within the studied parameter field, see Figures  5i,l,o,r).

Optimization of the Roller Face/Rib Geometry
Multi-objective optimization by means of the EA was performed based upon the MOP response surfaces in order to minimize the COFs or maximize the minimum lubricant gap. This led to a Pareto front of which the extrema in the respect to the optimization goals (lowest friction or highest film height) are summarized in Figure 6. The resulting maximum pressure, minimum film height and COF are illustrated in Figures 6a-c). It can be seen that both objectives led to contradicting results since the lowest friction (sphere/cone and torus/cone geometry in Figure 6c) also resulted in rather small lubricant gaps (sphere/cone and torus/cone geometry in Figure 6b) while higher lubricant gaps also led to higher COFs. The hydrodynamic pressure behaved similar to the minimal lubricant gap and thinner lubricant gaps resulted in higher hydrodynamic pressures. Generally, toroidal rib geometries generated lower hydrodynamic pressures. Accordingly, the optimized geometry parameters in Table 6 differed in dependency of the objective. When targeting a higher minimum film height, the optima were found near the extrema of the studied parameter field (largest possible radii) while the optima were found rather in the middle of the rib radius and the roller face radius field for the pairings torus/torus, sphere/torus concave and torus/torus concave when minimizing friction.

Verification of Optimized Roller Face/Rib Geometry
Finally, the optimization results based upon the MOP/EA were verified. Therefore, the parameter combinations from the optimization for each geometry pairing (Table 6) were recalculated by means of EHL contact simulations. The maximum hydrodynamic pressure, the minimum lubricant gap and the COF were compared to the prediction from the MOP/EA, see Table 7. Thereby, rather small deviations between the prediction and the EHL simulation for most variables and geometry combinations occurred. The minimum lubricant film height and COF in particular were excellently predicted with errors below 2%. Only the accuracy of the hydrodynamic pressure for the sphere/torus, torus/torus concave after minimizing the COF as well as for sphere/torus and torus/torus when maximizing the lubricant gap showed larger discrepancies.

Load Carrying Capracity Versus Friction-Influences of the Geometries
The roller face/rib contact was characterized by rather moderate loads and high sliding speeds. Hence, pressure and lubricant gap distribution ( Figure 4) showed typical characteristics for EHL contacts operating under such conditions [48] and similarities to other studies on that contact reported in literature [16,17,21]. As a result of the velocity differences between the two bodies, there were also considerable frictional losses depending on the respective geometry pairings ( Figure 5). Generally, larger radii tended to result in wider contact areas and thus also in higher fluid film heights as well as lower pressures required to carry the load. Correspondingly, larger radii exhibited a bigger load-bearing capacity, i.e., better protection against potential solid asperity contact and wear. However, these did not necessarily lead to lower COFs due to more fluid being sheared within the wider contact areas. In this respect, the findings corresponded well to the experimental studies from Korrenn [11]. Going beyond other studies [11,16,17,21], it was shown that, in addition to the radii, the basic geometric pairings also had a major influence on the tribological properties in terms of load-bearing capacity and frictional losses. This can be attributed to differences in the shapes and sizes of the contact areas as well as the velocity distributions, which is also substantiated by the experimental investigations from Jamison et al. [10]. Nevertheless, there remains a certain trade-off between high load carrying capacity and low friction (see Figure 6b,c). Thereby, spherical or toroidal geometries on the roller end face featuring a large radius paired with a tapered (cone) rib geometry without curvature were found to be advantageous in terms of low friction losses. For larger lubricant film heights to separate surface asperities or to support higher loads, however, spherical or toroidal roller on toroidal rib geometries with radii in the center of the analyzed parameter range were favorable (Table 6). Similar to [16,17], favorable ratios of the radii were withing the range of 0.8-0.9 for concave geometries. However, advantageous ratios were substantially larger for convex geometry pairings (~1-2) or not available after all due to one contacting partner having an uncurved surface, especially when optimizing for minimum friction.

Applicability and Limitations
Within the scope of this contribution, one TRB type, lubricant and load case with pure axial load were investigated. The conclusions are mainly valid in this respect. Changes in the boundary conditions may lead to an alteration or shift of presented results. For example, radial loads lead to the formation of a load zone and thus to variable forces and speeds for each individual roller end face/rib EHL contact during one revolution of the bearing, which would demand for time-transient calculation. In addition, the influences of the load carrying capacity of the secondary roller end face/rib contact on the behavior of the primary roller/raceway contacts were not further considered. For additional quantitative prediction accuracy regarding the frictional behavior, which would allow a comparison with corresponding experimental data, a more detailed characterization of the rheological fluid behavior as well as an adaptation of the used models might be essential [49]. Furthermore, as for most numerical EHL studies, the underlying simulations for this investigation were subject to certain assumptions (regarding the applicability of the Reynolds equation, fully flooded conditions, linear elastic material behavior etc.), which, however, can be considered to be relatively well satisfied for the studied contact. Surface roughness and thermal effects might have an influence on the tribological behavior. Taking them into account, however, also increases computational costs, which is why they were neglected considering the large number of calculations performed within the DoE. Therefore, we consider the implications to be conclusive and to have universalizability.

Prediction Capability of the ML Approach
For most geometries and objectives, the adopted ML approach based upon the MOP provided a very good prognosis, which was reflected in the high COP values mostly above 90% (Table 5) as well as the small deviations of the optimization based upon MOP/EA compared to EHL verification calculations mainly below 2% (Table 7). This indicated the possession of a sufficient database as well as the quality of the meta-models automatically selected by the MOP. Only torus rib geometries yielded somewhat larger errors, especially for very small rib radii, which can be attributed to instabilities in the EHL simulations (nonconverged solutions) and thus in the underlying database in that region. However, these areas were not pertinent for the optimization regarding friction or load carrying capacity, which is why the good prediction quality in most other areas prevailed. The approach to predict and optimize the tribological behavior with the chosen approach consisting of MOP and EA can thus be evaluated as suitable. This does not have to be limited to the macro-geometry of the roller end face/rib EHL contact but can also be applied to its micro-geometry, e.g., for tailoring roughness, grinding patterns or surface textures, or to other EHL contacts.

Conclusions
The tribological behavior of the roller end/face rib contact of tapered roller bearings as a function of different macro-geometries was investigated within the scope of this contribution. Therefore, geometric parameters describing the roller end face and rib geometries were sampled by a statistical design of experiments. The tribological behavior was predicted by means of EHL contact simulations taking into account the complex geometry pairings and kinematics. For each of the geometric pairings considered (sphere/cone, torus/cone, sphere/torus, torus/torus, sphere/torus concave and torus/torus concave), a database was generated. Key target variables such as COF, maximum pressure and minimal lubricant film height were approximated by a meta-model of optimal prognosis and optimization was carried out using a genetic algorithm. Finally, the predicted optima were verified again by EHL simulations. For the studied load case, contact conditions and geometries, the following conclusions could be drawn:

•
The introduced machine learning approach featured excellent prediction quality and was able to efficiently and effectively support geometry optimization of EHL contacts with respect to the tribological behavior.

•
The tribological behavior of the roller end face/rib contact is mainly determined by the basic geometric pairing and the radii; eccentricity is of subordinate role. • There is a certain trade-off between high load carrying capacity and low frictional losses.

•
If the bearing is subjected to rather low axial loads and/or higher velocities, i.e., there is moderate risk of mixed lubrication and wear, but energy losses are to be minimized, spherical or toroidal geometries on the roller end face featuring a large radius paired with a tapered (cone) rib geometry without curvature are favorable.
• If the bearing expects higher demands regarding axial load carrying capacity and/or lower speeds, spherical or toroidal roller on toroidal rib geometries with radii in the center of the analyzed parameter range are advantageous. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.