Analytical Formula for the Ratio of Central to Minimum Film Thickness in a Circular EHL Contact

Prediction of minimum film thickness is often used in practice for calculation of film parameter to design machine operation in full film regime. It was reported several times that majority of prediction formulas cannot match experimental data in terms of minimum film thickness. These standard prediction formulas give almost constant ratio between central and minimum film thickness while numerical calculations show ratio which spans from 1 to more than 3 depending on M and L parameters. In this paper, an analytical formula of this ratio is presented for lubricants with various pressure–viscosity coefficients. The analytical formula is compared with optical interferometry measurements and differences are discussed. It allows better prediction, compared to standard formulas, of minimum film thickness for wide range of M and L parameters.


Introduction
Prediction of film thickness represents an important step necessary for up-to-date design of lubricated concentrated contact.Analytical prediction formula can provide simple and rapid estimate of film thickness in elastohydrodynamically lubricated (EHL) contact.Minimum film thickness is used for calculation of film parameter to judge if machine will operate in full or mixed lubrication regime.
Since the pioneering work of Hamrock-Dowson in 1977 [1], numerous prediction formulas have been published in the literature [2][3][4][5][6].An overview of the prediction formulas can be found in References [7][8][9].A common approach is to make full numerical simulations of a contact for some range of conditions and solve regression analysis of film thickness results on operating conditions.These regressions were proven to be surprisingly accurate (considering the time of original formulation), if conditions of contact inlet are isothermal Newtonian.This is especially true for central film thickness, but generally not the case of minimum film thickness [10][11][12][13][14].
Several authors have showed that the ratio between central and minimum film thickness vary significantly with operating parameters [5,[15][16][17].Moreover, trend of the ratio is not monotonous depending on L parameter, i.e., reaching a local maximum at some value of L parameter.Classic prediction formula defined as a power function of operating conditions cannot represent such phenomenon.For better minimum film thickness formula, it is necessary to search for analytical form that can describe the ratio between central and minimum film thickness.Recently, a parameterization for the ratio of central to minimum film thickness was found for slender contacts (contact size in entrainment direction is much larger than in perpendicular direction) [18].
In this paper, the ratio between central film thickness and minimum film thickness is studied depending on M and L dimensionless parameters.New analytical formula allowing rapid prediction of minimum film thickness in circular point contact is suggested based on numerical simulations.
The new formula is compared to ratio that come from Hamrock-Dowson equations and optical interferometry measurements.

Numerical Calculations
An isothermal Newtonian numerical solution of elastohydrodynamically lubricated (EHL) point contact is assumed.Mathematical model consists of Reynolds equation with boundary conditions p(x a , y) = p(x b , y) = p(x, y a ) = p(x, y b ) = 0 and the cavitation condition p(x, y) ≥ 0, film thickness equation and force balance equation Roelands pressure-viscosity and Dowson-Higginson pressure-density relations have been considered ρ(p) = ρ 0 5.9•10 8 + 1.34p 5.9•10 8 + p The mathematical model in Equations ( 1)-( 3) was solved by the multilevel multi-integration technique [19].Second order discretization was assumed.Several mashes with number of grid points from 129 × 129 to 1025 × 1025 were used and differences between results were monitored to ensure results free of significant grid effect.
Pressure-viscosity response is important for film forming capability of a contact, which is described well by reciprocal asymptomatic isoviscous pressure coefficient α* or film pressure-viscosity coefficient α film which is proportional to the α* evaluated up to maximum inlet zone pressure [20].This pressure is, according to Bair et al. [20], considered to be equal to 3/α*.
The α film pressure-viscosity parameter is used in the present analytical model of central to minimum film thickness ratio.

Film Thickness Measurement
Film thickness measurements were done on a custom developed ball-on-disk optical tribometer (Brno University of Technology, Brno, Czech republic) [12].This device measures film thickness in an elastohydrodynamic contact between steel ball and glass (sapphire) disk based on optical interferometry principle.Interferograms are formed between light beams reflected from a ball surface and bottom side of the glass (sapphire) disk coated by thin layer of chromium to ensure high interference contrast.Interferograms are evaluated by Thin film colorimetric interferometry.
The measurement method uses robust calibration procedure which establishes film thickness/color calibration from monochromatic and chromatic interferogram of a static contact.The calibration procedure and method are described in References [21,22].The present configuration is without spacer layer which enables evaluation of film thickness in a range of 0-800 nm with resolution better than 1 nm.
Film thickness was measured for steel on glass and steel on sapphire configurations with reference fluid tri (2-ethylhexyl) trimellitate (TOTM, Sigma-Aldrich, Darmstadt, Germany) as a lubricant.Rheological properties of this fluid were measured in [23].Basic rheology parameters for test temperature are in Table 1.Temperature was measured in contact inlet by thermocouple calibrated to more accurate thermistor sensor.Properties of contact bodies are in Table 2.

Results and Discussion
A series of EHL contacts with parameters listed in Table 3 were numerically simulated.A convergence better than 10 −4 was required.The ratio of central to minimum film thickness h c /h min was evaluated.Contact simulations were done on square grids with different number of points in a range from 129 × 129 to 1025 × 1025.A difference between results on coarse and finer meshes were evaluated and accepted only when data were with relative difference < 1% for M < 750 and < 3% for M ≥ 750 to ensure they were free of mesh density influence.All calculated data of ratios are listed in the Tables A1-A3.Film thickness ratios were compared to published results in [17,18].The average difference (Table 4) is 4.4% and 3.5%, respectively.Because it was found in present results that film thickness ratio depends on pressure-viscosity coefficient, besides M and L parameters, certain part of the deviations can come from differences in this coefficient.The values of pressure-viscosity, for which the film thickness ratios were obtained, are not clear in the literature sources.Figure 1 shows a domain of current simulations, operating points of experiments and domain of parameters on which Hamrock-Dowson film thickness equations were established.These equations are one of the most employed in EHL field, therefore, it was chosen for comparison.According to [1], the central and minimum film thickness can be predicted by Figure 1 shows a domain of current simulations, operating points of experiments and domain of parameters on which Hamrock-Dowson film thickness equations were established.These equations are one of the most employed in EHL field, therefore, it was chosen for comparison.According to [1], the central and minimum film thickness can be predicted by Figure 2 shows a dependency of film thickness ratio hc/hmin−1 on M parameter for three fixed L values.Figure 3 presents a dependency of film thickness ratio hc/hmin on L parameter for three fixed M values.Both plots are for αfilm = 20.6 GPa −1 .The dependency of the film thickness ratio monotonically rises on M while there is a local maximum in dependency on parameter L. This local maximum is around L = 5 and seems to slightly change with increasing M parameter.The ratio approaches 1 for low L and M values.Since axes in Figure 2 are in log scale, the points follow power trend with parameter M. When x axis in Figure 3 is transformed to log scale, the points follow close to quadratic trend.Therefore, the complete dataset was fitted to Equation (11), where a, b, c and d are fitting constants.Figure 2 shows a dependency of film thickness ratio h c /h min −1 on M parameter for three fixed L values.Figure 3 presents a dependency of film thickness ratio h c /h min on L parameter for three fixed M values.Both plots are for α film = 20.6 GPa −1 .The dependency of the film thickness ratio monotonically rises on M while there is a local maximum in dependency on parameter L. This local maximum is around L = 5 and seems to slightly change with increasing M parameter.The ratio approaches 1 for low L and M values.
Figure 1 shows a domain of current simulations, operating points of experiments and domain of parameters on which Hamrock-Dowson film thickness equations were established.These equations are one of the most employed in EHL field, therefore, it was chosen for comparison.According to [1], the central and minimum film thickness can be predicted by ℎ = 2.69  . . .
1 − 0.61e .Figure 2 shows a dependency of film thickness ratio hc/hmin−1 on M parameter for three fixed L values.Figure 3 presents a dependency of film thickness ratio hc/hmin on L parameter for three fixed M values.Both plots are for αfilm = 20.6 GPa −1 .The dependency of the film thickness ratio monotonically rises on M while there is a local maximum in dependency on parameter L. This local maximum is around L = 5 and seems to slightly change with increasing M parameter.The ratio approaches 1 for low L and M values.Since axes in Figure 2 are in log scale, the points follow power trend with parameter M. When x axis in Figure 3 is transformed to log scale, the points follow close to quadratic trend.Therefore, the complete dataset was fitted to Equation (11), where a, b, c and d are fitting constants.Since axes in Figure 2 are in log scale, the points follow power trend with parameter M. When x axis in Figure 3 is transformed to log scale, the points follow close to quadratic trend.Therefore, the complete dataset was fitted to Equation (11), where a, b, c and d are fitting constants.The effect of different reduced radii of curvature, elastic properties or lubricant ambient viscosity on the ratio was checked by numerical simulations and no effect was found; therefore, the results are considered independent of chosen geometry and elasticity.On the other hand, it was found that film thickness ratio is sensitive to pressure-viscosity coefficient; simulations were done for three values (Table 3) in a range from 8.7 to 32.7 GPa −1 .It corresponds to the pressure-viscosity coefficient according definition in Equation ( 8) published by Bair [23].It was shown that it is able precisely capture relation of viscosity on pressure necessary for film thickness formation in EHL contacts.According Table 3, the values of αfilm are not far from α* coefficient.The analytical model fitted to the results is shown in Figure 4 in a form of contour plot on the left side and residuals from fitting in the right plots.Final form of analytical model of film thickness ratio which depends on M and L parameters and αfilm pressure-viscosity coefficient is given by Equation (12).Quality parameters of the fits shown in Figure 4 are listed in Table 5; root mean square of the error is 0.03-0.04which represents about 2% of average film thickness ratio.Equation (12) was plotted as a model in Figures 2-5.
where αfilm is in GPa −1 .The equation is repeated together with list of assumptions and conditions for which it was obtained in Appendix B. According Equation ( 12) and Figure 2, the pressure-viscosity coefficient modifies constant in power trend on M parameter, but the power slope remains the same.In terms of dependency on L parameter, the pressure-viscosity coefficient shifts the position of maximum peak, as it is shown in Figure 5.The effect of different reduced radii of curvature, elastic properties or lubricant ambient viscosity on the ratio was checked by numerical simulations and no effect was found; therefore, the results are considered independent of chosen geometry and elasticity.On the other hand, it was found that film thickness ratio is sensitive to pressure-viscosity coefficient; simulations were done for three values (Table 3) in a range from 8.7 to 32.7 GPa −1 .It corresponds to the pressure-viscosity coefficient according definition in Equation ( 8) published by Bair [23].It was shown that it is able precisely capture relation of viscosity on pressure necessary for film thickness formation in EHL contacts.According Table 3, the values of α film are not far from α* coefficient.The analytical model fitted to the results is shown in Figure 4 in a form of contour plot on the left side and residuals from fitting in the right plots.Final form of analytical model of film thickness ratio which depends on M and L parameters and α film pressure-viscosity coefficient is given by Equation (12).Quality parameters of the fits shown in Figure 4 are listed in Table 5; root mean square of the error is 0.03-0.04which represents about 2% of average film thickness ratio.Equation (12) was plotted as a model in Figures 2-5.
where α film is in GPa −1 .The equation is repeated together with list of assumptions and conditions for which it was obtained in Appendix B. According Equation ( 12) and Figure 2, the pressure-viscosity coefficient modifies constant in power trend on M parameter, but the power slope remains the same.In terms of dependency on L parameter, the pressure-viscosity coefficient shifts the position of maximum peak, as it is shown in Figure 5.       Next to the simulations, a film thickness was measured at three levels of Hertzian pressure.It was a contact of steel against glass for 26 and 112 N giving Hertzian pressure of 0.493 and 0.799 GPa, respectively, and a contact of steel with sapphire loaded by 63 N giving 1.186 GPa.Central film thickness was evaluated and corrected for refractive index change with pressure, as in [10], while minimum film thickness was left without correction since the place of minimum is expected to be in the area of low pressure.A range of speeds from 100 to 1200 mm/s was measured and h c /h min ratio was calculated.
In Figure 6, there is a comparison of measured ratios with predictions given by present model (Equation ( 12)) and prediction according Hamrock-Dowson formulas (Equations ( 9) and ( 10)).The Hamrock-Dowson formulas gives almost constant h c /h min ratio which is close to measured value only for 0.5 and 1.2 GPa and small range of speed around 500 mm/s.Present model gives good prediction of trends in all cases.Quantitatively, there is very good agreement for 1.2 GPa and not as good agreement for 0.5 and 0.8 GPa.Average differences between measurement and predictions are listed in Table 6.It shows that present formula predicts the ratio within 11%, while Hamrock-Dowson is off by up to 37%.Ratios calculated from measured data having uncertainties can produce a large error when thin films are studied.Therefore, conditions with h min film thickness above 60 nm and h c film thickness above 140 nm were considered.Measurements were repeated three times for each pressure; ratios were evaluated independently and averaged.Average standard deviation of measured ratios is 0.028 which corresponds to 1.3%.Measurement has another error which comes from the fact that minimum film thickness is evaluated as a global minimum in a contact.Every ball has a certain roughness (3 nm RMS in present measurement); therefore, the evaluated h min is practically always picked on a peak of roughness.As a result, measured minimum film thickness tends to be systematically evaluated as lower than the true value.Typical in-contact RMS value of roughness was 0.6 nm, therefore, distance of mid plane to highest peaks was about 1.8 nm which makes average impact on film thickness ratio 0.03, i.e., much less then observed difference.
Lubricants 2018, 6, x FOR PEER REVIEW 7 of 11 Next to the simulations, a film thickness was measured at three levels of Hertzian pressure.It was a contact of steel against glass for 26 and 112 N giving Hertzian pressure of 0.493 and 0.799 GPa, respectively, and a contact of steel with sapphire loaded by 63 N giving 1.186 GPa.Central film thickness was evaluated and corrected for refractive index change with pressure, as in [10], while minimum film thickness was left without correction since the place of minimum is expected to be in the area of low pressure.A range of speeds from 100 to 1200 mm/s was measured and hc/hmin ratio was calculated.
In Figure 6, there is a comparison of measured ratios with predictions given by present model (Equation ( 12)) and prediction according Hamrock-Dowson formulas (Equations ( 9) and ( 10)).The Hamrock-Dowson formulas gives almost constant hc/hmin ratio which is close to measured value only for 0.5 and 1.2 GPa and small range of speed around 500 mm/s.Present model gives good prediction of trends in all cases.Quantitatively, there is very good agreement for 1.2 GPa and not as good agreement for 0.5 and 0.8 GPa.Average differences between measurement and predictions are listed in Table 6.It shows that present formula predicts the ratio within 11%, while Hamrock-Dowson is off by up to 37%.

Case 1
Case 2 Case 3 Measurement pressure 0.5 GPa 0.8 GPa 1.2 GPa H&D formulas (Equations ( 9) and ( 10)) 12.9% 36.5% 9.7% Present formula (Equation ( 12)) 8.8% 10.6% 2.2% Ratios calculated from measured data having uncertainties can produce a large error when thin films are studied.Therefore, conditions with hmin film thickness above 60 nm and hc film thickness above 140 nm were considered.Measurements were repeated three times for each pressure; ratios were evaluated independently and averaged.Average standard deviation of measured ratios is 0.028 which corresponds to 1.3%.Measurement has another error which comes from the fact that minimum film thickness is evaluated as a global minimum in a contact.Every ball has a certain roughness (3 nm RMS in present measurement); therefore, the evaluated hmin is practically always picked on a peak of roughness.As a result, measured minimum film thickness tends to be systematically evaluated as lower than the true value.Typical in-contact RMS value of roughness was 0.6 nm, therefore, distance of mid plane to highest peaks was about 1.8 nm which makes average impact on film thickness ratio 0.03, i.e., much less then observed difference.8) and ( 9)) and present analytical model given by Equation ( 12).8) and ( 9)) and present analytical model given by Equation (12).
As shown in Figures 3 and 5, present model (Equation ( 12)) cannot fully capture the highest values of simulated ratios.The highest fit residuals are between L parameter 5 and 7 where experimental conditions of 0.5 and 0.8 GPa are localized (see Figure 1).Nevertheless, these local residuals of the fit are up to 0.08, thus they cannot fully explain differences in measurements.
In present study, Roelands model for pressure-viscosity relationship and Dowson-Higginson relation of density dependency on pressure were used.These models have limitations.Various lubricants can behave differently, especially under high pressure.Nevertheless, it is commonly expected that low pressure rheology is important for film forming capabilities.Furthermore, real lubricant compressibility influences central film thickness, which affects the film thickness ratio.

Conclusions
The ratio between central and minimum film thickness in a circular EHL contact varies between 1 and 3.16 for considered range of conditions.This trend cannot be captured by widely used prediction formulas employing (monotonous) power function, since dependency on L parameter has a local extreme at L = exp(3/α film 0.2 ).Isothermal Newtonian numerical simulations were used to compute the film thickness ratio for wide range of M and L parameters and three pressure-viscosity coefficients.A new analytical formula for fitting of the ratio is presented.

Figure 1 .
Figure 1.Present calculation domain compared to operating points of experiments and domain for regression of Hamrock-Dowson formulas.

Figure 1 .
Figure 1.Present calculation domain compared to operating points of experiments and domain for regression of Hamrock-Dowson formulas.

Figure 1 .
Figure 1.Present calculation domain compared to operating points of experiments and domain for regression of Hamrock-Dowson formulas.

Figure 2 .
Figure 2. Dependency of h c /h min ratio on M parameter for α film = 20.6 GPa −1 .

Figure 3 .
Figure 3. Dependency of h c /h min ratio on L parameter for α film = 20.6 GPa −1 .

Figure 5 .
Figure 5. Dependency of hc/hmin ratio on L for M = 750 and all three simulated αfilm.

Figure 5 .
Figure 5. Dependency of hc/hmin ratio on L for M = 750 and all three simulated αfilm.

Figure 5 .
Figure 5. Dependency of h c /h min ratio on L for M = 750 and all three simulated α film.

Table 2 .
Properties of contact bodies.

Table 3 .
Range of numerically simulated conditions.

Table 4 .
Comparison of present h c /h min results to published ones.

Table 5 .
Parameters of fitting goodness.

Table 5 .
Parameters of fitting goodness.

Table 6 .
Average difference of measurement and h c /h min ratio estimation by Hamrock-Dowson (H&D) and present analytical formula.

Table 6 .
Average difference of measurement and hc/hmin ratio estimation by Hamrock-Dowson (H&D) and present analytical formula.

Table A2 .
Final regression of numerical simulations has root mean square error of 0.036, i.e., 1.9%.A comparison with measurement showed good trend agreement and quantitatively smaller difference than Hamrock-Dowson formulas.This new formula together with published formula for central film thickness can be used for minimum thickness prediction.Ratios of central film thickness to minimum film thickness in a point contact for α film = 20.6 GPa −1 .

Table A3 .
Ratios of central film thickness to minimum film thickness in a point contact for α film = 32.7 GPa −1 .

Table A4 .
List of contact assumptions and the range of conditions in the numerical analyses.Isothermal Newtonian (without shear thinning effects) conditionsCircular point contact Smooth surface Lubricant rheology governed by Roelands pressure-viscosity and Dowson-Higginson pressure-density relationships M parameter from 2 to 1000 L parameter from 1 to 30 Pressure-viscosity coefficient α film from 8.7 to 32.7 GPa−1