Energy ‐ Based Modelling of Adhesive Wear in the Mixed Lubrication Regime

: Adhesive wear in dry contacts is often described using the Archard or Fleischer model. Both provide equations for the determination of a wear volume, taking the load, the sliding path and a set of material parameters into account. While the Fleischer model is based on energetic approaches, the Archard formulation uses an empirical factor—the wear coefficient—describing the intensity of wear. Today, a numerical determination of the wear coefficient is already possible and both approaches can be deduced to a local accumulation of friction energy. The aim of this work is to enhance existing energy ‐ based wear models into the mixed lubrication regime. Therefore, the pressure distribution within the contact area will be determined numerically taking real surface topographies into account. The emerging contact area will be divided into one part of solid and a second part of elastohydrodynamically lubricated (EHL) contacts. Based on the resulting pressure and shear stress distribution, the formation of micro cracks within the loaded volume will be described. Determining a critical number of load cycles for each asperity, a locally resolved wear coefficient will be derived and the local wear depth calculated.


Introduction
An increased application of low-viscosity lubricants in machine elements increasingly leads to operating conditions in the mixed lubrication regime along with the risk for the occurrence of wear. Since between highly loaded asperities small contact areas form, the material's yield stress can be reached quickly and micro welds may form [1]. The subsequent separation of surfaces can occur in planes different from the original ones and material can be transferred from one body to the other. Figure 1 illustrates this process, referred to as adhesive wear, for a pair of contacting spherical shaped asperities-before (a), during (b) and after (c) the contact. The asperities are assumed to be subjected to a (local) normal load and a relative velocity . Neglecting chemical effects, the occurring wear volume can be described inter alia by probabilistic models, e.g., the Archard theory [3]. Archard as well assumes spherically shaped asperities with a radius of , such that at position a circular contact area of π between a pair of contacting asperities forms. He furthermore assumes this area to be so small that the occurring local pressure immediately reaches the yield stress of the surface material. According to [4], the yield stress can be equated with the surface hardness , which leads to the relation / .
Since each asperity in contact is loaded by the yield stress, this relation holds for all contact spots. The volumes of evolving wear particles can each be represented by the volume of a detached asperity, which equals a hemisphere: 2/3 π . Defining additionally a micro sliding path in the length of the asperity diameter, 2 , a local relative wear volume can be defined as: By the accumulation of all local relative wear volumes for the entire contact area and considering the force equilibrium for all asperities, ∑ (where is the total normal force working on the considered system), the total relative wear volume can be given as: The variable is the total sliding path. The factor takes into account that an instant formation of wear does not occur for every contact. It thus constitutes a parameter for the probability of the formation of wear or-in other words-for the intensity of wear. Finally, using the simplification /3, the total wear volume can be defined as: Equation (3) is known as the Archard equation and the factor as the Archard wear coefficient. The wear coefficient is usually determined empirically-for metals it extends in general over a range from 10 to 10 [5].
Another model describing the adhesive wear mechanism was proposed by Fleischer [6]. The underlying theory is based on an energy perspective and takes the performed frictional work into account: The factor is the friction coefficient. The intensity of wear is described by the apparent frictional energy density * , which leads to the Fleischer wear equation: 1 * Comparing Equations (3) and (5) it can trivially be derived, that the apparent frictional energy density can be given as a function of the Archard wear coefficient and vice versa: This work focusses on the numerical determination of the Archard wear coefficient and the numerical calculation of adhesive wear in local resolution (as it was already presented in [2]). The pressure distribution, taking rough surfaces into account and being one of the main drivers for the occurrence of wear, will be calculated for both dry and mixed lubricated contacts. In the latter case, the force transmission from one body to the other through the asperity contacts is supported by an evolving pressurized elastohydrodynamic fluid film.

Numerical Contact Model
The contact model is based on the elastic half-space theory [7] and the Hertz-Signorini-Moreau condition. The latter describes the relationship between an arbitrary gap between two surfaces within the considered area and the corresponding pressure . For dry contacts, the condition takes the following form [5]: This condition states firstly, that an overlap of volumes cannot occur ( 0) and secondly, that in case of a mechanical contact ( 0) a pressure will evolve ( 0). The relationship between a local pressure at position (respectively a local force ) and the resulting elastic deformation at position k can be expressed using the influence matrix , which provides the relationship between both variables [4]: (8) The number of asperities is . The number of elements contained in the influence matrix can quickly become very large ( ), especially for a finely resolved mesh. Due to the need for a high resolution (if the surface roughness shall be taken into account) this issue is tackled by means of the Method of Combined Solutions [8].
For mixed lubricated contacts, Equation (7) has to be adapted by the condition that occurring gaps will be filled with pressurized lubricant (underlying an elastohydrodynamic pressure ). The pressures within the gaps support the asperity contacts in transmitting forces form one body to the other, while at the same time separating the surfaces partwise from each other. Consequently, the number of contacting asperities decreases as well as the mechanical pressures for each contact do. According to the elastohydrodynamically lubricated (EHL) theory, a central film thickness ℎ for point contacts can be given in the following form [9]: is a geometry parameter, a velocity parameter, a lubrication parameter and a load parameter. describes the ellipticity of the contact area. For rough surfaces, the film thickness is treated as the distance between the median planes of the deformed asperities of respective contacting surfaces. The case of mixed lubrication is defined by a fluid film, which is not sufficiently thick to separate the surfaces completely from each other. Thus, parts of the surface will remain in contact but, however, being supported by the elastohydrodynamic pressure. For contacting asperities, Equation (7) keeps being valid, while for non-contacting parts (the gaps), the formation of an elastohydrodynamic pressure is postulated. The lubricant pressure is assumed to follow the shape of the Hertzian pressure distribution for smooth surfaces [10,11] and can be given for any specific point within the contact area as: The resulting mean value ̅ of all gaps ( 0) within the contact area equals the film thickness and thus defines the distance between a pair of contacting bodies with rough surfaces: Satisfying Equations (7), (10) and (11), a force equilibrium between the external normal force and the sum of all solid ( ) and EHL ( ) forces prevails: The number of contacting asperities is and of non-contacting asperities (gaps) . The force equilibrium and the film thickness ℎ are depicted in Figure 2.

Numerical Wear Model
In order to develop a numerical wear model, the asperities are modelled as cubes, having an edge length of , Figure 3 (in accordance with a rectangular discretized half-space). As described before, the wear coefficient can be understood as the probability for the detachment of a single wear particle. It thus forms the reciprocal value of a number of sustainable load cycles and can be written as [3,12]: For a numerical calculation of , the theory of Continuum Damage Mechanics is used, describing the progressive weakening of solid materials under a cyclic shear stress [13][14][15][16]. In Figure 4 a control volume, loaded by a cyclic shear stress is shown. The damaged area consists of several micro cracks, depicted as black areas (d ), while the control area d is represented in light gray.
It is assumed that any load above a specific endurance limit will cause micro cracks within the loaded volume. The ratio of the damaged area ∑ d , defined by the normal vector , to a cross section d can be expressed as [17]: The factor is a probability factor. For isotropic materials, the damage ratio is independent of the normal direction , which leads to . The damage propagation increases exponentially and can be expressed as a function of the performed load cycles [17]. If the endurance limit is not reached, no damage propagation will occur: The parameter describes the amount of damage propagation and can be given as follows [18]: The parameter is the hardening coefficient, the hardening exponent and the true failure stress [18]. If after cycles a critical value is reached, the number of load cycles performed equals the critical number: . The validity of the model was proved by Beheshti et al. in [18] by a comparison of predicted wear coefficients and experimental results as a function of the friction coefficient. In Figure 5, the comparison for aluminum 6061-T6 is depicted, while the experiments where conducted on a pin-ondisk test. The tests were conducted with a pin (diameter 8 mm) made of aluminum in contact with a disk made of stainless steel (SS 304, diameter 100 mm). The sliding velocity was set within a range of 30-120 mm/s, the normal load in a range from 3-30 N. The test duration was 2-5 h. The amount of wear was obtained by weighing the pin before and after the test.
(a) (b) Figure 5. Pin-on-disk tribometer (a) and comparison of experimental and simulation results (b) for the prediction of the wear coefficient. [18] Applying the aforementioned model to all contacting asperities, the wear coefficient can be determined in local resolution. Therefore it is assumed, that the asperity is loaded by a cyclic shear stress , which is the product of the local pressure and the friction coefficient: (17) Thus, the local wear coefficient is instantly a function of the local shear stress. Assuming an incremental wear path in the length of (element edge length) and an incremental wear volume of (volume of a cubic element), the local relative wear volume now takes the form: Here, the probability (or wear-) coefficient is considered in local resolution. Rearranging the equation and using the relationship , the local wear depth can be expressed as . Taking additionally the number of load cycles into account, the cumulative wear path can be considered. The locally resolved wear depth can finally be given in the form of a vector:

Calculation Algortihm
The simulation follows the algorithm depicted in Figure 6. The left side shows the calculation of the pressure distribution, which follows an iterative approach. The first step is to read the material, geometry and load parameters. Based on this, the Hertzian pressure distribution is calculated and the discretized mesh is defined. For the determination of the pressure distribution the influence matrix is calculated. Depending on the fluid and the EHL parameters, a film thickness ℎ is calculated. Before the calculation starts, an initial guess of the solid body overlap is estimated and imprinted on the system. Taking the micro geometry and the solid body overlap into account, the gap distribution is calculated based on the elasto-plastic half space. Within the Hertzian contact area a mean value ̅ of the gap distribution is determined. If the deviation between ̅ and the film thickness ℎ exceeds a specific limit , the solid overlap is updated (increased in case of ̅ /ℎ 1 and decreased in case of ̅ /ℎ 1) and the calculation will be repeated. If the deviation is small enough, the mean gap is assumed to be equal the EHL film thickness and the iteration is finished. Finally, the calculation of the solid pressure distribution follows, the local wear coefficient is calculated and the wear distribution is determined (right side of Figure 6).

Model System
Although the presented contact and wear models have been developed to describe wear in machine elements like bearings or gears, for investigations into the wear mechanism it is more suitable to abstract the geometry to a simplified model system. Thus, characteristic contact geometries and loading situations of rolling element bearings have been transferred to a model contact consisting of two discs. The material is bearing steel 100Cr6 (AISI52100). The discs are pressed together by a normal force of 6850 N. Each disc has a radius of 60 mm and is driven independently from the other with different circumferential velocities ( and ). Thus, in the contact area an entrainment velocity of 1/2 5 m/s and a relative velocity of 0.5 m/s evolve, depicted in Figure 7. the surface geometries of both discs are projected to one resulting (combined) model surface, which contacts to an opposing flat elasto-plastic half-space surface. In Table 1 the material, geometry and simulation parameters are listed.

Simulation of Gap and Pressure Distributions
The results of the simulation of the gap distribution for a dry and a mixed lubricated contact are depicted in the following Figures 8 and 9. In Figure 8 the side view (XZ view) of the deformed surface geometry between the combined surface (black) and the half-space surface (zero line) for a dry contact is shown. The gap height (in -direction) is given in µm on the ordinate, while the abscissa shows the -direction along the long half-axis. The latter is normalized to the Hertzian contact length . On the top right corner, the corresponding undeformed surface is depicted keeping the aspect ratio constant. The numerically calculated film thickness is ℎ 587 nm (for mixed lubrication), which is marked by a blue dotted horizontal line. Because of numerical deviations the calculated film thickness differs slightly from the analytical solution of 552 nm.  Figure 9 shows the corresponding results for mixed lubrication. It can be seen, that due to the supporting effects of the EHL pressure the deformation and the number of contacting asperities has significantly decreased. In the following Figures 10-12 the simulation results for the pressure distribution are depicted. In each figure, the results for a dry contact are shown on the left side, results for mixed lubrication on the right. In both cases, only the occurring pressures for contacting asperities are drawn, meaning that the EHL pressures are blanked out. Figure 10 shows the top (XY) view of the contact zone. The abscissa again shows the -direction along the long half-axis ( ) of the contact zone, while the ordinate shows the -direction along the short half-axis ( ). Both axes are normalized to the respective Hertzian contact length. It can be seen, that the dry contact fills out almost the entire Hertzian contact area (denoted in red dotted line), while the mixed contact area does not. Furthermore, a cyclic waviness can be identified, which corresponds to a waviness of the surface along the -direction due to the manufacturing process. , an EHL part does not exist (see Equation (12) Thus, the percentage of the force, transmitted through asperity contacts is much lower, than in the case of a dry contact: In Figure 11 the pressure distribution is shown in the XZ view. The abscissa now shows the long half-axis ( ) and the ordinate shows the normalized pressure / . The red dotted line depicts the corresponding Hertzian pressure distribution for smooth surfaces. In both cases (dry and mixed), a multitude of asperities reaches the yield stress, why the maximum achievable normalized pressure is limited to 3.8. The pressure density is much higher in the case of the dry contact, as it also can be observed in Figure 10. Finally, Figure 12 shows the pressure distribution in YZ view. The abscissa now shows the short half-axis ( ).

Simulation of Wear Distribution
The simulation of the adhesive wear follows the model presented in chapters 2.2. and 2.3. The shear stress and thus the friction coefficient and local pressure have large impacts on adhesive wear. While the local pressure can be calculated numerically as shown in the previous chapter, the determination of the (solid) friction coefficient is still difficult. To illustrate the influence of both parameters on the wear coefficient, a parameter study was carried out. Varying the (local) pressure in a range of 3500 to 8000 MPa and the friction coefficient in a range of 0.10 to 0.20 leads to wear coefficients of up to 1.3 • 10 (while the minimum value for the wear coefficient was set to 1.0 • 10 , according to [5]). In Figure 13 the corresponding curves are depicted. The factor ̅ shall here be referred to as the wear intensity factor. The ratio between both, the mixed lubricated and the dry factor, is 8.77% and thus of a similar magnitude to the percentage of force, transmitted through asperity contacts (8.80%). Figure 14 depicts the simulation results for the wear depth (according to Equation (19)) for the dry contact. It shows a section through the XZ plane, while the total wear volume is distributed evenly over the entire circumference of the discs. The wear depth (depicted in black color) for both contacting bodies is projected to the half-space, depicted in gray color. The abscissa is showing the -direction along the long half-axis ( ) of the contact zone, while the ordinate shows the wear depth. It can be seen, that the occurring wear spreads over the entire contact zone, although areas of high work density (i.e., high product of relative velocity and pressure) show a significantly higher density of wear (here: in the central area, where the amount of pressure is the largest). The maximum peaks reach values of  In Figure 15 the corresponding wear for the mixed lubricated contact is shown. It is trivial, that a reduction of solid forces due to the formation of a supporting fluid film leads to a reduction of wear. As already shown, the wear intensity factor is significantly lower in case of the mixed regime. The occurring wear decreases in general, while the maximum wear in particular decreases to a value of , , 10.6 µm.

Conclusion and Outlook
Within the presented work, simulation results for a dry and a mixed lubricated contact have been shown and compared. The wear depth has been determined by calculating a pressure distribution, taking a rough surface geometry into account, and subsequently calculating a locally resolved and energy-based wear coefficient. It has been shown that the cumulative force, transmitted by contacting asperities can be reduced by an evolving elastohydrodynamic film, followed by a significant reduction of wear. The results show a high accumulation of wear for specific spots of single sharp asperities, while it is assumed not to correspond to real systems. It should be noted, that it is possible to avoid this numerical fault by updating the contact geometry after each load cycle, such that the pressure distribution will adapt the new geometry and the high pressure regime will relocate to less-worn areas. The disadvantage of this approach is the strongly increasing computing time. The calculation of the pressure distribution takes by far the majority of the entire calculation time. In addition, the creation of wear-resistant chemical layers, built up by chemical additives, may be taken into account in future works.
Besides the pressure, the coefficient of friction is the main driver for the amount of the calculated wear coefficients. In this work, the coefficient of friction is assumed to be 0.12, which involves huge uncertainty. Further work is therefore focused on the numerical and experimental determination of the friction coefficient.
Furthermore, an experimental validation of the wear model on real machine elements will be the subject of future works. As already shown in [2], wear in machine elements shows several effects, which are not considered in this paper. Nevertheless, experiments were already carried out on thrust roller bearings in [2] and show a basic agreement of experimental and simulation results. The experimental part was conducted using an FE8 test rig [20], Figure 16. In this test rig, two test bearings (2) are mounted on a common shaft, while the housing washers are fixed in the housing. The shaft is radially supported by two cylindrical roller bearings and driven by an electric motor. A static axial load is applied through a cup spring assembly (1) such that both test bearings are subjected to the same load. The test lubricant was a mineral oil (viscosity 186 cSt at 40 °C). The investigations were carried out at low relative speeds between the contacting bodies, so that the film thickness ratio was also low ( 0.07). Three load cases with axial loads of 25, 37.5 and 50 kN were defined, resulting in HERTZian pressures of 1028, 1216 and 1383 MPa, respectively. For each load case (LC) a test with a constant shaft speed of 7.5 min was conducted. Each test duration was 80 h (36000 shaft revolutions). The resulting experimental wear volume was determined by tactile measurement of the profile before and after each test. The corresponding simulations were conducted without consideration of an EHL film and assuming a friction coefficient of 0.18.
In Figure 17 the results of the experiments are shown on the left side (a) and for the simulations in the middle (b). The left half of each picture is the inner part of the washer/roller contact. In this region negative slip occurs, while on the outer part (right half) positive slip occurs. In load case 1, little to no wear can be observed. In load case 2, a wear depth of around 5 μm in the area of the negative slip occurred, while in the area of positive slip considerably less wear can be observed. In load case 3, the wear depth is around 10 µm. The simulation results show that the amount of wear can be calculated in the right order of magnitude, but the correct prediction of the local amounts is still difficult. Two of the reasons may be the missing consideration of additional wear mechanisms, like abrasive wear, and the missing consideration of the influence of positive/negative slip. Funding: The publication of this article was funded by the Open Access Fund of the Leibniz University Hannover, Germany.

Conflicts of Interest:
The authors declare no conflict of interest.