the characteristics of the running-in process. The contact geometry includes the radii of curvature of the bodies, waviness and the surface roughness and structure. Kinematics and load can include, for example, the rolling motion, the normal and frictional force and the contact duration. Due to the running-in, a change in the coefficient of friction with the variation in speed can be determined in friction force tests. Depending on the minimum speed or lubricant film thickness, a more or less pronounced hysteresis results for the coefficient of friction. Knowledge of the running-in hysteresis can contribute to targeted optimization when designing tribological systems. Science and Research 24 Tribologie + Schmierungstechnik · volume 72 · issue 1/ 2025 DOI 10.24053/ TuS-2025-0004 1 Introduction, Motivation und Objective Highly stressed rolling-sliding contacts, such as the tooth flank contact, are tribological systems, for short tribosystem [CZIC10]. The properties of the tribological system are determined by the interactions between the contact partners, the lubricant and the ambient medium [POPO10]. For the tribosystem “gear”, the base and counter body can be defined as pinion and wheel, and the intermediate and ambient medium as gear oil and air. Depending on the initial state and the operating conditions, the tribological system changes continuously over time [PREX90, VOLG91, LÖPE15, LOHN16]. The change in tribological contact conditions is most pronounced during the initial loading of the tribosystem. This initial change is referred to as the running-in phase [ASMI92, CZIC10]. On the one hand, local deformation of the active surface takes place during the run-ning-in process (contact mechanics), which can be attributed to a local exceeding of the material strength or yield stress and can be quantified by recording the geometric surface properties. On the other hand, interactions of the triboelements can take place at the atomic and molecular level during runningin, which can be responsible for the formation of boundary layers as a result of adhesion, physisorption and chemisorption, for example [SCHU10, BREC16]. Figure 1-1 shows the variables influencing the running-in process in the rolling contact. The running-in process combines the initial state of the rolling surfaces (after production) with the state during operation. In addition to the kinematics and the lubrication, the contact geometry, the load and the material properties of the surface zone are among the central influences on Elastic-plastic micro contact calculation for large scale, lubricated rolling-sliding contacts Dieter Mevissen, Sebastian Sklenak, Christian Brecher, Thomas Bergs* Presented at GfT Conference 2024 The surface roughness and structure of rolling-sliding contacts have a decisive influence on the running behavior in terms of friction, wear and fatigue. Due to running-in effects, there is a continuous change of the surface roughness during the first load cycles. For this reason, the aim of this work was a calculation method to predict the friction force hysteresis during running-in. The basis of the method is a large-scale micro contact calculation based on the elastic halfspace theory in combination with an analytical stress calculation. The material behavior is modeled elastoplastically, whereby a limiting stress is used as a yield criterion instead of the limiting pressure usually used in the elastic half-space. The influence of the lubricant pressure is implemented by a lubricant boundary condition for the contact calculation in the elastic halfspace according to B OUSSINESQ , L OVE and H ARNETT . Finally, the calculation method is validated by means of friction force tests. For this purpose, S TRIBECK curves were recorded on a disc-on-disc test rig and compared with the calculation results. Keywords Micro contact calculation, half-space theory, Boussinesq, Love, Harnett, rolling-sliding contact, friction force measurement, lubricant boundary condition, macro and micro scale Abstract * Dr.-Ing. Dieter Mevissen Sebastian Sklenak M.Eng. Prof. Dr.-Ing. Christian Brecher Werkzeugmaschinenlabor der RWTH Aachen Campus-Boulevard 30, 52074 Aachen Prof. Dr.-Ing. Thomas Bergs Manufacturing Technology Institute MTI der RWTH Aachen Campus-Boulevard 30, 52074 Aachen This article presents a calculation method for predicting the friction force hysteresis during running-in (Figure1-2). The basis of the method is a large-scale microcontact calculation in the elastic half-space in combination with an analytic stress calculation. The material behavior is modelled elastic-plastic, whereby a limiting stress is used as a yield criterion instead of the limiting pressure usually used in the elastic half-space theory. The influence of the lubricant pressure is implemented by a lubricant boundary condition for the contact calculation in the elastic half-space according to B OUSSINESQ , L OVE and H ARTNETT . Finally, the calculation method is validated by means of friction force tests. For this purpose, S TRIBECK curves were recorded on a two-disk test rig. Each series of measurements started with a speed run-down in order to achieve continuous smoothing of the surface. This was followed by a speed run-up to determine the hysteresis of the friction coefficient with the smoothed surface. The experimental procedure was simulated on the basis of the measured surfaces and finally compared. 2 Elastic-plastic and stress-based micro-contact calculation The main challenge in designing the calculation method is the realization of a large-scale contact calculation with sufficient resolution to capture the surface roughness and structure. Figure 2-1 shows the schematic structure of the calculation method. The starting point of the calculation is the surface topography of the contact bodies Science and Research 25 Tribologie + Schmierungstechnik · volume 72 · issue 1/ 2025 DOI 10.24053/ TuS-2025-0004 Figure 1-1: Influence factors on running-in in tribological systems © WZL Figure 1-2: Objective and approach © WZL the plastic material behavior on the basis of a yield criterion. In the third step, the elastic stress tensor is converted into a plastic strain tensor using the plasticity model according to PRANDTL-REUSS [LEE12]. Finally, the deformation in the normal direction to the surface is calculated using the semianalytical approach according to JACQ (4 th step) [JACQ02]. A detailed description of the individual calculation steps is summarized in [MEVI21]. Key parameters and factors influencing the calculation result are discussed below. During the stress calculation in the micro-contact, interactions occur between the stress fields of neighboring roughness peaks, see Figure 2-2. To analyze the interactions, patterns with different distances are generated based on the pressure distribution for the reference model, Science and Research 26 Tribologie + Schmierungstechnik · volume 72 · issue 1/ 2025 DOI 10.24053/ TuS-2025-0004 in the initial state. The calculation method is based on an iterative solution algorithm with four main calculation steps. The iterative solution procedure is necessary due to the pronounced non-linearities in the elastic-plastic micro-contact calculation. The output variable of the calculation is the surface topography in the running-in state with a sufficient size for a standard-compliant roughness evaluation. In the first step, the pressure distribution for the purely elastic contact is calculated (Figure 2-1). The two-dimensional pressure distribution within the contact surface is then converted into a three-dimensional stress state on the basis of AIRY’s stress functions [AIRY62]. The stress tensor for each volume element forms the basis for calculating an equivalent stress and for modelling Figure 2-1: Structure of the Calculation Method © WZL Figure 2-2: Difference of pressure and volume-stress based plasticity modell © WZL in which the pressure distribution is duplicated. This creates a model contact. The maximum pressure of p H = 3139 MPa is the same for all calculations, so that the variants do not differ in the case of a pressure-based yield criterion. Figure 2-2 above shows the influence of the distance between two neighboring point contacts on the hydrostat and the MISES equivalent stress. The results are evaluated at the location of the maximum equivalent stress and normalized to the variant without direct neighbors (infinite distance). As the distance decreases, the hydrostat increases in magnitude and consequently the equivalent stress (deviator) for the contact point in the center decreases. The change is due to the superposition of stress components caused by radiation effects from the neighboring micro-contacts. In addition to the distance, the number of contact points has an influence on the hydrostat and the equivalent stress, see Figure 2-2 below. If the number of direct neighbors is further increased by a cross pattern, the hydrostat increases in magnitude. The increase converges to a limiting value, which is physically to be expected because the influence of highly distant contact points continuously decreases with increasing distance. The maximum reduction of the MISES equivalent stress is approx. 16 % for the examples presented here. In addition to the interaction between individual roughness contacts, knowledge of the permissible evaluation range in the depth direction is a key requirement. In micro-contact calculation, the minimum evaluation depth as a function of the meshing resolution on the surface plays a particularly important role. Figure 2-3 on the left shows the stress component σ x as a function of depth for three different crosslinking resolutions (dy, dz) on the surface (reference model, dx = 0.1 µm). While at greater depths (x > 0.01 mm) the stress curves are identical for the different crosslinking resolutions, a difference is recognizable at shallower depths, with the stress component σx converging towards negative infinity in all cases for x = 0 mm. According to the analytical calculation, this course is not to be expected, as the magnitude of the stress component σx forms its maximum directly at the surface at x = 0 with a magnitude equal to the contact pressure according to H ERTZ p H . It can also be seen that the stress curves reach the amount of H ERTZ pressure at different depths. The higher the meshing resolution, the lower the depth at which the H ERTZ stress is first reached. To quantify the minimum evaluation depth, a convergence analysis is carried out as a function of the surface crosslinking in Figure 2-3 on the right. The minimum evaluation depth is the numerical limit up to which the semi-analytical stress calculation can be evaluated depending on the meshing resolution. The criterion for determining the minimum evaluation depth is the depth at which the normal stress σ x (normal to the surface) assumes the value of the H ERTZ stress. In the case of square meshing, the minimum evaluation depth decreases almost linearly with the meshing resolution at the surface, see Figure 2-3 on the right. For rectangular surface elements, the minimum evaluation depth is determined by the longer edge length, so that an optimization of the calculation time for the two-dimensional pressure calculation with rectangular elements is at the expense of the minimum evaluation depth. To define the resolution dy, dz in the finely meshed contact area by, bz, a convergence study is carried out for the stress and deformation calculation. For the convergence analysis, a contact calculation is also carried out for a section of a rough surface. In the stress calculation, the individual microstress fields are resolved with suffi- Science and Research 27 Tribologie + Schmierungstechnik · volume 72 · issue 1/ 2025 DOI 10.24053/ TuS-2025-0004 Figure 2-3: Minimum evaluation depth depending on the surface mesh size © WZL re 3-1 shows the solution of the elastic half-space taking into account a lubricant boundary condition. With the lubricant boundary condition, the pressure in the areas without solid contact is not set to zero, but to the value of the lubricant pressure. As the lubricant pressure causes deformation cross influences on the surrounding elements, the right-hand side of the linear system of equations changes for all elements. It is therefore no longer possible to simply delete the rows and columns, as is done with the adhesion boundary condition due to setting the contact pressure to zero. Therefore, the linear system of equations is solved in the first step and then sorted by size. A negative contact pressure exists if the element has to be contracted for contact determination. This is precisely the case for the roughness valleys, which lie behind the roughness peaks and do not come into contact despite the elastic flattening. Finally, the lubricant boundary condition is applied in this step, in which a negative contact pressure is set equal to the lubricant pressure. In the third step, the linear system of equations is separated into a variable and a constant term at the transition from positive to negative contact pressure, see Figure 3-1. The constant term no longer changes in the further elastic contact calculation because the average lubricant pressure in the roughness valleys is specified as constant. In the last step, the linear system of equations is reduced and the constant term on the left-hand side is subtracted from the vector on the right-hand side. It is precisely at this point that the deformation cross influences caused by the lubricant pressure are taken into account by the subtraction when solving the linear system of equations. The solution steps are repeated until no more negative contact pressures occur after the first iteration step. Science and Research 28 Tribologie + Schmierungstechnik · volume 72 · issue 1/ 2025 DOI 10.24053/ TuS-2025-0004 cient accuracy from a resolution of dy = dz = 0.8 µm, see Figure 2-4 left. To compare different meshing resolutions, the stress depth curves for a contact point in the microcontact are also evaluated. It can be seen that the stress curves converge, whereby deviations can be seen for a resolution of dy = dz = 2.0 µm, particularly at greater depths compared to a finer resolution. Using the same procedure, the local changes in the contact geometry for different meshing resolutions are evaluated in the deformation calculation, see Figure 2-4 on the right. A positive deformation value means that the surface has risen at this point. The amounts of deformation on the surface converge with increasing meshing resolution. If the mesh is too coarse (dy = dz = 2.0 µm), the deformation tends to be calculated too small. From an element size of dy = dz = 0.8 µm, the calculation results for the example contact are comparable. Based on the convergence analysis for the stress and deformation calculation, the high-resolution contact area (by, bz) must therefore be meshed with an element size of at least dy = dz = 0.8 µm. The meshing resolution at the surface results in a minimum evaluation depth of x min = 1.3 µm. 3 Lubricant pressure within elastic half-space theory The elastic half-space offers an efficient way for solving the contact problem of two bodies with arbitrary geometry and is suitable for predicting the geometric surface change in fluid-free contact. In this section, a model for the lubricant pressure in the elastic half-space is presented and a boundary condition for the lubricant volume equilibrium is explained according to [MEVI19]. Figu- Figure 2-4: Convergence Study for Plastic Deformation © WZL In addition to the lubricant boundary condition, the calculation of mixed friction states requires the coupling of the micro-contact calculation based on the elastic half-space theory with the knowledge from the contact calculations of the EHD theory with ideally smooth surfaces. In contrast to the established procedure in the literature, the contact calculation in the elastic half-space is used as the main model, which draws on the results of the EHD contact calculation (auxiliary model) [BART10, SAUE18]. The basic idea of the calculation approach is that a certain volume of lubricant is conveyed into the lubrication gap by the hydrodynamics of the lubrication gap and is distributed evenly in the contact geometry deformed under load. Based on this basic idea, a further boundary condition can be defined for the micro-contact calculation in the elastic half-space during iterative solution finding, which takes into account the balance between the hydrodynamically conveyed lubricant volume and the pocket volume of the roughness valleys under load, see Figure 3-2 center. The boundary condition of the lubricant volume equilibrium is the link between the calculation results of the EHD theory at macro level and the micro-contact calculation in the elastic half-space. In iterative solution finding, the fulfilment of this new equilibrium condition is made possible by the degree of freedom gained in the specification of the lubricant pressure in the area of the roughness valleys. Science and Research 29 Tribologie + Schmierungstechnik · volume 72 · issue 1/ 2025 DOI 10.24053/ TuS-2025-0004 Figure 3-1: Solving algorithm for Lubricant Boundary Condition for elastic half-space © WZL Figure 3-2: Determination of Lubricant Pressure © WZL intersects the calculated curve beforehand, so that a potential for reducing the Ra value remains with a further reduction in speed. The comparison with the experiment is based on the measured and calculated roughness values at the end of the speed run-down, see Figure 4-1 bottom right. For the increasing normal force, the calculation correctly maps the increasing smoothing. In addition, the degressive course of the smoothing with increasing normal force is captured by the calculation. Overall, the change in the Ra value for all three normal forces is overestimated by the calculation, leaving room for optimizing the elastic-plastic material behaviour and modelling the material distribution. The basis for the calculation of the friction coefficient is the solid contact ratio, which is a central output variable from the large-scale micro-contact calculation. Figure 4-2 shows the change in the solid contact ratio for the different normal forces plotted against the average lubrication pocket height. According to the experiment, an elastic-plastic microcontact calculation (speed downrun) is carried out in the first step. This is followed by a purely elastic contact calculation (speed run-up). The reversal point between speed run-down and run-up depends on the lubricant volume conveyed at minimum speed or the minimum average lubricant film thickness h EHD,min . As the normal force increases, the reversal point shifts to lower solid contact ratios. The maximum possible smoothing is therefore not achieved for the higher normal forces for the operating parameters, see also Figure 4-1. Each time the operating conditions change, the surface roughness is reduced again, resulting in hysteresis for all Science and Research 30 Tribologie + Schmierungstechnik · volume 72 · issue 1/ 2025 DOI 10.24053/ TuS-2025-0004 4 Prediction of running-in hysteresis in disc-on-disc contact During the determination of S TRIBECK curves by a speed run-down and speed run-up, a hysteresis occurs for the friction coefficient. In the following section, the calculation method is applied to the tests with different normal forces and a prediction of the hysteresis is attempted. The calculation parameters and the material model correspond to the documentation in [MEVI21] in Chapter 7. Figure 4-1 shows the calculated change in the Ra value for the three different normal forces. Due to the elasticplastic contact calculation and the consideration of the lubricant pressure, the Ra value is continuously reduced as the average lubrication pocket height decreases (speed deceleration). The leftmost point of each curve marks the state with a solid contact ratio of ψ load = 100 % and thus the maximum possible smoothing for the respective normal force. Due to the higher total deformation with increasing normal force, the maximum possible smoothing is highest for the largest normal force. The vertical lines in the diagram in Figure 4-1 represent the link between experiment and calculation. The minimum mean lubricant film thickness h EHD,min is equivalent to the minimum lubricant volume conveyed, which occurred in the experiment at a speed of n = 200 rpm. The mean lubricant film thickness was calculated according to the approximation formula of H AMROCK / D OWSON for the line contact and is dependent on the normal force, see Figure 4-1. For the low normal force, the maximum possible smoothing was achieved as a result of the calculation result with a speed of n = 200 1/ min. For the two higher normal forces, the amount of smoothing increases, but the vertical line (minimum mean lubricant film thickness h EHD,min in the experiment) Figure 4-1: Prediction of surface roughness and change of surface roughness © WZL three normal forces. The hysteresis is most pronounced for the initial load at the lowest normal force and hardly recognizable at the highest normal force. In addition, the curve for a purely elastic contact calculation with the initial geometry is shown for comparison at the low normal force (Figure 4-2 left). The elastic-plastic contact calculation reduces the average lubrication pocket height, particularly with high solid contact ratios, which can be attributed to the smoothing of the roughness. To validate the contact calculation, the coefficients of friction from the experiment are finally compared with the results of the calculation, see Figure 4-2. The coefficient of friction is calculated on the basis of the solid contact ratio in the same way as in [MEVI19] with the same values for the solid and liquid coefficients of friction. The speed in the calculation results from the lubricant volume equilibrium and is equivalent to the mean lubrication pocket height. The elastic-plastic micro-contact calculation maps the renewed hysteresis between the speed run-down and run-up with each increase in the normal force. For the low normal force of FN = 1086 N, the hysteresis is overestimated by the calculation, resulting in a steeper drop in the friction coefficient. One reason for the deviation is the differences between the calculated and measured surface topography after the speed run, which are present after the deformation calculation. The differences in the surface topography can influence Science and Research 31 Tribologie + Schmierungstechnik · volume 72 · issue 1/ 2025 DOI 10.24053/ TuS-2025-0004 Figure 4-3: Friction Force Coefficients in Calculation and Experiment © WZL Figure 4-2: Calculation Results of Method © WZL [CZIC10] Czichos, H.; Habig, K.-H.: Tribologie-Handbuch. Tribometrie, Tribomaterialien, Tribotechnik. 3. Aufl. Wiesbaden: Vieweg+Teubner, 2010 [JACQ02] Jacq, C.; Nélias, D.; Lormand, G.; Girodin, D.: Development of a Three- Dimensional Semi-Analytical Elastic-Plastic Contact Code. In: J. Tribol., 124. Jg., 2002, Nr. 4, S. 653-667 [LEE12] Lee, Y.; Barkey, M.; Kang, H.: Metal Fatigue Analysis Handbook. Practical Problem-Solving Techniques for Computer-Aided Engineering. 1. Aufl. Waltham, MA: Butterworth-Heinemann, 2012 [LOHN16] Lohner, T.: Berechnung von TEHD Kontakten und Einlaufverhalten von Verzahnungen. Diss. TU München, 2016 [LÖPE15] Löpenhaus, C.: Untersuchung und Berechnung der Wälzfestigkeit im Scheiben- und Zahnflankenkontakt. Diss. RWTH Aachen University, 2015 [MEVI19] Mevissen, D.; Löpenhaus, C.; Bergs, T.: Calculation of mixed friction conditions in large-scale rolling-sliding contacts for different surface structures. In: Forsch Ingenieurwes, 83. Jg., 2019, Nr. 3, S. 351-366 [MEVI21] Mevissen, D.: Vorhersage der geometrischen Oberflächenveränderung im Wälzkontakt. Diss. RWTH Aachen University, 2021 [POPO10] Popov, V.: Kontaktmechanik und Reibung. Von der Nanotribologie bis zur Erdbebendynamik. Berlin: Springer, 2010 [PREX90] Prexler, F.: Einfluss der Wälzflächenrauheit auf die Grübchenbildung vergüteter Scheiben im EHD-Kontakt. Diss. TU München, 1990 [SAUE18] Sauer, B.; Magyar, B.; Oehler, M.: Schneckengetriebewirkungsgrade. Abschlussbericht zum Forschungsvorhaben FVA-Nr. 729 I, Heft 1226, Forschungsvereinigung Antriebstechnik e.V., Frankfurt a.M., 2018 [SCHU10] Schulz, J.; Holweger, W.: Wechselwirkung von Additiven mit Metalloberflächen. Renningen: Expert, 2010 [VOLG91] Volger, J.: Ermüdung der oberflächennahen Bauteilschicht unter Wälzbeanspruchung. Diss. RWTH Aachen University, 1991 Science and Research 32 Tribologie + Schmierungstechnik · volume 72 · issue 1/ 2025 DOI 10.24053/ TuS-2025-0004 the course of the solid contact ratio and thus the friction coefficient. The lower hysteresis with increasing normal force is correctly captured by the calculation. In addition, the friction coefficient curves intersect with decreasing speed, so that the calculation reflects the change from an increasing friction coefficient to a decreasing friction coefficient with increasing normal force. In summary, it can be stated that the running-in hysteresis of the friction coefficient can be predicted by the calculation method. In addition to the selection of the meshing parameters, the use of a stress-based yield criterion and the consideration of the lubricant pressure during the deformation of the half-space are key elements for the successful application. Acknowledgement The authors would like to thank the the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) [407625150 and 390969378] for providing the financial means to carry out the research project on which the results presented here are based. Literatur [AIRY62] Airy, G.: On the Strains in the Interior of Beams. In: Phil. Trans. R. Soc. London, 153. Jg., 1862, S. 49-79 [ASMI92] ASM International Handbook Committee: ASM Handbook. Friction, Lubrication, and Wear Technology. Bd. Nr. 18, Materials Park, OH: American Society for Metals, 1992 [BART10] Bartel, D.: Simulation von Tribosystemen. Grundlagen und Anwendungen. Wiesbaden: Vieweg+ Teubner, 2010 [BREC16] Brecher, C.; Löpenhaus, C.; Greschert, R.: Interaktion von fertigungs- und betriebsbedingten Grenzschichten imWälzkontakt. In: Tagungsband zur 57. Arbeitstagung “Zahnrad- und Getriebeuntersuchungen”. Aachen, 11./ 12.05.2016. Eigendruck des WZL-Getriebekreises, 2016