Abstract
A new model was developed for the simulation of the friction coefficient in lubricated sliding line contacts. A halfspacebased contact algorithm was linked with a numerical elastohydrodynamic lubrication solver using the loadsharing concept. The model was compared with an existing asperitybased friction model for a set of theoretical simulations. Depending on the load and surface roughness, the difference in friction varied up to 32 %. The numerical lubrication model makes it possible to also calculate lightly loaded contacts and can easily be extended to solve transient problems. Experimental validation was performed by measuring the friction coefficient as a function of sliding velocity for the stationary case.
Introduction
Friction, lubrication and wear are highly related phenomena. The presence of lubricating oil between contacting surfaces is often beneficial, as it reduces the shear stresses and thereby friction and wear. By contrast, breakdown of the oil film leads to adverse effects. A film failure is associated with surface roughness. In the case of insufficient lubrication, two sliding bodies under normal load come in direct contact first at the highest peaks (asperities), and the lubricant film breakdown is initiated at these spots. Wear is initiated at these contact spots as well. In line with the efforts to reduce energy consumption, significant efforts have been made in the development of friction models. The literature reveals the existence of two general approaches for the calculation of friction in lubricated contacts: rough elastohydrodynamic lubrication (EHL) theorybased methods and loadsharing concept models [1, 2]. In the rough EHL approach, the governing lubricant flow and surface deformation equations are solved numerically assuming a rough surface in both the Reynolds equation as well as for the solid deflection. Despite the significant progress achieved in the development of such models [3–6], certain problems still exist in the “thin film lubrication regime,” such as convergence, accuracy and meshdependence [7]. The selection of the film breakdown criteria in mixed lubrication models has a significant effect on the calculated friction values. So far, there is no consensus on this topic.
On the other hand, loadsharing models, first introduced by Johnson [8], offer advantages of robustness and a relatively simple methodology for the estimation of friction. In this approach, separate smooth surface lubrication and dry rough contact models are employed which simplifies the problem. Asperity contact conditions are estimated using dry rough surface contact models, such as the Greenwood and Williamson model [9]. The lubrication models apply the EHL theory under the assumption of smooth surfaces [10, 11] and consequently are more robust than the rough EHL algorithms, as they are typically based on function fits. The contact and lubrication models are linked through the proportionality of the loads carried by liquid and surface contacts [8].
Since Johnson introduced the loadsharing concept, a large number of researchers employed this approach for the calculation of the friction coefficient in lubricated contacts, see for example [12–14], and many attempts were undertaken to reduce the assumptions, making the models more widely applicable. Most of the researchers addressed the limitations of using the dry rough contact models. Greenwood and Williamson [9] pioneered with the asperitybased statistical approach for predicting surface contact interaction. The surface was represented as a set of fully elastic, independent, hemispherical asperities with equal radii of curvature with a Gaussian distribution of heights. Zhao et al. [15] extended the theory to the case of elastic/plastic and fully plastic deformation regimes. Arbitrary height distributions and individual radii of curvature were introduced in deterministic models [16, 17]. Mentioned models can be categorized as uncoupled [18], since the real surface roughness is replaced by a set of independent asperities of simplified geometry. An overview of existing contact models can be found in Ref. [19].
Based on the uncoupled approach, a number of researchers developed friction and wear models. Gelinck et al. [20] extended Johnson’s model to calculate friction in all lubrication regimes in line contacts. Akbarzadeh et al. [21] included a thermal reduction in the film thickness and viscosity of the lubricant due to heat generated by both the lubricant and asperities using the statistical approach. The statistical elastoplastic asperity contact model was utilized by Masjedi et al. [13] along with the numerical solution of EHL equations for smooth surfaces to construct a curvefit formula of the traction coefficient. Recently, Chang et al. [22] introduced the influence of boundary film tribochemistry on the asperity level to account for variation of the shear strength of boundary films.
On the other hand, also coupled models for the solution of dry rough contact problems exist. These models take the measured surface topography as it is and do not require any additional assumptions on the geometry of asperities or their height distribution. Interactions due to deflections are taken into account automatically. In this type of models, the halfspace approximation theory is widely used [23]. The solution in this case is obtained numerically [24]. Numerical methods for the solution of dry rough contact problems based on halfspace approach were developed by Polonsky [25], Liu [26] and Tian [27].
So far, the loadsharing lubricated friction methodology mostly relied on uncoupled models due to their very limited computational effort. Some estimations state that the assumption of independent asperities employed in such models does not hold when the ratio of the central film thickness to the standard deviation of asperity heights is less than 0.5 [2]. A model based on coupled contact approach was used by Bobach et al. [28] to calculate friction in mixed lubrication. The halfspace theory was employed to precalculate the dry mean contact pressure as a function of deformed gap height by assuming nominally flat surfaces in contact. And further, it was incorporated into the loadsharing concept to find the fraction of the load carried by solid contacts and lubricant. Multilevel multiintegration was used to facilitate the calculations. A loadsharing conceptbased approach which employs the halfspace theory was also implemented by MoralesEspejel et al. [29]. They used Hertz’ theory [23] to calculate the average pressure, and the mean film thickness was obtained from “central film thickness formula,” such as the Dowson and Higginson fit [10] be it with the extension of an equivalent roughness using the perturbation method. However, this makes the approach applicable only to highly stationary, loaded cases. This method was further developed in this paper so that it can also be applied for low loads and transient problems.
The increase in computational power and the development of efficient numerical methods made it possible to consider not only coupled contact models, using the halfspace theory, but also numerical smooth EHL theory, to tackle mentioned limitations. Therefore, in the current paper, the loadsharing conceptbased models are further improved by taking a nonuniform film into account. The model combines the halfspace approximationbased dry contact model with a numerical EHL model for smooth surfaces in the scope of the loadsharing concept. It will be demonstrated that the approach can be employed for the calculation of friction, real contact area and contact pressures in lubricated line contacts. Several simulation tests will be performed, and the results will be compared with existing models and experimental data to validate the developed model.
Problem Formulation and Solution Methods
As mentioned above, in the loadsharing method, separate dry contact and lubrication models are needed. Therefore, these models are described individually followed by a brief description of the algorithm for the friction calculation. Detailed overviews on algorithms of the load sharingbased models can be found elsewhere [14, 20, 30]. In addition, an artificial rough surface generation technique is discussed in this section.
Dry, Rough, HalfSpaceBased Contact Model
An alternative to the classical description of the surface by means of asperities is the direct use of the (measured) height profile. The halfspace approximation can be used to calculate the surface deflection, described by [24]:
where \(u\left( {x,y} \right)\) is the deflection, \(E^{\prime}\) is the composite elastic modulus, and \(p\left( {x^{\prime},y^{\prime}} \right)\) is the contact pressure.
If the pressure is known, then the surface deflection can be found according to Eq. (1). In the loadsharing concept, however, the situation is inverse: The deflection is given and the pressure has to be found. For a given separating distance H _{S} between the surfaces, the deflection is:
where A _{C} is the contact area, \(z\left( {x,y} \right)\) is the surface roughness profile, \(\delta \left( {x,y} \right) = z\left( {x,y} \right)  h_{s} \left( {x,y} \right)\). The contact area A _{C} is actually not known a priori. The boundary condition reads \(p\left( {x,y} \right) \ge 0\) [27]. A schematic diagram of the contact is shown in Fig. 1.
The corresponding system of equations to be solved then reads:
where P _{U} is the hardness of the material, beyond which the yielding is initiated. Once the pressure \(p\left( {x,y} \right)\) is known, the load carried by surface contacts F _{C} can be obtained:
In the system (3), \(p\left( {x,y} \right)\) and A _{C} are unknowns, whereas \(h_{\text{s}} \left( {x,y} \right)\) is related with the hydrodynamic film thickness, as shown in the next section. The pressure \(p\left( {x,y} \right)\) outside of the contact area is zero. The set of equations cannot be solved analytically, and numerical methods have to be used. It also should be noted that the last equation in the system (3) is used to model perfectly plastic behavior of the material [25]. In the case of pure elastic simulations, this condition should be eliminated from the system of the equations.
Computational time necessary for the solution of the system (3) is dominated by the calculation of the integral Eq. (1). This is a convolution integral, and thus, the computational burden can be decreased from \(O\left( {N^{2} } \right)\) to the \(O\left( {N{ \log }N} \right)\) by discrete fast Fourier transformation technique (DCFFT) [26]. In this case, Eq. (1) is rewritten in the following form:
where K is the influence matrix, as described in Appendix 1.
The foundation of the numerical solver of the system (3) is the method developed by Polonsky and Keer [25] and is based on the conjugate gradient method [31]. Their algorithm was developed to solve Eq. (1) for pressure with known applied load. In the current approach, the applied load is unknown. However, there is additional equation for the deflection in the contact area, which is given by the second equation of the system (3). This equation links the dry rough contact solver with the lubrication model through the separating distance \(h_{\text{s}} \left( {x,y} \right)\), which is related to the hydrodynamic film. The separating distance must fulfill the condition that the total lubricant volume between the smooth surfaces is equal to that between rough surfaces, as it was defined by Johnson [8]. Therefore, the method of Polonsky and Keer was adapted to solve the considered. The details of the algorithm are given in Appendix 2.
Numerical Lubrication Model
As mentioned in the introduction, existing numerical methods and increased computational power make it possible to compute the smooth EHL solution in acceptable time. Implementation of the smooth EHL model, rather than a film thickness equation, can extend the applicability of the concept both to a wide range of loads and to transient problems.
The classical system of smooth transient line contact EHL equations is given in the following form:
where μ is the kinematic viscosity, R is the radius of the cylinder, and H and P _{ H } are unknown film thickness and pressure, respectively. F _{H} is the load carried by the lubricant film and is considered to be given. \(U = \left( {U_{1} + U_{2} } \right)/2\), where U _{1} and U _{2} are the velocities of the moving surfaces. With boundary conditions, this system can be solved numerically to obtain H and P _{H}. Once this is solved, the bulk deflection of the surfaces due to hydrodynamic pressure P _{H} is also known and therefore can be taken into account in the calculation of the separation by subtracting it from the initial surface profile \(z'\left( {x,y} \right) = z\left( {x,y} \right)  u_{h} \left( {x,y} \right)\), where \(u_{h} \left( {x,y} \right)\) is defined in the next paragraph.
In a smooth line contact, the hydrodynamic film thickness \(h\left( x \right)\) does not vary in Y direction. However, this does not apply to the roughness \(z\left( {x,y} \right)\). It is therefore necessary to have the hydrodynamic solution defined in Y direction as well. This can be accomplished by defining \(h\left( {x,y} \right) = h\left( x \right),\forall y\). The same is applied to the deflection due to hydrodynamic pressure: \(u_{h} \left( {x,y} \right) = u_{h} \left( x \right), \forall y\) and \(u_{h} \left( x \right) =  \frac{4}{{\pi E^{\prime}}}\smallint P_{\text{H}} \left( {x^{\prime}} \right){ \ln }\left( {x  x^{\prime}} \right){\text{d}}x^{\prime}\).
The most widely used numerical solution techniques in EHL are the multigrid method [32–34] and the differential deflection algorithm [35]. The latter approach solves the system in a coupled manner, which makes it robust and fast to converge. This method was used in the current work to calculate the EHL film thickness.
The hydrodynamic film thickness H is used to calculate the separating distance \(h_{\text{s}} \left( {x,y} \right)\):
where \(\varepsilon\) is a constant and A is the total area. So, the separation \(h_{\text{s}} \left( {x,y} \right)\) has the same form as the film thickness (Fig. 10), with the offset value \(\varepsilon\) which is found from:
where \(\widetilde{A}\) is the lubricated area. Equation (8) follows the definition of the separating distance of Johnson [36], but it was modified here to take into account the nonuniform hydrodynamic film and deflection of the surface due to contact and hydrodynamic pressures. It should be mentioned that the total displacement of the surfaces is a sum of the deflection \(u\left( {x,y} \right)\) due to contact pressure and the deflection \(u_{\text{H}} \left( {x,y} \right)\) due to hydrodynamic pressure. The lefthand side of the equation represents the total volume of the lubricant between smooth surfaces, which is equal to the volume occupied by the pockets formed by the noncontacting part of the surfaces, represented by the righthand side of this equation. It should be emphasized that the separation value \(h_{\text{s}} \left( {x,y} \right)\) depends on the deflection \(u\left( {x,y} \right)\). When the surface is loaded, the surface deflects and the volume formed by noncontacting parts between surfaces changes. Therefore, to preserve the lubricant volume according to Eq. (8), the separation values must be changed. On the other hand, as it can be seen from the system (3), the deflection \(u\left( {x,y} \right)\) depends on \(h_{\text{s}} \left( {x,y} \right)\). Therefore, an equilibrium between \(u\left( {x,y} \right)\) and \(h_{\text{s}} \left( {x,y} \right)\) must be established. This can be accomplished by iteratively solving the system (3) and Eq. (8), similarly to the algorithm of MoralesEspejel et al. [29]. However, it was found to be timeconsuming and a different singleloop approach was developed, which is described in Appendix 2.
Calculation of the Friction Coefficient
For the calculation of the friction coefficient, the load carried by the surface contacts F _{ C } and the hydrodynamic film F _{H} are calculated, with the boundary condition that their sum should be equal to the applied load F _{T} = F _{C} + F _{H}. Description of the iterative procedure for finding loads satisfying such a boundary condition can be found in Ref. [14, 37]. The friction coefficient F _{C} can then be obtained from
where the friction coefficient in boundary lubrication regime F ^{B}_{C} has to be obtained experimentally and F _{sh} is the lubricant shear force, which can be calculated as
where \(\tau\) is the shear stress of the lubricant:
and \(v = v\left( {x,\widetilde{z}} \right)\) is the lubricant speed, given by [38]:
The shear stress was calculated on the upper surface \(\left( {\widetilde{z} = h} \right)\). For the calculation of the pressure–viscosity dependency, the Roelands [38] equation was used.
Generation of Surface Roughness
To illustrate the performance of the model, rough surfaces with different properties were created digitally, using an approach based on Fourier analysis from Hu [39]. The simulation allows variation of two major parameters, i.e., standard deviation of the surface heights σ and the autocorrelation length L. Examples of generated surfaces are shown in Fig. 2. It is clear that with the increase in the autocorrelation length, the surface roughness becomes smoother in the sense that the curvature of peaks increases and that transitions from one peak to another become less abrupt. This parameter can therefore be used as a rough indicator for the characteristic distance between asperities, as well as the sharpness of surface features.
Results and Discussion
The model was employed for the calculation of friction coefficients in stationary lubricated contacts. A set of simulations was performed to validate and explore the behavior of the developed model.
Artificial Surface Simulation Results
For the first simulation case, surfaces consisting of uniformly distributed spherical asperities of equal radius of curvature (30 μm) and height (500 nm) were used, see Fig. 3. Computational grid size contained \(N_{x} \times N_{y} = 2000 \times 2000\) data points with the horizontal spacing A _{X} = A _{Y} = 115 nm. The distance λ between the asperities was varied: 14, 20, 30 and 40 μm. The separating distance H _{S} was considered to be constant throughout the contact, so the nominally flat surfaces in contact were assumed. In the loadsharing concept, the load carried by asperities is of interest, as it determines the fraction of the dry and hydrodynamic friction in the total friction coefficient. This was further explored. When the mentioned surface is loaded against a flat and the halfspace theory is applied, a boundary effect rises. As the asperities interact, the inner asperities are influenced by pressures developed on the surrounding asperities from all directions; however, boundary asperities are influenced only from the inner parts of the surface. Hence, to study the interaction effect, it is necessary to exclude the mentioned effect by setting periodic boundary conditions. The considered surface is shown in Fig. 3.
According to the left picture in Fig. 4, the distance between the asperities does not influence the load carried by a single asperity in the deterministic model. This is as expected since the asperities are assumed to deform independently according to Hertz’ theory. However, this does not hold for the halfspace model. The right picture of Fig. 4 shows that for a fixed separation level, the load carried by the asperity decreases with decreasing distance between the asperities. This is caused by the increase in the interaction between the asperities, as the pressures developed at single asperities deflect the whole surface. The increase in the distance between asperities moves the load curves obtained by the halfspace model toward the deterministic solution due to the decreasing asperity interaction.
Generally, the friction in the fluid film is smaller than the friction in the contacting fraction of the surface (where the load is carried by asperities). Following the loadsharing concept, the friction will be smaller if the ratio of load carried by the fluid film and load carried by the asperities increases.
This leads to the following conclusion: The friction values obtained using the halfspace approximation will in general be lower or equal to those calculated by the deterministic asperitybased approach.
As a next step, more complex surfaces generated digitally were analyzed. The autocorrelation length L was varied to explore the differences between deterministic asperitybased and halfspace based models. In this set of simulations, the film thickness (and consequently, separating distance) was considered constant in the contact area in the new model as it is done in the deterministic approach. This makes it possible to compare the two models. After all, in this case, observed differences between calculated friction values will originate from the difference in the contact models. Periodic boundary conditions were used in the simulations.
Example friction curves are shown in Fig. 5. There is a clear difference between the frictional values calculated with both models. The maximum deviation is observed at approximately 0.1 m/s. This can be explained by means of the loadsharing concept. In the boundary lubrication regime, which corresponds to low velocities, the load carried by the film is small. Since the coefficient of friction in boundary lubrication is predefined, both models will give similar results. With the increase in the velocity, the load carried by the film increases and the load carried by the asperities decreases. Correspondingly, in the halfspace approach, the asperities will carry less load while the same separation is imposed compared to the deterministic model, and thus, the lubricant has to carry a larger fraction of the total load. Therefore, the friction coefficient is lower in the halfspace approach. With a further increase in the sliding velocity, the friction developed by the film starts to dominate the total friction coefficient, and the influence of the direct contact between asperities diminishes. Hence, the difference between deterministic and halfspace approaches becomes smaller as both models use the same EHL model. Finally, in the full film regime, there is no asperity contact, and thus the friction coefficients will match.
In order to further explore the behavior of the halfspace approach and quantify its deviation from the deterministic model, the maximum relative difference was introduced:
where F _{D} and F _{C} are the friction coefficient obtained by the deterministic and new model, respectively. The F ^{M}_{D} is the deterministic friction coefficient at the place where the difference \(\left {f_{\text{d}}  f_{\text{c}} } \right\) has reached a maximum.
The friction curves were calculated for two surfaces with equal standard deviation of heights (240 μm), but with autocorrelation lengths L = 34 μm and L = 57 μm. The maximum relative difference D was then calculated as a function of load and autocorrelation length, Fig. 6. Computational grid size contained \(N_{\text{x}} \times N_{\text{y}} = 816 \times 816\) data points with the horizontal spacing A _{ X } = A _{ Y } = 2 μm. First of all, it is clear that the maximum relative difference D can reach almost 32 % at a contact load of 1 GPa, which is considerable. Secondly, the difference depends on the load. Obviously, by increasing the load, more asperities come into contact, increasing the degree of the deformation, and therefore, most likely, the distance between neighboring asperities in contact will decrease. According to previous considerations, this leads to an increase in “neighboring” effect. Finally, there is a clear influence of the autocorrelation length parameter. With the increase in L, the deviation between the models diminishes. As the autocorrelation length roughly characterizes the distance between the surface peaks, such result is in agreement with conclusions obtained earlier for uniformly distributed asperity cases.
Comparison with Experimental Measurements
Further validation of the new model was accomplished using friction measurements in cylinder on disk sliding tests with mineral oil as a lubricant. The viscosity grade of the oil is 150 cSt at 40° C. The tests were performed at room temperature 25° C. Friction curves as a function of sliding velocity were measured for a number of samples and then averaged. A line contact problem was considered with the properties shown in Tables 1 and 2, where B is the contact width and α is the pressure–viscosity coefficient. The nominal contact pressure in this case is 102 MPa. It should be mentioned that since the pressure is not high, the pressure–viscosity coefficient plays a minor role in the lubricant friction. For the simulation, periodic boundary conditions were set in the direction perpendicular to sliding only, since in this case, a cylindrical geometry is considered.
The surface roughness is shown in Fig. 7. As it can be seen, the surface is quite smooth. Computational grid size of the contact model contained \(N_{x} \times N_{y} = 200 \times 200\) data points with the horizontal spacing A _{ X } = A _{ Y } = 0.45 μm. Computational grid size of the hydrodynamic solver contained N ^{H} = 2000 data points with a^{h} = 0.8 μm. The hydrodynamic grid was taken to be considerably larger than the dry contact area in order to exclude possible starvation effects. It should be noted that since the hydrodynamic and contact models grids are different, the former is interpolated to the contact grid.
Simulated friction curves along with experimental measurement data are shown in Fig. 8. Three simulation curves are shown. The most deviating curve corresponds to fully elastic, deterministic, asperitybased contact model linked with numerical shear stress calculation (thus, in full film regime, measurement and simulation match). The closest to experimental data simulation curves were obtained using the discussed halfspacebased models with fully elastic and elastic, fully plastic, surface deformation regimes.
As it can be seen, the halfspacebased models are in a good agreement with experimental data. Slight effect of plastic regime also can be noted. Assumption of the fully elastic deformation leads to the same effect as the assumption of independent asperity deformation—it overestimates the loads carried by contact spots and as a result overestimates the friction coefficient. Due to smooth surface and light load, the effect of plastic deformation is limited.
From the width of the measurement 90 % confidence interval, it can be noticed that the variance of the experimental data is almost negligible in the full film regime, i.e., the right side of the friction curve, and becomes larger toward the boundary lubrication, i.e., the left side. This variation is caused by the difference in roughness profiles of the samples. This difference does not play any role in the full film regime, as the surfaces are completely separated by the lubricant film, but is influential in the boundary lubrication.
The dry contact pressure profile developed in the boundary lubrication regime is shown in Fig. 9. It can be noticed that the pressure profile reaches the hardness values at some contact spots. At these areas, the surfaces deform in the fully plastic regime, which is modeled using a cutoff pressure in the CGM algorithm as is discussed [25]. The total number of asperities in contact is low, since the applied load is low.
Further, in Fig. 10, the cross sections of the film thickness, separation and the surface are shown in boundary lubrication and full film regimes, corresponding to the speed U = 0.001 m/s and U = 0.3 m/s, respectively. Clearly, in boundary lubrication, the separation cannot prevent the surfaces from contacting, as shown in Fig. 10a. The contact spots generate high total friction coefficient, as shown in Fig. 8. On the other hand, in the full film regime, Fig. 10b, there is enough lubricant film to completely separate surfaces from contact. Under these conditions, the friction is generated only by the shear of the lubricant film and therefore is low compared to the boundary value, see Fig. 8.
Conclusion
A robust and fast numerical algorithm for the calculation of the friction coefficient in lubricated line contact problems was developed based on a halfspace approximation approach and numerical smooth EHL theory using the loadsharing concept. This model does not require additional assumptions on the geometry of asperities and takes into account their interaction during deformation. The use of a numerical smooth EHL solver, rather than a film thickness formula, makes it possible to solve lightly loaded and transient problems in the future.
The model was validated using an existing deterministic asperitybased friction model, artificial surface simulation tests and experimental measurements. The results show that the independent asperity deformation assumption employed in many existing friction models can cause a deviation in predicted friction up to 32 %. This deviation depends on the applied load and surface roughness characteristics. It leads to an overestimation of the loads carried by the asperities, resulting in an overestimation of the friction coefficient. Comparison of the simulation data with experimental measurements showed a good agreement and illustrates the accuracy of this friction model.
Abbreviations
 \(x,y, \widetilde{z}\) :

Spatial coordinates (m)
 \(p\left( {x,y} \right)\) :

Dry contact pressure (Pa)
 \(u\left( {x,y} \right)\) :

Deflection due to pressure p (x, y) (m)
 E _{1}, E _{2} :

Young’s modulus of the cylinder and the substrate (Pa)
 \(\nu_{1} ,\nu_{2}\) :

Poisson’s ratio of the cylinder and the substrate (–)
 \(E^{\prime}\) :

Composite Young’s modulus (Pa) \(\frac{2}{{E^{\prime}}} = \frac{{1  \nu_{1}^{2} }}{{E_{1} }} + \frac{{1  \nu_{2}^{2} }}{{E_{2} }}\)
 \(z\left( {x,y} \right)\) :

Measured surface roughness (m)
 \(h_{\text{s}} \left( {x,y} \right)\) :

Separation distance between the bodies (m)
 \(A_{\text{c}} ,\widetilde{A}, A\) :

Dry, lubricated and nominal Hertzian contact areas, respectively (m^{2})
 \(F_{\text{C}}\) :

Load carried by surface contacts (N)
 K :

Influence matrix (m/pa)
 N :

Numerical grid size (–)
 \(h\left( {x,y} \right)\) :

Hydrodynamic film thickness (m)
 \(P_{\text{H}} \left( {x,y} \right)\) :

Hydrodynamic pressure (Pa)
 \(\mu\) :

Kinematic viscosity of the lubricant (Pa s)
 R :

Radius of the cylinder (m)
 \(F_{\text{H}}\) :

Load carried by the hydrodynamic film (N)
 \(F_{\text{T}}\) :

Total applied load (N)
 \(u_{\text{h}} \left( {x,y} \right)\) :

Elastic deflection due to hydrodynamic pressure \(P_{\text{H}} \left( {x,y} \right)\) (m)
 F _{C} :

Friction coefficient (–)
 F ^{B}_{C} :

Friction coefficient in boundary lubrication (–)
 \(\lambda\) :

Distance between the neighboring asperities (m)
 L :

Autocorrelation length (m)
 B :

Width of the cylinder (m)
 \(\alpha\) :

Pressure–viscosity coefficient (1/GPa)
 P _{U} :

Hardness of the substrate (Pa)
References
 1.
Lugt, P.M., MoralesEspejel, G.E.: A review of elastohydrodynamic lubrication theory. Tribol. Trans. 54, 470–496 (2011)
 2.
Spikes, H.A.: Mixed lubrication—an overview. Lubr. Sci. 9(3), 221–253 (1997)
 3.
Elcoate, C.D., Evans, H.P., Hughes, T.G.: On the Coupling of the elastohydrodynamic problem. J. Mech. Eng. Sci. 212, 307–318 (1998)
 4.
Evans, H.P.: Modelling Mixed Lubrication and the Consequences of Surface Roughness in Gear Applications. Paper presented at the 13th Nordic Symposium on Tribology (Nordtrib 2008), 2008
 5.
Evans, H.P., Snidle, R.W., Sharif, K.J.: Deterministic mixed lubrication modelling using roughness measurements in gear applications. Tribol. Int. 42, 1406–1417 (2009)
 6.
Hu, Y.Z., Zhu, D.: A full numerical solution to the mixed lubrication in point contacts. J. Tribol. 122, 1–9 (2000)
 7.
Holmes, M.J.A., Evans, H.P., Snidle, R.W.: Discussion on elastohydrodynamic lubrication in extended parameter ranges, Part II: load effect. Tribol. Trans. 46, 282–288 (2003)
 8.
Johnson, K.L., Greenwood, J.A., Poon, S.Y.: A simple theory of asperity contact in elastohydrodynamic lubrication. Wear 19(1), 91–108 (1972) doi:10.1016/00431648(72)904450
 9.
Greenwood, J.A., Williamson, G.P.B.: Contact of nominally flat surfaces. Proc. Roy. Soc. Lond. Series A Math. Phys. Sci. 295, 300–319 (1966)
 10.
Dowson, D., Higginson, G.R.: New rollerbearing lubrication formula. Engineering (London) 192, 158–159 (1961)
 11.
Nijenbanning, G., Venner, C.H., Moes, H.: Film thickness in elastohydrodynamically lubricated elliptic contacts. Wear 176, 217–229 (1994)
 12.
Chang, L., Jeng, Y.R.: A mathematical model for the mixed lubrication of nonconformable contacts with asperity friction, plastic deformation, flash temperature, and tribochemistry. J. Tribol. 136, 1–9 (2014)
 13.
Masjedi, M., Khonsari, M.M.: Theoretical and experimental investigation of traction coefficient in linecontact EHL of rough surfaces. Tribol. Int. 70, 179–189 (2014)
 14.
Faraon, I.C., Schipper, D.J.: Stribeck curve for starved line contacts. J. Tribol. 129, 181–187 (2007)
 15.
Zhao, Y., Maietta, D.M., Chang, L.: An asperity microcontact model incorporating the transition from elastic deformation to fully plastic flow. J. Tribol. 122, 86–93 (2000)
 16.
Faraon, I.C.: Mixed lubricated line contacts. PhD thesis., University of Twente (2005)
 17.
Popovici, R.: Friction in wheelrail contact. PhD thesis. In. University of Twente (2010)
 18.
Adams, G.G., Nosovsky, M.N.: Contact modeling—forces. Tribol. Int. 33, 431–442 (2000)
 19.
Bhushan, B.: Contact mechanics of rough surfaces in tribology: multiple asperity contact. Tribol. Lett. 4, 1–35 (1998)
 20.
Gelinck, E., Schipper, D.: Calculation of Stribeck curves for line contacts. Tribol. Int. 33, 175–181 (2000)
 21.
Akbarzadeh, S., Khonsari, M.M.: Thermoelastohydrodynamic analysis of spur gears with consideration of surface roughness. Tribol. Lett. 32, 129–141 (2008)
 22.
Chang, L., Jeng, Y.R.: A mathematical model for the mixed lubrication of nonconformable contacts with asperity friction, plastic deformation, flash temperature, and tribochemistry. J. Tribol. 136, 1–9 (2014)
 23.
Johnson, K.L.: Contact Mechanics. Cambridge University Press, Cambridge (1985)
 24.
Timoshenko, S.P., Goodier, J.N.: Theory of Elasticity, 3rd edn. McGrawHill, New York (1970)
 25.
Polonsky, I.A., Keer, L.M.: A numerical method for solving rough contact problems based on the multilevel multisummation and conjugate gradient techniques. Wear 231, 206–219 (1999)
 26.
Liu, S.: Thermomechanical Contact Analysis of Rough Bodies. PhD thesis., Northwestern University (2001)
 27.
Tian, X., Bhushan, B.: A numerical threedimensional model for the contact of rough surfaces by variational principle. J. Tribol. 118, 33–42 (1996)
 28.
Bobach, L., Beilicke, R., Bartel, D., Deters, L.: Thermal elastohydrodynamic simulation of involute spur gears incorporating mixed friction. Tribol. Int. 48, 191–206 (2012)
 29.
MoralesEspejel, G.E., Wemekamp, A.W., FelixQuinonez, A.: Microgeometry effects on the sliding friction transition in elastohydrodynamic lubrication. Journal of Engineering Tribology 224(7), 621–637 (2010)
 30.
Popovici, R., Schipper, D.J.: Stribeck and traction curves for elliptical contacts: isothermal friction model. Int. J. Sustain. Constr. Des. 4(2) (2013)
 31.
Shewchuk, J.R.: An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. http://www.cs.cmu.edu/~quakepapers/painlessconjugategradient.abstract (1994)
 32.
Venner, C.H., ten Napel, W.E.: Multilevel solution of the elstohydrodynamically lubricated circular contact problem. Part 2: smooth surface results. Wear 152, 369–381 (1992)
 33.
Zhu, D., Hu, Y.Z.: A computer program package for the prediction of EHL and mixed lubrication characteristics, friction, subsurface stresses and flash temperatures based on measured 3D surface roughness. Tribol. Trans. 44, 383–390 (2008)
 34.
Zhu, D., Liu, Y., Wang, Q.: On the numerical accuracy of rough surface EHL solution. Tribol. Trans. 57, 570–580 (2014)
 35.
Hughes, T.G., Elcoate, C.D., Evans, H.P.: Coupled solution of the elastohydrodynamic line contact problem using a differential deflection method. J. Mech. Eng. Sci. 214, 585–598 (2000)
 36.
Johnson, K.L., Spence, D.I.: Determination of gear tooth friction by disc machine. Tribol. Int. 24(5), 269–275 (1991)
 37.
Liu, Q., Ten Napel, W., Tripp, J.H., Lugt, P.M., Meeuwenoord, R.: Friction in highly loaded mixed lubricated point contacts. Tribol. Trans. 52(3), 360–369 (2009)
 38.
Szeri, A.Z.: Fluid Film Lubrication. Cambridge University Press, Cambridge (1998)
 39.
Hu, Y.Z., Tonder, K.: Simulation of 3D random rough surface by 2D digital filter and fourier analysis. Int. J. Mach. Tools Manuf 32(1), 83–90 (1992)
Acknowledgments
The authors would like to thank SKF Engineering & Research Centre, Nieuwegein, the Netherlands, and Bosch Transmission Technology BV, Tilburg, the Netherlands, for providing technical and financial support. We would like to thank Mr. A.J.C de Vries, SKF Director Group Product Development for his permission to publish this work. This research was carried out under project number M21.1.11450 in the framework of the Research Program of the Materials innovation institute (M2i).
Author information
Affiliations
Corresponding author
Additional information
This article is part of the Topical Collection on STLE Tribology Frontiers Conference 2014.
Appendices
Appendix 1
After discretization of the governing equations on a numerical grid Ω, the influence matrix K is given by a following expression [23]:
where A _{ X } and A _{ Y } are the half mesh size in direction X and Y correspondingly, X _{ I }, \(y_{j} \in\Omega\).
Appendix 2
The process starts with the initialization of the pressure \(p\left( {x,y} \right)\), which can be arbitrary or obtained from the previous iteration step. The area of contact A _{C} is then defined as the sum of those areas where \(p\left( {x,y} \right) > 0\). The hydrodynamic area is the total area subtracted the contact area \(\widetilde{A} = A \cap \overline{A}_{\text{c}}\). The corresponding deflection \(u\left( {x,y} \right)\) is found according to Eq. (1) using the DCFFT technique. The initial separation value H _{ S } is calculated through the Eq. (8). Auxiliary variables G _{OLD} = 1 and \(\delta = 0\) are defined.
Iterations start, and first the gap distribution is obtained:
This gap is positive in areas of contact and negative in lubricated areas. Note that here U _{ IJ } are positive. For the gap G _{ IJ }, new variable G is introduced:
And a new conjugate direction is computed:
Convolution of the influence matrix and T _{ IJ } is then obtained:
Mean value of R _{ IJ } is then adjusted as
where N _{C} is the number of points in the contact area A _{C}. The length of the step in the direction T _{ IJ } is then:
The pressure is stored as P ^{OLD}_{ IJ } = P _{ IJ } and the new pressure profile can be found from:
The negative values of P _{ IJ } are set to zero. Further, it is necessary to include into consideration the points, where pressure is zero, but surfaces overlap. The fact that they overlap indicates that they are in contact, and thus, the pressure must be positive. The surfaces are overlapping, if G _{ IJ } > 0. The following area is introduced
If \(A_{ol} = \emptyset\), \(\delta\) is set to unity; otherwise, it is set to 0. In this area, the pressure is then corrected accordingly as:
The pressure obtained using this equation is always positive as \(\tau\) is always negative [25]. After that, the deflection and the new value for separation are calculated according to the first of the Eqs. (3) and (8) correspondingly. The iterations continue until the difference in separation level between two successive iterations is less than the desired value \(\varepsilon_{0}\):
In rough surface contact problems, the plastic limit is usually set as an upper boundary for the pressure, where the material starts yielding. In the described algorithm, it is possible to impose this condition the same way as the condition on the positivity of the contact pressure on each iteration, \(p_{ij} \le P_{\text{u}}\), where P _{U} is the plastic limit.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Akchurin, A., Bosman, R., Lugt, P.M. et al. On a Model for the Prediction of the Friction Coefficient in Mixed Lubrication Based on a LoadSharing Concept with Measured Surface Roughness. Tribol Lett 59, 19 (2015). https://doi.org/10.1007/s112490150536z
Received:
Accepted:
Published:
Keywords
 Mixed lubrication
 Friction modeling
 Halfspace approximation