CFD modeling compared to temperature and friction measurements of an EHL line contact

In this paper, predictions from CFD modeling are compared against measurements of surface temperatures and friction for an EHL line contact lubricated with the fluid Santotrac 50. Two slide-to-roll-ratios (SRR), 50% and 100%, and entrainment velocities ranging from 0.211 to 1.13m/s are considered. Very good agreement is shown for the 50% SRR cases, with only a 3% deviation in friction coefficient values. At 100% SRR, the deviation in friction increases to 3–7% which is attributed to deficiencies in the modeling approach with regard to shearthinning. The temperature profiles agree reasonably well at 50% SRR and show larger deviations at 100% SRR. For all cases, the formation of a shear-band in the center of the fluid film is predicted. This is very pronounced for 100% SRR, although likely to be over-estimated by this CFD-approach. The data presented here serve as a basis from which further refinements in the modeling and measurements shown can be made.


Introduction
Elastohydrodynamic lubrication (EHL) occurs in non-conformal contacts in which component surfaces deform elastically, to give an appreciable area of contact, while lubricant viscosity increases as a result of very high local pressure. Considerable research has focused on modeling this form of lubrication and methods most commonly involve applying Reynolds' equation while taking into account effects such pressure on viscosity, surface deformation and lubricant temperature rise. Cavitation is also accounted for; usually by preventing the pressure falling below zero when solving Reynolds' equation. This approach has been highly successful in predicting EHL behavior; however it relies on a number of assumptions and therefore has certain limitations.
An alternative approach for solving EHL is to use computational fluid dynamics (CFD) and this has a number of advantages. Specifically, properties such as pressure, viscosity and temperature can vary across the film and also the whole contact can be modeled including the far inlet and the cavitation zone, rather than just the approximately flat central contact region. Important examples where CFD has been used to solve EHL include work by Schäfer et al. [1] and Almquist and Larson [2] [3], who validated the CFD method against Reynolds solvers and introduced thermal and transient behavior. Hartinger et al. [4] also validated a CFD-approach against a Reynolds solver and presented thermal EHL solutions. Bruyere et al. [5] presented a CFD approach for the fluid and an FEM approach for the solid, using previously published CFD-data as validation [4] and found satisfying agreement (maximum film thickness deviation between the two series of results lower than 2%) for a non-Newtonian, thermal case. One case they studied with zero entrainment velocity showed interesting results, explaining the dimple geometry of the film-thickness. Lohner et al. [6] introduced an approach implemented in COMSOL with a generalized Reynolds equation, discretizing the film-thickness with several layers that coupled with an FEM approach for the solid. They validated their solution against previously published CFD-data [4] and found very good agreement for non-Newtonian thermal and isothermal cases. Furthermore they extended their code to transient cases and validated it against moving surface indentations. A comparison against friction measurements was also presented showing good correlation for the friction coefficient. Feldermann et al. [7] also published an CFD approach validated against a Reynolds solver while Hajishafiee et al. [8] developed a CFD approached for fluid and FEM for solid in the OpenFOAM package. They also validated this against Hartinger et al.'s previously published CFDdata [4] and found good agreement for thermal and isothermal test cases. They achieved high numerical stability with contact pressures of over 3 GPa and showed flow details in these high-pressure regimes.
None of the above work validated CFD predictions against measured temperature predictions; all such validations have been made against measured film thickness or overall friction. In this paper we present such an experimental validation of the CFD approach developed by Hartinger et al. [4]. Although this solver is now somewhat dated, especially by still using analytical models for the solid for deformation and temperature, we believe that it can still serve as a well-validated base for further improvements in this field. On the experimental side, surface temperatures are measured by infrared microscopy and the corresponding friction factors are also measured. These measurements are compared against the CFD solutions and selected cases are analyzed in more detail.

Experimental details
This section outlines the technique used to measure the surface temperatures within the EHL contact. It should be noted that the technique has recently been refined to enable the temperature of the oil film itself to be measured three-dimensionally [9]. However this advanced technique has not been used in the current study, since the aim was only to compare surface temperatures with those predicted from CFD.

Apparatus
The experimental setup for this study was similar to that used to measure contact temperatures previously [10], [11]. The main experimental difference regards the contact geometry, which is elliptical in the current work, compared to a point contact studied previously. To produce an elliptical EHL contact, a steel spherical roller was loaded against a sapphire disc using a conventional EHL ring (manufactured by PCS Instruments, Acton, UK). Sapphire (Al 2 O 3 ) was used for the disc material due to its high transmittance of infrared. Using this equipment, both the rotation of the roller and disc can be driven independently to produce a range of slide-roll ratios. The roller was partially submerged in a temperature-controlled lubricant bath ( ± 0.5°C) so that oil was entrained into the contact during rotation. As shown in Fig. 1, this equipment was located beneath an infrared microscope, which focussed through the disc onto the contact to detect emitted radiation. The microscope consists of an infrared camera (manufactured by Flir Systems, Oregon, USA), fitted with a custom built 5 × infrared lens. The 320 × 256 focal plane array of the camera, combined with this magnification, gives a resolution of 6 μm per pixel.
The radiation reaching the camera is emitted from four sources: 1) the roller surface 2) the lower surface of the disc 3) the bulk of the sapphire disc (background radiation) 4) the oil film In this study measurements of the temperature of the roller surface and the lower surface of the disc were required. Therefore the radiation from these components had to be separated from each other and from radiation emitted from the oil itself and from the bulk of the sapphire disc. Separation was carried out as follows. The radiation component from the oil itself was eliminated by locating a cut-off filter in front of the camera lens to remove wavelengths below 3.75 μm. This cut-off wavelength was used since hydrocarbon lubricants emit primarily over a short band of wavelengths below this value [12], [13]. In order to distinguish between the three remaining radiation components, three separate sapphire discs were used. These discs have different surface coating applied to their lower surfaces as shown in Fig. 2.
An aluminium coated disc was used so that the radiation from the bulk of the sapphire could be isolated. This is possible since aluminium, having an emissivity of ∼0.04, is a very poor emitter of infrared. Therefore, when this disc is used, the radiation reaching camera originates almost entirely from the bulk of the sapphire disc itself. Background radiation obtained in this way was subtracted from subsequent measurements so that radiation components from roller and disc surfaces could be isolated. A chromium coating was used in order to measure the temperature of the disc's lower surface. The chromium, with a thickness of ∼120 nm, is completely opaque to infrared and has an emissivity of ∼0.4. Finally, an uncoated disc was used in order to detect radiation from the steel surface of the roller. For each test condition, measurements were made with each of the three sapphire discs. This differs from previous work, where a single disc, divided into three sections, was used [10].
A calibration procedure was carried out to correlate emission measurements with contact temperature. This was achieved by setting up a ball on disc contact in pure rolling, and monitoring the level of radiation received by the camera as the temperature of the lubricant bath was gradually increased. Under such pure rolling conditions, the temperature rise due to lubricant entrainment can be considered negligible [11]; so that the contact temperature equals that of the bulk lubricant within the bath.
For these calibrations, the background radiation (measured using the aluminium disc), was subtracted from measurements of the uncoated and chromium coated discs to give curves of roller and disc radiation components versus contact temperature. During subsequent tests, carried out under sliding conditions, these calibration curves were then used to obtain temperature values from radiation measurements.
By using moving heat source theory, the measured temperature maps can be converted into shear stress, and in this way it has been possible to validate this infrared microscopy method against friction measurements [10]. The method is sensitive to temperature increases as low as 0.1°C [11].

Test lubricants and conditions
In all tests, a spherical roller was used with a radius of 9.525 mm in the rolling direction and 70.0 mm in the transverse direction. The roller was made from AISI 52100 bearing steel, with a R a roughness of < 11 nm. The Al 2 O 3 disc used had a diameter of 100 mm, thickness of 3.5 mm and roughness of 7 nm. For all tests the applied load was 20 N and the lubricant bath as maintained at 40°C. According the Hertz theory, the elliptical contact formed under these conditions has half widths of approximately 260 and 70 μm. Therefore, given the resolution of the camera, 84 × 14 temperature measurements were possible within the contact. The traction fluid Santotrac 50 was used as lubricant in all of the tests. Santotrac 50 is formulated synthetic fluid based on a bicyclohexylalkane and is believed also to contain a polymethacrylate viscosity modifier, the antiwear additive zinc dialkyldithiophosphate, and an oxidation inhibitor. This type of lubricant has been designed for use in traction drives, which require high shear stress to transmit power without excessive normal contact stress [12]. As a result of this high shear stress and thus friction, high contact temperatures are generated [14].

Numerical method
The freely available CFD package, Openfoam, was used to solve the EHL problem numerically. Details of the numerical method used in this paper have been described previously by Hartinger et al. [4]. Here the governing equations are repeated and an extension of the viscosity calculation is presented.
The approach is based on treating the lubricant as a fluid/vapor two-component mixture. For this mixture the continuity equation is: The vapour fraction is given by:    where ρ l,sat is the saturated liquid density and ρ v,sat is the vapour density at the cavitation pressure p sat to be specified by the user. The compressibility of the mixture is given by where ψ l and ψ v are the compressibilities of the liquid and vapour. The mixture density is The Dowson liquid pressure density is given by where ρ l,0 is the density at zero pressure and p is given in Pas [15]. The momentum equation reads The mixture dynamic viscosity is assumed to be: Subsequently the temperature equation reads A model for the dependency of viscosity on pressure and temperature was proposed by Ref. [16] and developed by Ref. [17]: The dependence of viscosity on shear rate (shear thinning) is described by the Eyring model with the Eyring stress τ 0 depending on temperature and pressure outlined by Evans and Johnson [14].
where τ 0p is the Eyring stress at zero pressure and T τ,ref , α τ,p is the influence coefficient of pressure α τ,p , and β t1 and β t2 are the temperature influence coefficients. To avoid division by zero, the shear-thinning part is only used if the shear-rate is larger than γṁ in = 10 −8 s −1 , leading to the final formulation of the non-Newtonian viscosity.      The elastic deflection is modelled following the Hertzian contact theory in two dimensions given by Ref. [18]: For two contacting bodies, E 1,2 denote the Young moduli and σ 1,2 the Poisson ratios. The reduced elastic modulus E r is defined as       The geometry of a deforming half-cylinder with radius R is defined by where ω ref is the deflection of the fixed reference point and h c0 is the undeformed film thickness at the origin. The temperature boundary condition for moving solids are based on the work of Carslaw and Jaeger [19] [20], assuming a moving point heat-source and a semi-infinite body. For a one-dimensional line contact the temperature of the solids is given by The applicability of this integral depends on the Peclet number, which is the ratio of advection to the rate of diffusion and is defined for thermal diffusion as =

Pe
Lu α T (20) where L is the characteristic length scale, u is the velocity and α T is the thermal diffusivity defined as α k ρC T p (21) where k is the thermal conductivity, ρ is the density and C p is the heat capacity at constant pressure. The characteristic scales for an EHL contact are the contact radius a, the velocity of the moving surface u s and the thermal diffusivity of the solid α T,s . According to Johnson [21], at large Peclet numbers, Pe > 5, the heat will diffuse only a short distance into the solid in the time taken for the surface to move through the heated zone. The heat flow will then be approximately perpendicular to the surface. For Peclet numbers lower than Pe = 5, equation (19) should not be used.    caused by the curvature of the surfaces. This is not an issue, since, in the central contact region, where the surfaces are approximately parallel.

Experimental vs. computational results
The computation domain is assumed to be a two-dimensional cylinder in contact with the disc. In order to obtain results that are comparable to the experimental three-dimensional roller, a similar boundary condition has to be derived. For this, the maximum of the Hertzian contact pressure, p max , is calculated and then a load for the line contact is chosen to give the same maximum pressure of 0.54 GPa. This load corresponds to 61700 N/m, (see Table 1).
The common numerical parameters for all cases are listed in Table 2.
The parameters describing the viscosity for Santotrac 50 are taken from Evans and Johnson [14]. The pressure coefficient Z and temperature coefficient b are determined by the GRG-algorithm of Microsoft Excel to fit the data of Evans and Johnson. Only a relevant subset of the original data is chosen to fit the coefficients, spanning temperatures 40, 60, 80, 100°C and pressures up to 0.6 GPa. The parameters for the Eyring stress are also fitted to the data of Evans and Johnson using the GRG-algorithm; see Fig. 4 and Fig. 5.
One important deficiency is the implementation of the temperature boundary condition of the solid after Carslaw-Jaeger. This is not valid for all conditions used since the Peclet number should be greater than five and a lower Peclet number results in underprediction of sideways conduction. An overview of Peclet numbers is given in Table 3. To summarize the validity of the thermal boundary condition: -SRR = 0.5: disc surface is valid for all cases -SRR = 0.5: ball surface valid for entrainment velocities greater 413 mm/s -SRR = 1.0: disc surface valid for all cases -SRR = 1.0: ball surface valid for entrainment velocities greater 579 mm/s In the following section, experimental results are compared to computational results for SRRs of 0.5 and 1.0 and 810 mm/s entrainment velocity. Fig. 6 shows the temperature rise above ambient at the surface of the roller and disc for SRR 1.0 and entrainment velocity of 810 mm/s. The other cases with SRR 1.0 show all very similar behaviour. For a detailed analysis of the results, the entrainment velocity 810 mm/s is chosen as the thermal boundary conditions are valid and the effects are more pronounced.

SRR = 100% and entrainment velocity = 810 mm/s
In Fig. 6, the temperature distribution of ball and surface from experiment and CFD is displayed. The experimental results are taken from the centreline where the highest temperatures occur. The fluid flow is from negative to positive x-values. The x-coordinate is set to zero for the peak-temperature of the disc for experiment and CFD. There is very good agreement on the peak temperatures and x = [-0.1, 0] mm for both surfaces. In the inlet region (x = [-0.2,-0.1] mm), the experimental results are 1-4°C higher than the numerical results. The reason for this difference is not known. Possible reasons could be the under-prediction of side-conduction of the Carslaw-Jaeger approach or self-heating of the roller. For x = [0, 0.1] mm the experimental results are up to 5°C higher than the CFD. Again, the reason for this difference in not known. For x > 0.1 mm the experimental results on the ball are not valid due to optical interference in the cavitated region described above.
It should be noted that the CFD and IR measurements were carried out independently and no adjustments were made to either sets of data. It is therefore surprising that there is such good agreement between modelling and experiment is found -especially considering the shortcomings associated with both techniques (discussed below) and the fact that the Santotrac rheological properties used were taken from tests by Evans and Johnson in 1986. However, as later is shown the agreement drops significantly when the underlying assumptions of the temperature boundary conditions are violated. Fig. 7 shows the film thickness, pressure and velocity distributions across the film-thickness. The x-coordinate is zero at the centre of rotation. The pressure is constant across the film-thickness, in agreement with Reynolds equation. However, a shear-band is developing in the centre of the film-thickness in the region of high pressure. Fig. 8 shows the corresponding shear-rate. Fig. 9 shows the temperature rise of up to 46°C in the centre of the film thickness. Fig. 10 shows the corresponding viscosity distribution resulting from the combined effects of pressure, temperature and shear-thinning. This shows higher viscosities occurring at high pressures and low temperatures and a band of lowviscosity caused partly by temperature but more by shear-thinning. Fig. 11 shows the pressure-and temperature-dependent Eyring-stress, τ 0 , which is important for determining the onset of shear-thinning in the centre. The reduction of the viscosity and the high shear-rate due to shear-thinning is likely to be overestimated, which will be argued later when looking at the comparison with temperature and friction measurements in all cases. This could be caused by a too low Eyring-stress or modelling assumptions for shear-thinning.

SRR = 50% and entrainment velocity = 810 mm/s
Figs. 12 to 17 show corresponding data sets at a SRR of 50% instead of 100%. The temperature distributions show similar pattern of behaviour as the 100% SRR case but lower temperatures. As can be seen in Fig. 12, the peak temperatures again agree well. However, the experiments show higher temperatures than numerical results by up to 4°C in the inlet and outlet regions. Contrary to the SRR = 100% case, there is only a mild shear-band developing, as displayed in Fig. 13, Fig. 15 and Fig. 16. The temperature rise at the centre of the film-thickness is 26°C, lower than at 100% SRR, see Fig. 14. The viscosity, shown in Fig. 15, is also not constant across the film-thickness. However, the differences are not as pronounced as at SSR = 100% (Fig. 17). Fig. 18 compares the maximum temperature rise of the disc and roller from experiments and CFD for SRR = 50%. The temperatures in the range of validity of the Carslaw-Jäger condition agree very well with differences in the range of 1-2°C and with the gradients of the curves in agreement. The largest difference of 2°C can be observed at the lowest entrainment velocity of 211 mm/s. Fig. 19 compares the maximum temperature rise of disc and ball from experiments and CFD for SRR = 100%. The temperatures in the range of validity of the Carslaw-Jäger condition agree reasonably well, with differences of up to 3°C. The gradients of the curves do not agree with experiments, showing higher temperature with higher speed, whereas the CFD results are levelling off. This may be attributed to an overestimation of the reduction of viscosity due to shear-thinning effect shown in the cases with SRR = 100% leading to the levelling of the gradient and also to an increase of peak temperature as heat generation is more concentrated in the centre. For the cases below the range of applicability CFD cannot predicted the hotter surface accurately and the discrepancies with the experiment results are up to 6°C. This large deviation at low entrainment speeds may be attributed to the Carslaw-Jäger boundary condition, which is further outside its application bounds for the roller surface where it shows the highest deviations.

Comparison of friction coefficient
In Fig. 20, friction results from CFD are compared to experimental measurements. The latter are obtained in two different ways, (i) directly from strain gauges and (ii) inferred from IR-measurements using moving heat source theory. At a SRR of 50%, the friction coefficient from CFD and strain gauges agree well with a maximum deviation of 3%, which is consistent with the very good agreement between the temperature rises. At a SRR of 100%, there is a systematic offset of around 3-7% between CFD and strain gauge measurements, but the gradient is largely preserved. This is in better agreement than one would expect from the difference in temperatures shown above. The offset is likely to be caused by an overestimation of shear-thinning. As for SRR = 50% the friction coefficient and temperature are in good agreement and the principal difference between these cases is the shearrate.

Conclusions
CFD is a powerful method of predicting EHL behavior. In this work, CFD results are compared to IR microscopy and strain gauge friction measurements obtained for the same contact conditions. Reasonably good agreement is seen between the three approaches. This lends weight to the predictions of lubricant propertiessuch as viscosity, temperature and strain ratethrough the thickness of the film. The prediction of shear-thinning however seems to be overestimated.
The surface temperatures measured in this study arise due to complex, non-Newtonian, behavior occurring within the contact whereas film thickness is determined primarily from conditions in the contact inlet. For this reason, surface temperature measurements provide a more effective/stringent validation of the CFD approach than film thickness ones.
A number of simplifications exist in both the experimental and the CFD technique. For instance, using the most modern infrared microscopy approach it is now possible to map the temperature of the oil within the film in three dimensions. On the modelling side, recent advances include fully discretized modelling for solid displacement and solid temperature and also more complex models for thermal conductivity. The purpose of this paper has been to use the basic versions of these techniques in order to provide an initial comparison. Three dimensional validation using state of the art techniques is the subject of ongoing study.