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

Abstract

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.

1. 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 $f^i$ and a relative velocity $v_{\mathrm{rel}}$.

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 $r$, such that at position $i$ a circular contact area of $A_{\mathrm{R}}^i\approx \mathsf{\pi }{r}^2$ between a pair of contacting asperities forms. He furthermore assumes this area to be so small that the occurring local pressure $p^i$ immediately reaches the yield stress $p_{\lim }$ of the surface material. According to [4], the yield stress can be equated with the surface hardness $H$, which H leads to the relation $H=f^i/A_{\mathrm{R}}^i$. 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: $V_{\mathrm{w}}^i=2/3\mathsf{\pi }{r}^3$. Defining additionally a micro sliding path in the length of the asperity diameter, $S_{\mathrm{w}}^i=2r$, a local relative wear volume can be defined as:

(1)$\frac{V_{\mathrm{w}}^i}{S_{\mathrm{w}}^i}=\frac{1}{3}\mathsf{\pi }{r}^2=\frac{1}{3}{A}_{\mathrm{R}}^i=\frac{1}{3}\frac{f^i}{H}$

By the accumulation of all local relative wear volumes for the entire contact area and considering the force equilibrium for all $Z$ asperities, ${{\displaystyle \sum }}_{i=1}^{\mathrm{Z}}{f}^i=F_{\mathrm{N}}$ (where $F_{\mathrm{N}}$ is the total normal force working on the considered system), the total relative wear volume can be given as:

(2)$\frac{V_{\mathrm{w}}}{S_{\mathrm{w}}}={{\displaystyle \sum }}_{i=1}^{\mathrm{Z}}\frac{V_{\mathrm{w}}^i}{S_{\mathrm{w}}^i}=\frac{1}{3}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{Z}}\frac{f^i}{H}=\frac{1}{3}{K}_1\frac{F_{\mathrm{N}}}{H}$

The variable $S_{\mathrm{w}}$ is the total sliding path. The factor $K_1$ 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 $K=K_1/3$, the total wear volume $V_{\mathrm{w}}$ can be defined as:

(3)$V_{\mathrm{w}}=K\frac{F_{\mathrm{N}}{S}_{\mathrm{w}}}{H}$

Equation (3) is known as the Archard equation and the factor $K$ as the Archard wear coefficient. The wear coefficient is usually determined empirically—for metals it extends in general over a range from $K={10}^{-8}$ to ${10}^{-1}$ [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 $W_{\mathrm{R}}$ into account:

(4)$W_{\mathrm{R}}=\mu F_{\mathrm{N}}{S}_{\mathrm{w}}$

The factor $\mu$ is the friction coefficient. The intensity of wear is described by the apparent frictional energy density $e_{\mathrm{R}}^{*}$, which leads to the Fleischer wear equation:

(5)$V_{\mathrm{w}}=\frac{1}{e_{\mathrm{R}}^{*}}{W}_{\mathrm{R}}$

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:

(6)$e_{\mathrm{R}}^{*}=\frac{1}{K}\mu H \leftrightarrow K=\frac{1}{e_{\mathrm{R}}^{*}}\mu H$

This work focusses on the numerical determination of the Archard wear coefficient $K$ 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.

2. Calculation Method

2.1. 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 $g$ between two surfaces within the considered area and the corresponding pressure $p$. For dry contacts, the condition takes the following form [5]:

(7)$g\ge 0, p\ge 0, p·g=0$

This condition states firstly, that an overlap of volumes cannot occur ($g\ge 0$) and secondly, that in case of a mechanical contact ($g=0$) a pressure will evolve ($p\ge 0$). The relationship between a local pressure $p^i$ at position $i$ (respectively a local force $f^i$) and the resulting elastic deformation $w^k$ at position k can be expressed using the influence matrix $C=\left[C_{ik}\right]$, which provides the relationship between both variables [4]:

(8)$w^k={{\displaystyle \sum }}_{k=1}^Z{C}_{ik}{f}^i$

The number of asperities is $Z$. The number of elements contained in the influence matrix can quickly become very large ($Z^2$), 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 $p_{\mathrm{EHL}}^j$). 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 $h_{\mathrm{c}}$ for point contacts can be given in the following form [9]:

(9)$h_{\mathrm{c}}=2.69R{U}^{0.67}{G}^{0.53}{W}_{}^{-0.067}\left(1-0.61{\mathrm{e}}^{-0.73\kappa }\right)$

$R$ is a geometry parameter, $U$ a velocity parameter, $G$ a lubrication parameter and $W$ a load parameter. $\kappa$ 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 $j$ within the contact area as:

(10)$g^j>0 \to p^j\approx p_{\mathrm{Hertz}}^j$

The resulting mean value $\overline{g}$ of all gaps ($g^j>0$) within the contact area equals the film thickness and thus defines the distance between a pair of contacting bodies with rough surfaces:

(11)$\overline{g}=h_{\mathrm{c}}$

Satisfying Equations (7), (10) and (11), a force equilibrium between the external normal force $F_{\mathrm{N}}$ and the sum of all solid ($f_{\mathrm{sol}}^i$) and EHL ($f_{\mathrm{EHL}}^j$) forces prevails:

(12)$F_{\mathrm{N}}=F_{\mathrm{sol}}+F_{\mathrm{EHL}}={\displaystyle \sum }_{i=1}^{n_{\mathrm{sol}}}{f}_{\mathrm{sol}}^i+{\displaystyle \sum }_{j=1}^{n_{\mathrm{EHL}}}{f}_{\mathrm{EHL}}^j$

The number of contacting asperities is $n_{\mathrm{sol}}$ and of non-contacting asperities (gaps) $n_{\mathrm{EHL}}$. The force equilibrium and the film thickness $h_{\mathrm{c}}$ are depicted in Figure 2.

2.2. Numerical Wear Model

In order to develop a numerical wear model, the asperities are modelled as cubes, having an edge length of $d$, Figure 3 (in accordance with a rectangular discretized half-space).

As described before, the wear coefficient $K$ can be understood as the probability for the detachment of a single wear particle. It thus forms the reciprocal value of a number $N_{\mathrm{c}}$ of sustainable load cycles and can be written as [3,12]:

(13)$K=\frac{1}{N_{\mathrm{c}}}$

For a numerical calculation of $N_{\mathrm{c}}$, the theory of Continuum Damage Mechanics is used, describing the progressive weakening of solid materials under a cyclic shear stress $\tau$ [13,14,15,16]. In Figure 4 a control volume, loaded by a cyclic shear stress $\tau$ is shown. The damaged area consists of several micro cracks, depicted as black areas ($\mathrm{d}{A}_{D_m}$), while the control area $\mathrm{d}A$ is represented in light gray.

It is assumed that any load above a specific endurance limit ${\tau }_{\mathrm{e}}$ will cause micro cracks within the loaded volume. The ratio $D$ of the damaged area ${\displaystyle \sum }\mathrm{d}{A}_{D_n}$, defined by the normal vector $\mathit{m}$, to a cross section $\mathrm{d}A$ can be expressed as [17]:

(14)$D\left(\mathit{m}\right)=k\frac{{\displaystyle \sum }\mathrm{d}{A}_{D_m}}{\mathrm{d}A}\in \left[0,1\right]$

The factor $k$ is a probability factor. For isotropic materials, the damage ratio is independent of the normal direction $\mathit{m}$, which leads to $D\left(\mathit{m}\right)=D$. The damage propagation increases exponentially and can be expressed as a function of the performed load cycles $n$ [17]. If the endurance limit ${\tau }_{\mathrm{e}}$ is not reached, no damage propagation will occur:

(15)$D\left(n\right)=D_n=\{\begin{array}{c}1-\left(1-D_{n-1}\right)B_n ,\mathrm{if}\tau >{\tau }_{\mathrm{e}}\\ D_{n-1},\mathrm{else}\end{array}$

The parameter $B_n$ describes the amount of damage propagation and can be given as follows [18]:

(16)$B_n=\frac{{\left(\frac{{\tau }_{\mathrm{e}}}{L}\right)}_{}^{1+M}+\left(\frac{3{\tau }_{\mathrm{f}}}{4L}\right){\left(\frac{{\tau }_{\mathrm{e}}}{L}\right)}_{}^{1+M}}{\frac{\tau }{L\left(1-D_{n-1}\right)} + \left(\frac{3{\tau }_{\mathrm{f}}}{4L}\right){\left(\frac{{\tau }_{\mathrm{e}}}{L}\right)}_{}^{1+M}}$

The parameter $L$ is the hardening coefficient, $M$ the hardening exponent and ${\tau }_{\mathrm{f}}$ the true failure stress [18]. If after $n$ cycles a critical value $D_n=D_{\mathrm{c}}$ is reached, the number of load cycles performed equals the critical number: $n=N_{\mathrm{c}}$.

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-on-disk 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\mathrm{mm}/\mathrm{s}$, the normal load in a range from $3–30\mathrm{N}$. The test duration was $2–5\mathrm{h}$. The amount of wear was obtained by weighing the pin before and after the test.

Applying the aforementioned model to all contacting asperities, the wear coefficient can be determined in local resolution. Therefore it is assumed, that the asperity $i$ is loaded by a cyclic shear stress ${\tau }^i$, which is the product of the local pressure and the friction coefficient:

(17)${\tau }^i=\mu \frac{f_{\mathrm{sol}}^i}{A_{\mathrm{R}}^i}=\mu p_{\mathrm{sol}}^i$

Thus, the local wear coefficient is instantly a function of the local shear stress. Assuming an incremental wear path in the length of $S_{\mathrm{w}}^i=d$ (element edge length) and an incremental wear volume of $V_{\mathrm{w}}^i=d^3$ (volume of a cubic element), the local relative wear volume now takes the form:

(18)$V_{\mathrm{w}}^i=K^i{d}^3 \to \frac{V_{\mathrm{w}}^i}{d^2}=K^{i}d$

Here, the probability (or wear-) coefficient $K^i$ is considered in local resolution. Rearranging the equation and using the relationship $V_{\mathrm{w}}^i=z_{\mathrm{w}}^i{d}^2$, the local wear depth can be expressed as $z_{\mathrm{w}}^i=K^{i}d$. Taking additionally the number of load cycles $n_{\mathrm{load}}$ into account, the cumulative wear path can be considered. The locally resolved wear depth can finally be given in the form of a vector:

(19)$z_{\mathrm{w}}^i=K^{i}d{n}_{\mathrm{load}} \to z_{\mathrm{w}}=\mathit{K}d{n}_{\mathrm{load}}$

2.3. 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 $h_{\mathrm{c}}$ 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 $\overline{g}$ of the gap distribution is determined. If the deviation between $\overline{g}$ and the film thickness $h_{\mathrm{c}}$ exceeds a specific limit $\epsilon$, the solid overlap is updated (increased in case of $\overline{g}/h_{\mathrm{c}}>1$ and decreased in case of $\overline{g}/h_{\mathrm{c}}<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).

3. Simulation Results

3.1. 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 $F_{\mathrm{N}}=6850\mathrm{N}$. Each disc has a radius of $60\mathrm{mm}$ and is driven independently from the other with different circumferential velocities ($v_1$ and $v_2$). Thus, in the contact area an entrainment velocity of $v=1/2 \left(v_1+v_2\right)=5\mathrm{m}/\mathrm{s}$ and a relative velocity of $v_{\mathrm{rel}}=v_2-v_1=0.5\mathrm{m}/\mathrm{s}$ evolve, depicted in Figure 7.

Disc 1 is additionally equipped with a crowning radius of $100\mathrm{mm}.$ The contact and load parameters lead to a Hertzian pressure of $p_{\mathrm{Hertz}}=1875\mathrm{MPa}$. The root mean squared (RMS) roughness $R_{\mathrm{q}}$ of the discs is $458\mathrm{nm}$. The dynamic viscosity of the lubricant is $\eta =44.2\mathrm{mPas}$. The entrainment velocity is responsible for the formation of an EHL film, which takes a value of $552\mathrm{nm}$, respectively a specific film thickness of $\lambda =h_{\mathrm{c}}/R_{\mathrm{q}}=1.21$. The contact area is discretized into $3010\times 1410$ square elements with an edge length of $d=2\mu \mathrm{m}$. According to the half-space theory, 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.

3.2. 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 Figure 8 and Figure 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 $g$ (in $z$-direction) is given in $\mu \mathrm{m}$ on the ordinate, while the abscissa shows the $x$-direction along the long half-axis. The latter is normalized to the Hertzian contact length $a_{\mathrm{Hertz}}$. On the top right corner, the corresponding undeformed surface is depicted keeping the aspect ratio constant. The numerically calculated film thickness is $h_{\mathrm{c}}=587\mathrm{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\mathrm{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 Figure 10, Figure 11 and Figure 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 $x$-direction along the long half-axis ($a_{\mathrm{Hertz}}$) of the contact zone, while the ordinate shows the $y$-direction along the short half-axis ($b_{\mathrm{Hertz}}$). 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 $x$-direction due to the manufacturing process.

For the case of a dry contact, the sum of all local forces leads to $F_{\mathrm{N},\mathrm{dry}}=F_{\mathrm{N},\mathrm{sol},\mathrm{dry}}\approx F_{\mathrm{N}}$, an EHL part does not exist (see Equation (12)):

(20)$F_{\mathrm{N},\mathrm{sol},\mathrm{dry}}={{\displaystyle \sum }}_{i=1}^{n_{\mathrm{sol}}}{f}_{\mathrm{sol}}^i=\mathrm{6,847}\mathrm{N}, F_{\mathrm{N},\mathrm{EHL},\mathrm{dry}}=0\mathrm{N}$

(21)$F_{\mathrm{N},\mathrm{dry}}=\mathrm{6,847}\mathrm{N}$

Slight deviations from the input values occur due to numerical approximations.

For the mixed lubricated contact, the sum of all local solid forces is:

(22)$F_{\mathrm{N},\mathrm{sol},\mathrm{mix}}={\displaystyle \sum }_{i=1}^{n_{\mathrm{sol}}}{f}_{\mathrm{sol}}^i=603 \mathrm{N}$

Thus, the percentage $P$ of the force, transmitted through asperity contacts is much lower, than in the case of a dry contact: $P=F_{\mathrm{N},\mathrm{sol},\mathrm{mix}}/F_{\mathrm{N}}=8.80\%$. The remaining part of the total force, $F_{\mathrm{N},\mathrm{EHL},\mathrm{mix}}$, has to be transmitted by the lubricant film:

(23)$F_{\mathrm{N},\mathrm{EHL},\mathrm{mix}}=F_{\mathrm{N}}-F_{\mathrm{N},\mathrm{sol},\mathrm{mix}}=6247\mathrm{N}$

In Figure 11 the pressure distribution is shown in the XZ view. The abscissa now shows the long half-axis ($a_{\mathrm{Hertz}}$) and the ordinate shows the normalized pressure $p/p_{\mathrm{Hertz}}$. 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 $P_{\mathrm{flow}}\approx 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 ($b_{\mathrm{Hertz}}$).

3.3. Simulation of Wear Distribution

The simulation of the adhesive wear follows the model presented in Section 2.2 and Section 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\times {10}^{-3}$ (while the minimum value for the wear coefficient was set to $K_{\min }=1.0\times {10}^{-8}$, according to [5]). In Figure 13 the corresponding curves are depicted.

Within this work a friction coefficient of $\mu =0.12$ is assumed such that the calculated wear coefficient lies within an interval of $K^i\in \left[1\times {10}^{-8},1.5\times {10}^{-5}\right]$. The simulated test duration is $60\min$. While the mean values of all local wear coefficients do not differ significantly for the dry and the mixed lubricated case (${\overline{K}}_{\mathrm{dry}}=9.05\times {10}^{-6}$, ${\overline{K}}_{\mathrm{mix}}=8.56\times {10}^{-6}$), the number of asperities in contacts does: $n_{\mathrm{sol},\mathrm{dry}}=\mathrm{18,846}$, $n_{\mathrm{sol},\mathrm{mix}}=1766$. In order to derive a characteristic parameter, describing the intensity of wear for each considered system, the mean value $\overline{K}$ can be multiplied by the respective number of contacting asperities:

(24)${\overline{I}}_{\mathrm{dry}}={\overline{K}}_{\mathrm{dry}}{n}_{\mathrm{sol},\mathrm{dry}}=0.171$

(25)${\overline{I}}_{\mathrm{mix}}={\overline{K}}_{\mathrm{mix}}{n}_{\mathrm{sol},\mathrm{mix}}=0.015$

The factor $\overline{I}$ 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 $z_{\mathrm{w}}$ (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 $x$-direction along the long half-axis ($a_{\mathrm{Hertz}}$) 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 $z_{\mathrm{w},\max ,\mathrm{dry}}=-40.2 \mu \mathrm{m}$, which are not shown within the scale due to a better resolution of the figure.

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 $z_{\mathrm{w},\max ,\mathrm{mix}}=-10.6 \mu \mathrm{m}$.

4. Conclusions 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 ${\nu }_{40}=186\mathrm{cSt}$ at $40°\mathrm{C}$). The investigations were carried out at low relative speeds between the contacting bodies, so that the film thickness ratio was also low ($\lambda =0.07$). Three load cases with axial loads of $25$, $37.5$ and $50\mathrm{kN}$ were defined, resulting in Hertzian pressures of $1028$, $1216$ and $1383\mathrm{MPa}$, respectively. For each load case (LC) a test with a constant shaft speed of $7.5 {\min }^{-1}$ was conducted. Each test duration was $80\mathrm{h}$ ($36,000$ 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 $\mu =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\mathsf{\mu }\mathrm{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\mu \mathrm{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.

Author Contributions

Conceptualization, J.T.T. and G.P.; methodology, J.T.T. and G.P.; software, J.T.T. and M.A.F.; validation, J.T.T. and M.A.F.; formal analysis, J.T.T. and M.A.F.; investigation, J.T.T. and M.A.F.; resources, G.P. and F.P.; data curation, M.A.F.; writing—original draft preparation, J.T.T.; writing—review and editing, M.A.D., F.P. and G.P.; visualization, J.T.T. and M.A.F.; supervision, G.P.; project administration, J.T.T. and F.P.; funding acquisition, J.T.T., F.P. and G.P. All authors have read and agreed to the published version of the manuscript.

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.

List of Symbols

$A_{\mathrm{R}}^i$	Micro contact area at pos. $i$ $\left({\mathrm{m}}^2\right)$	$K^i$	Archard Wear coefficient at pos. $i$
$a_{\mathrm{Hertz}}$	Hertzian contact length $\left(\mathrm{m}\right)$	$P$	Normalized pressure
$a_1,\dots ,a_4$	Parameter (film thickness)	$P_{\mathrm{flow}}$	Normalized flow pressure
$B_n$	Damage propagation after $n$ load cycles	$p$	Pressure $\left(\mathrm{Pa}\right)$
$b_{\mathrm{Hertz}}$	Hertzian contact width $\left(\mathrm{m}\right)$	$p^i$	Pressure at pos. $i$ $\left(\mathrm{Pa}\right)$
$C$	Influence matrix $\left(\mathrm{m}/\mathrm{N}\right)$	$p_{\mathrm{Hertz}}$	Hertzian pressure $\left(\mathrm{Pa}\right)$
$D$	Damage variable	$p_{\mathrm{Hertz}}^i$	Pressure acc. to Hertz at pos. $i$ $\left(\mathrm{Pa}\right)$
$D_{\mathrm{c}}$	Critical damage	$p_{\lim }$	Flow pressure $\left(\mathrm{Pa}\right)$
$D_n$	Damage variable after $n$ load cycles	$p_{\mathrm{sol}}^i$	Solid pressure at pos. $i$ $\left(\mathrm{Pa}\right)$
$d$	Element length $\left(\mathrm{m}\right)$	$R$	Geometry parameter
$\mathrm{d}A$	Cross section $\left({\mathrm{m}}^2\right)$	$R_{\mathrm{q}}$	RMS roughness $\left(\mathrm{m}\right)$
$E$	Young’s modulus $\left(\mathrm{Pa}\right)$	$r$	Asperity radius $\left(\mathrm{m}\right)$
$e_{\mathrm{R}}^{*}$	App. frict. energy density $\left(\mathrm{Nm}/{\mathrm{m}}^3\right)$	$r_{1,1}$, $r_{1,2}$	Radius disc 1, plane 1/2 $\left(\mathrm{m}\right)$
$f^i$	Normal force at pos. $i$ $\left(\mathrm{N}\right)$	$r_{2,1}$ $r_{2,2}$	Radius disc 2, plane 1/2 $\left(\mathrm{m}\right)$
$f_{\mathrm{EHL}}^j$	Solid force at pos. $i$ $\left(\mathrm{N}\right)$	$S_{\mathrm{w}}$	Sliding path $\left(\mathrm{m}\right)$
$f_{\mathrm{sol}}^i$	EHL force at pos. $i$ $\left(\mathrm{N}\right)$	$S_{\mathrm{w}}$	Sliding path vector $\left(\mathrm{m}\right)$
$F_{\mathrm{N}}$	Normal force $\left(\mathrm{N}\right)$	$S_{\mathrm{w}}^i$	Micro sliding path at pos. $i$ $\left(\mathrm{m}\right)$
$F_{\mathrm{sol}}$	Solid force $\left(\mathrm{N}\right)$	$t$	Duration $\left(\mathrm{h}\right)$
$F_{\mathrm{N},\mathrm{sol},\mathrm{dry}}$	Solid force, dry contact $\left(\mathrm{N}\right)$	$U$	Velocity parameter
$F_{\mathrm{N},\mathrm{sol},\mathrm{mix}}$	Solid force, mixed contact $\left(\mathrm{N}\right)$	$V_{\mathrm{w}}$	Wear volume $\left({\mathrm{m}}^3\right)$
$F_{\mathrm{EHL}}$	EHL force $\left(\mathrm{N}\right)$	$V_{\mathrm{w}}^i$	Micro wear volume at pos. $i$ $\left({\mathrm{m}}^3\right)$
$F_{\mathrm{N},\mathrm{EHL},\mathrm{dry}}$	EHL force, dry contact $\left(\mathrm{N}\right)$	$v$	Entrainment velocity $\left(\mathrm{m}/\mathrm{s}\right)$
$F_{\mathrm{N},\mathrm{EHL},\mathrm{dry}}$	EHL force, mixed contact $\left(\mathrm{N}\right)$	$v_1$	Circumferential velocity body 1 $\left(\mathrm{m}/\mathrm{s}\right)$
$G$	Lubrication parameter	$v_2$	Circumferential velocity body 2 $\left(\mathrm{m}/\mathrm{s}\right)$
$g$	Gap $\left(\mathrm{m}\right)$	$v_{\mathrm{rel}}$	Relative velocity $\left(\mathrm{m}/\mathrm{s}\right)$
$g^i$	Gap at pos. $i$ $\left(\mathrm{m}\right)$	$w^i$	Deformation at pos. $i$ $\left(\mathrm{m}\right)$
$\overline{g}$	Mean gap $\left(\mathrm{m}\right)$	W	Load parameter
$H$	Surface hardness $\left(\mathrm{Pa}\right)$	$W_{\mathrm{R}}$	Frictional work $\left(\mathrm{Nm}\right)$
$h_{\mathrm{c}}$	Film thickness $\left(\mathrm{m}\right)$	Z	Number of asperities
$\overline{I}$	Wear intensity factor	$z_{\mathrm{w}}^i$	Wear depth at pos. $i \left(\mathrm{m}\right)$
$i,j,k,l$	Position $i,j,k,l$	$z_{\mathrm{w}}$	Wear depth vector $\left(\mathrm{m}\right)$
$L$	Hardening coefficient	$\eta$	Dynamic viscosity $\left(\mathrm{Pas}\right)$
$M$	Hardening exponent,	$\lambda$	Specific film thickness
	Number of elements in $x$-direction	$\kappa$	Ellipticity
$m$	Normal vector	$\mu$	Friction coefficient
N	Number of elements in $y$-direction	ν	Poisson’s ratio
$N_{\mathrm{c}}$	Critical number of load cycles	$\tau$	Shear stress $\left(\mathrm{Pa}\right)$
$n$, $n_{\mathrm{load}}$	Number of load cycles	${\tau }_{\mathrm{e}}$	Endurance limit
$n_{\mathrm{EHL}}$	Number of non-contacting asperities	${\tau }_{\mathrm{f}}$	True failure stress
$n_{\mathrm{sol}}$	Number of contacting asperities	${\tau }^i$	Shear stress at pos. $i$ $\left(\mathrm{Pa}\right)$
$K$	Archard Wear coefficient	${\displaystyle \sum }\mathrm{d}{A}_{D_n}$	Damaged area $\left({\mathrm{m}}^2\right)$
$K$	Archard Wear coefficient vector
K1, k	Probability factor

Figures and Table

Figure 1. Adhesive wear mechanism: (a) before contact; (b) during contact; (c) after contact [2].

Figure 2. Equilibrium of forces and central film thickness.

Figure 3. Wear mechanism for cubic asperities: (a) before contact; (b) during contact; (c) after contact [2].

Figure 4. Micro cracks and weakened cross section [2,17].

Figure 5. Pin-on-disk tribometer (a) and comparison of experimental and simulation results (b) for the prediction of the wear coefficient [18].

Figure 6. Simulation algorithm. (Left): Simulation of pressure. (Right): Simulation of wear.

Figure 7. Model system.

Figure 8. Gap distribution for dry contact (XZ view).

Figure 9. Gap distribution for a mixed lubricated contact (XZ view).

Figure 10. Pressure simulation results—XY view. (a) Dry contact. (b) Mixed lubricated contact.

Figure 11. Pressure simulation results—XZ view. (a) Dry contact. (b) Mixed lubricated contact.

Figure 12. Pressure simulation results—YZ view. (a) Dry contact. (b) Mixed lubricated contact.

Figure 13. Wear coefficient as a function of the friction coefficient and the (local) pressure.

Figure 14. Wear simulation results for a dry contact (XZ section).

Figure 15. Wear simulation results for a mixed lubricated contact (XZ section).

Figure 16. FE8 test rig (a) [20] and test bearing (b).

Figure 17. Experimental (a) and simulation results (b), FE8 test rig [2]. Kinematics in test bearings (c).

Table 1

Material, geometry and simulation parameters [18,19].

Material Parameters			Crowning disc 2	$r_{2,2}$	$\infty$
Young’s Modulus	$E$	$210\mathrm{GPa}$	Root mean squared (RMS) roughness	$R_{\mathrm{q}}$	$458\mathrm{nm}$
Poisson’s ratio	$\nu$	$0.3$	Simulation Parameters
Flow pressure	$p_{\mathrm{flow}}$	$7.13\mathrm{GPa}$	Normal force	$F_{\mathrm{N}}$	$6850\mathrm{N}$
Cycl. hard. exponent	$M$	$6.22$	Hertzian pressure	$p_{\mathrm{Hertz}}$	$1875\mathrm{MPa}$
Hard. modulus	$H$	$3.44\mathrm{GPa}$	Entrainment velocity	$v$	$5\mathrm{m}/\mathrm{s}$
True failure stress	${\tau }_{\mathrm{f}}$	$2.59\mathrm{GPa}$	Relative velocity	$v_{\mathrm{rel}}$	$0.5\mathrm{m}/\mathrm{s}$
Endurance limit	${\tau }_{\mathrm{e}}$	$0.77\mathrm{GPa}$	Duration	$t$	$60\min$
Critical Damage	$D_{\mathrm{c}}$	$0.5$	Film thickness	$h_{\mathrm{c}}$	$552\mathrm{nm}$
Dyn. viscosity	$\eta$	$44.2\mathrm{mPas}$	Rel. film thickness	$\lambda$	$1.21$
Geometry Parameters			Friction coeff. (solid)	$\mu$	$0.12$
Radius disc 1	$r_{1,1}$	$60\mathrm{mm}$	Element width	$d$	$2\mu \mathrm{m}$
Crowning disc 1	$r_{1,2}$	$100\mathrm{mm}$	Elements in $x$-dir.	$M$	$3010$
Radius disc 2	$r_{2,1}$	$60\mathrm{mm}$	Elements in $y$-dir.	$N$	$1410$