A Thermal Equilibrium Analysis of Line Contact Hydrodynamic Lubrication Considering the Influences of Reynolds Number, Load and Temperature

Thermal effects such as conduction, convection and viscous dissipation are important to lubrication performance, and they vary with the friction conditions. These variations have caused some inconsistencies in the conclusions of different researchers regarding the relative contributions of these thermal effects. To reveal the relationship between the contributions of the thermal effects and the friction conditions, a steady-state THD analysis model was presented. The results indicate that the contribution of each thermal effect sharply varies with the Reynolds number and temperature. Convective effect could be dominant under certain conditions. Additionally, the accuracy of some simplified methods of thermo-hydrodynamic analysis is further discussed.


Introduction
Hydrodynamic lubrication has been used extensively to reduce friction loss between metal surfaces and the corresponding wear. In 1886, Osborne Reynolds proposed the famous Reynolds equation, which marked the beginning of theoretical hydrodynamic lubrication analysis. Gradually, some researchers [1,2] recognized the importance of including thermal effects in hydrodynamic lubrication analysis. In recent years, to achieve the higher lubrication and cooling demands of engineering design, research on thermal effects has become meaningful and important because of its great influence on lubrication performance [3][4][5][6].
In 1962, Dowson [7] improved the previous studies and proposed the thermo-hydrodynamic (THD) theory, in which variations of temperature, viscosity and heat transfer within an oil film were considered. Through further studies, Dowson et al. [8] and some other researchers [9,10] found that the steady-state THD condition includes three forms of thermal effects: viscous dissipation produced by the shear-rate within the oil film; conduction from or to the metal surfaces; and convection within the lubricant. When considering the unsteady-state THD conditions, the change of internal energy of the lubricant should be taken into account.
During the last few decades, numerous studies have reported on the factors that influence the thermal effects in THD lubrication. Pierre et al. [11] studied the influences of geometric parameters and operating conditions on the THD lubrication performance of journal bearings and showed that oil film thickness, friction loss and temperature distribution are affected significantly by the bearings' diameter, load and rotational speed; and also showed that the viscous dissipation could not be neglected in THD analysis. Chauhan et al. [12] concluded that load carrying capacity, temperature and pressure distribution are highly dependent on viscosity. Wang et al. [13] derived modified Reynolds and energy equations to investigate both thermal and cavitation effects in a journal bearing using micropolar fluids as the lubricant. Chun et al. [14,15] compared numerically the different effects of real and constant density and specific heat parameters values on high-speed journal bearing THD analysis, and concluded that the influences of real density and specific heat values on the calculation of the temperature, load and frictional loss could not be ignored. Boncompain et al. [16] and Hatakenaka et al. [17] studied the influences of reverse flow and showed that the temperature gradient across and along the oil film has a significant effect on the lubrication behavior. Harigaya et al. [18,19] studied the effects of viscosity and temperature on piston ring-cylinder liner lubrication of a diesel engine. The results showed that the oil film thickness increases with engine speed. An unreasonable estimation of viscosity and density would lead to a calculation error that cannot be ignored. Meanwhile, the calculated film thickness varies due to the temperatures of the ring and liner surfaces.
However, because there are numerous factors that influence the thermal effects and because the physical objects studied differ greatly, the previous reports are inconsistent regarding the relative contributions of viscous dissipation, conduction and convection to hydrodynamic thermal effects. Boncompain et al. [16] concluded that most of the heat generated was removed by convection. In contrast, Chun et al. [14] held the opinion that heat convection played only a small role in determining the friction loss and load and that convection could even be ignored in some cases. Conversely, Harigaya et al. [18,19] showed that degrees of heat transfer due to conduction and convection vary significantly with the boundary conditions and that the amounts of conduction and convection were similar, which meant that neither should be ignored in THD lubrication analysis. In reference [20][21][22], the mechanisms of heat transfer in pad trust bearings have been studied in detail, and the results in reference [20] show that heat dissipation distribution varies remarkably from load and speed in pad trust bearings, which might be helpful to heat transfer analysis in journal bearings. To address the above disagreements, the variations of the proportions of dissipation, conduction and convection thermal effects with changes in the influencing factors such as temperature, load and relative speed of the metal surfaces are worthy of further investigation.
In this paper, the thermal equilibrium of line contact THD lubrication is investigated by a combination of simulation and experiment. A steady-state THD analysis model of line contact lubrication has been developed to study the influences of relative speed (Reynolds number) of motion, load and the temperatures of the oil inlet and of the metal surfaces on the thermal effects, especially the proportions of each thermal effect. In parallel, a test rig corresponding to the simulation model has also been established to measure the metal surfaces' temperatures, the viscous heating rate and the minimum oil film thickness. The experiment data verified the credibility of simulation model.

Simulation Model
As shown in Fig 1, the simplified model has been used to study the thermal equilibrium of line contact hydrodynamic lubrication. The cylindrical surface can rotate around its geometric center, and its radius is 0.15 m (>>oil film thickness). The stationary block is fixed in the x axis and y axis directions, and it is pressed by a force F onto the rotating cylindrical surface in the z axis direction (along the oil film thickness direction). The oil is supplied from the left side. Because of the wedge-shaped structure and the relative motion of the metal surfaces, a stable hydrodynamic lubrication film can be formed. The linear speed of the cylindrical surface varies from 4.71 to 21.60 m/s, and the inlet oil temperature varies from 60 to 90°C. The block surface temperature varies from 60 to 100°C. The applied load F varies from 220 to 380 N. The nomenclature of symbols are listed in Table 1.
Generalized Reynolds equation. The Reynolds equation to describe the steady-state line contact lubrication behavior may be expressed as [15,23] @ @x where δ and U represent the oil film thickness and linear speed of the cylindrical surface. The viscosity η is allowed to be temperature dependent. The exponential model [19] is assumed: where a 0 = 7.36×10 −5 Pa s, T 1 = 1103.11°C and T 2 = 113.48°C. Other data of oil 15W40 are shown in Table 2.
The oil film thickness δ in this paper can be defined by the expression where δ 0 , B and R f represent the minimum oil film thickness, width of the block and radius of the flywheel.
The Half Sommerfeld Boundary Condition are used to solve the Reynolds equation: p 0 represents the atmospheric pressure. Force balance equation. Because stable oil film lubrication has been established, the integral value of the pressure distribution should be equal to the applied load F: where L represents the length of block. This equality represents the applied load of the line contact THD lubrication. CFD control equations. The mass conservation equation when density is constant:\ r Á @u @x þ @v @y þ @w @z ¼ 0 ð6Þ The velocity distribution is calculated by using the three-dimensional Navier-Stokes equations: Pressure inlet and outlet BC are used at the inlet and the outlet, and the pressure values are the same as Eq (4).
Considering the viscous dissipation heating effect in the oil film, the three-dimensional steady-state energy equation is expressed as follows: where F represents the dissipation term. The expression can be written as follows: The local viscosity in Eqs (6,7,8 and 10) is estimated by using Eq (2). The thermal boundary conditions are as follows: Thermal equilibrium equations. The thermal balance of viscous dissipation Q vis , conduction Q b , Q f and convection Q conv in the oil film can be written as follows: where Because the load F is not very large and all of the calculations in this paper are within the range of laminar flow, the elastic deformation and turbulence effects are neglected. Furthermore, the presented simulation model includes the influences of reverse flow at the oil inlet and outlet on the THD lubrication performance.

Numerical Procedure
To calculate the oil film thickness δ and the pressure distribution within it, a program is written using MATLAB to solve the Reynolds equation by the finite difference method. The main stages in the numerical procedure are summarized as follows: 1. The minimum oil film thickness δ 0 was assumed, and the viscosity η was calculated according to the oil temperature T n (in the first calculation, T 0 = T inlet ). Then the two-dimensional pressure distribution in the oil film was determined using the program mentioned above.
2. The integral value of the pressure distribution was calculated and compared with the applied load F. If they were not equal, the value of δ 0 was modified, and step (1) was iterated until the calculation converged (residuals of F < 0.5N).
3. The value of δ 0 was used to determine the domain of the oil film thickness, and the threedimensional CFD grid was established.
4. The Navier-Stokes equations and energy equations were solved by the Fluent 13.0 solver.

5.
The calculated volume-weighted average temperature T n+1 are compared with the T n . If the difference between them are larger than 0.1°C, step (1)-(4) were repeated. Usually 4-5 iterations were needed to reach convergence.
6. The velocity and temperature distributions Q f , Q b , Q vis , and Q conv were calculated.
For example, the pressure, temperature and velocity distribution at the mid-plane(y = L/2) of oil film, with boundary conditions of U = 11.40 m/s, T inlet = 70°C, F = 300 N, and T b = 100°C, are shown in Fig 2. The scale in the z axis direction(thickness direction) was enlarged 20 times for display. In the contours of pressure, pressure in the diverging area is set to 0 by applying the Half Sommerfeld Boundary Condition. In the contours of temperature, the convective effect can be seen at left side of the minimum oil film area.

Grid Independence Analysis
In this paper, 10 structured grid schemes have been used to reach grid independency, and the details of these grids are shown in Table 3. With boundary conditions of U = 11.40 m/s,

Experimental apparatus
An experimental apparatus was developed to validate the credibility of the calculated results. The schematic of the test system are shown in Fig 4 and Fig 5. The flywheel was used to imitate the cylindrical surface's rotation. Four electric heating bars were buried in the block to adjust T b . A spring-rail mechanism was designed to restrict the block to slide only in the z axis direction and to generate the applied load F at the same time. The force sensor and the capacitance type sensor (Fig 5) were used to measure the applied load F and the oil film thickness δ 0 , respectively, and the parameters of these sensors are shown in Table 4.
The surface temperatures T b and T f were measured by k-type thermocouples. The first-layer thermocouples were placed 1.5 mm from the metal surfaces, and the second-layer thermocouples were 6 mm from those surfaces. The temperature signals on flywheel side were saved by an storage measurement system and uploaded into the computer after experiments. The storage measurement system can ensure the time-synchronization between the temperature signals of both sides. Then, T b and T f were acquired by an interpolation method. The shaft torques M and M 0 were also measured. Because the measurements were under steady-state conditions, the motor power consumption was considered to balance the viscous dissipation which could be written as follows: 15W40 oil was used as the test lubricant in this paper. The oil could be heated within the pump station, and the oil temperature was controlled at a stable value (the maximum temperature fluctuation 0.94°C).

Comparison between experimental data and numerical results
The experimental conditions, which were also applied in the calculation, are shown in Table 5. Fig 6 shows the comparison between the experimental data and the numerical results for δ 0 and Q vis . The trends and the values agree well. The difference between the experimental data and the numerical results is less than ±13%, except for #1, #4 and #7, which might be caused by vibration of the working flywheel at specific speeds. The experiment data verified the credibility of the simulation model, thus this simulation model was adopted in the following studies.

Effect of the Reynolds number
The Reynolds number is defined by the expression [24]: Where The cases shown in Table 5 were used to study the effect of the Reynolds number on the thermal effects of line contact THD lubrication. T f data in Tables 5, 6   In Fig 7, Q b and Q f decrease when the Re increases, whereas Q conv and Q vis increase remarkably. Under low Re conditions, most of the heat is transferred by conduction. For example, the proportions of Q b and Q f in the thermal effects are over 80% when Re = 17.7, and therefore conduction plays a primary role in heat transfer under low Re conditions, meanwhile Q conv plays a very small role (4.3%) in bring heat into the oil film. However, when the Re is increased to 111.8, the contribution of the convection heat transfer rate Q conv to the heat dissipation is more than 50%. The proportion of the viscous heating rate Q vis increases dramatically with the Re.
To explain these results, the variation of the minimum oil film thickness δ 0 is shown in Fig  6. The value of δ 0 increases with the linear speed U, which agrees with the results of Dowson et al. [8]. Because δ 0 increases, the temperature gradient in the z axis direction decreases near both sides of the metal surfaces leading to a decrease of Q b and Q f . Conversely, the oil shear rate between these two metal surfaces increases according to Eq 10, which results in a significant increase of the viscous heating rate Q vis . Because the oil flow rate also increases with the linear speed U, the convection heat transfer rate Q conv increases from 3.6% to 51.7%, which means convection becomes significant in heat dissipation.
According to the discussion above, we can conclude that convection and viscous heating play important roles in the thermal effects of line contact THD lubrication, and the influence is Thermal Analysis of Line Contact Hydrodynamic Lubrication greater for higher relative speed of the metal surfaces. The influence cannot be neglected during line contact THD lubrication analysis; otherwise, the calculation will include large errors and cannot reflect reality, especially under high Re conditions.

Effect of temperature
This section discusses the influences of the oil inlet temperature T inlet and the surface temperatures T b and T f on the line contact THD lubrication as studied using the cases shown in Table 6. The speed U and the applied load F remained constant, and T b was set to 60, 80 and 100°C, respectively.  Thermal Analysis of Line Contact Hydrodynamic Lubrication The results are shown in Fig 8. With increasing T inlet , the effect of Q conv in the thermal effects changes from dissipating heat to bringing heat into the oil film. The point of conversion depends on the difference of T inlet and T b : If T inlet >> T b , Q conv brings heat; If T inlet << T b , Q conv dissipates heat; however, if the difference between T inlet and T b is small, it is difficult to deduce the direction (in or out) of heat flow without careful calculations. The calculation results for Q b reflect the same situation. With increasing T inlet , the value of Q vis decreases slightly, but the proportion of Q vis in the thermal effects changes dramatically. These results can be explained as follows: According to the research of Harigaya et al. [19], the minimum oil film thickness δ 0 decreases with increasing T inlet due to decreasing oil viscosity as shown in Fig  9. However, the decrease of δ 0 results in an increase of the shear rate in the oil film. Combining these two factors, the total amount of Q vis fluctuates slightly. Meanwhile, other effects fluctuate rapidly with variation of δ 0 , resulting in the fluctuation of the proportion of Q vis with temperature. From Fig 8 it can be concluded that Q b , Q f , Q conv and Q vis are coupled. The quantity and the direction of heat flux near the metal surfaces are determined not only by the temperature difference between the metal surfaces but also by the inlet temperature. Thus, Q b , Q f and the flow directions should be calculated for specific cases, and the proportions of Q b , Q f , Q conv and Q vis vary greatly with the boundary conditions. In some cases in the literature [25][26][27], the convection heat transfer effect was neglected and replaced with a one-dimensional thermal resistance model in THD lubrication analysis, and the quantity and the direction of the heat flux were determined only by the difference between the temperatures of the metal surfaces. This type of simplified method may not represent the real physical processes, and the calculation results need to be further modified.

Effect of applied load
The dimensionless load is defined by the expression [28]: The cases for analysis of the influence of applied load are shown in Table 7, where the block surface temperature T b and oil inlet temperature T inlet remained constant, and U was set separately to 4.71, 19.20, and 21.60 m/s. Fig 10 shows the calculated results of Q b , Q f , Q conv and Q vis for various values of applied load F. The results show that the values of Q b and Q f increase, and their proportions of the thermal effects change slightly with increasing applied load F. For low linear speed U (4.71 m/s), Q vis does not vary obviously with increasing F, but at high linear speed, Q vis increases far more obviously with applied load F. These results can be explained through the calculated variation of oil film thickness, which is shown in Fig 11 and is consistent with the previous studies [29]. The minimum oil film thickness δ 0 decreases when the applied load F increases regardless of the cylindrical surface speed U, but the degree of variation of δ 0 under high speed conditions is much larger than under low speed conditions. This behavior of δ 0 leads to two consequences: 1) a higher temperature gradient near the metal surfaces and 2) a higher shear rate in the oil film, which causes Q b , Q f and Q vis to all increase. However, regarding the proportions, Q vis changes more than the others, so the proportions of Q b and Q f decrease slightly. The value and proportion of Q conv both increase, but not obviously, with the applied load F. The calculation results for different linear speed U conditions show the same tendencies, with slight differences in the magnitudes and proportions. The results indicate that the proportions of Q b , Q f , Q conv and Q vis in the thermal effects remain almost unchanged with variation of the applied load F, whereas the cylindrical surface speed U has a greater effect on these proportions. From the view of dimensionless load, the proportion of conduction decreases while the proportions of convection and viscous dissipation increases with the increasing of F Ã .

Conclusion
To address the disagreements on the proportions of conduction, convection and viscous dissipation effects at various friction conditions in former studies, a simplified steady-state THD analysis model of line contact lubrication was developed to analyze the influences of Reynolds number, applied load, oil inlet temperature and temperatures of the metal surfaces.The calculation results agree well with the experimental data. The following conclusions can be drawn: The convective heat transfer rate and viscous dissipation heating rate was found increase with the Re. The proportion of convection in the thermal effects of line contact lubrication could be more than 50% under high Re conditions. It proves that the effect of convection on high Re THD lubrication analysis is significant and should not be neglected.
In some situations wherein the temperatures of the inlet oil and the friction pair surfaces are quite close, the quantity and direction of conductive and convective heat flow should be determined through calculation with specific boundary conditions, instead of being estimated from the difference of the friction pair surface temperatures. The one-dimensional thermal resistance model can not reflect the real heat transfer process.
The applied load has little impacts on thermal effect. The viscous heating rate increases with applied load slightly. The conductive and convective heat transfer rate remains almost constant with varying applied load. From the view of dimensionless load, the proportion of conduction In further studies, the influence of the temporal variation of the relative speed and roughness of surfaces would be taken into account. Additionally, turbulence and cavitation effects could be included in the coming studies.