Temperature Variation at Solid-Fluid Interface of Thin Film Lubricated Contact Problems

Many calculating methods have been already developed for solving contact problems of parts such as gears, cams, and followers under fluid film lubrication conditions considering the temperature and pressure dependence. Similarly, the determination of the elasto-hydrodynamic pressure distribution the processes taking place in the lubricant and the contacting bodies, as well as in their environment, have to be dealt with simultaneously for the determination of the temperature field. A system of equation for the modelling of thermo-elastohydrodynamic lubrication between two contacting bodies containing hydrodynamic, thermodynamic, and strength problems is a highly non-linear system which becomes even more so if the temperature and pressure dependence of the material properties are considered. To solve this system, scientists started to use the finite element formulation in the 1960s and it was found to be a promising and reliable method. Earlier, the lubrication analysts used only the h-version finite element method (h-FEM) till 1991, when the first usage of the p-version finite element method (p-FEM) was published in the literature. In order to reduce the problem, in case of point or line contact, the contact bodies can be handled as semi-infinite ones. Following this simplification that had been successfully applied for the gap size determination, a substructure model was defined using analytical solution of the moving heat source. Instead of an iterative way between the solid and fluid problem, in this paper we present an efficient solution when thermal model for lubricant and surfaces were coupled and solved by a direct numerical method in one step.


Introduction
For rolling sliding parts such as cams and followers, gears, and bearings often operate under high loads, high speed, and high slip; the local or global temperature rise caused by the heat dissipation generated by the pressure distribution acting on the surfaces and the tangential stresses developing in the lubricant, respectively, may reach a level resulting in a non-negligible deformation of the surfaces as well as influencing the lubricant properties. Sternlicht et al. [1] (1961) investigated thermal effects in elasto-hydrodynamic lubricated (EHL) contacts. Cheng and Sternlicht [2] (1965) assumed that the viscosity is constant through the oil film and solved the thermal EHL model. Later, Cheng [3] (1965) obtained a refined solution in which the TEHL model got variable viscosity through the film thickness. Murch and Wilson [4] (1975) proved the significant effect of thermal conditions on minimum film thickness applying high rolling speeds.
Sadeghi [5] (1990) published a complete numerical solution for TEHD (thermal elasto-hydrodynamic) problem using different slip effects on the contact area with the same conclusions as Murch and Wilson.

Theoretical Description of the Elastohydrodynamic Lubrication
A general case of contacting surface pairs in the status of TEHL (thermal elastohydrodynamic lubrication) can be seen in Figure 1. The gap between the mating bodies is filled with lubricant because of to the relative velocity of the bodies. The developed hydrodynamic pressure is affected by the motion of the lubricant. The relative velocity difference between the mating surfaces results in tangential (shear) stress in the lubricant which cause the flow of the fluid material. Thus, at the proper kinematic condition of the contacting bodies and the pressure distribution in the lubricant can maintain balance with the clamping force on the contacting bodies and with this prevent a body-tobody contact. Pressure distribution acting on the surfaces may reach a level resulting in a nonnegligible deformation of the surfaces. The pressure and the tangential stresses developing in the lubricant generate heat dissipation which cause local temperature rise influencing the lubricant properties. If the circumstances of contact developing underneath thermal elasto-hydrodynamic conditions of lubrication are needed to be modelled, then hydrodynamic, thermodynamic, and solidstructural problems must be solved simultaneously with a modelling system which is highly nonlinear, even by itself, but also because of material properties. However, these three main areas may be separated clearly from each other in respect of their basic equations. For the contact problem of lubrication theory-due to its nature-it is convenient to employ in general a Gauss coordinate system with axis z perpendicular to the center contact surface.
Because of the nature of the contact problem of lubrication theory, it is generally advisable to use a Gauss coordinate system of which the z axis is perpendicular to the ideal center of the contact area. However, in the case of point or spot contact the center contact surface may be considered as a plane with its normal being parallel with the line of action of the force pressing the surfaces to each other, i.e., of the contact pressure. Consequently, it is most convenient to use for the investigation and description of the phenomenon an orthogonal coordinate system with its axis z being coincident with the line of action of the contact pressure developed by the compressive force. For the contact problem of lubrication theory-due to its nature-it is convenient to employ in general a Gauss coordinate system with axis z perpendicular to the center contact surface.
Because of the nature of the contact problem of lubrication theory, it is generally advisable to use a Gauss coordinate system of which the z axis is perpendicular to the ideal center of the contact area. However, in the case of point or spot contact the center contact surface may be considered as a plane with its normal being parallel with the line of action of the force pressing the surfaces to each other, i.e., of the contact pressure. Consequently, it is most convenient to use for the investigation and description of the phenomenon an orthogonal coordinate system with its axis z being coincident with the line of action of the contact pressure developed by the compressive force.

Contact Pressure and Film Shape
The generalized Reynolds equation of Dowson (1961) [14] is highly non-linear partial differential equation with variable viscosity and density extents through the film thickness which is used to calculate the contact pressure generated by fluid film lubrication.
The gap height can be calculated as a sum of the sizes of the original geometry and deformation of a half-space to which the displacement of a rigid surface has been added [15]: where ∆ rigid is the displacement of rigid surface, h g is the original gap height, the δ i deformation of the contacting solid bodies. The exact details of the generalized Reynolds and displacement field equations and calculation with them can be found in literature [15,16] in order to determine the pressure and flow field of the liquid lubricant.

Penalty Cavitation
During the description of lubrication problems both the boundary conditions defined hypothetically and the limitations applicable to the free edge forming the boundary of the actual lubrication zone which developed due to the cavitation, respectively, were encountered simultaneously.
In the cavitation zone the lubricant flows adhered to the surfaces and were broken up into strips. For this reason, the magnitude of the gas filling factor must be taken into account also for the thermodynamical equation.
At the hypothetically assumed edge the primary equation definable for the pressure field: can be identified on the one hand and the secondary boundary conditions defining the flow rate of lubricant exiting at the edge, on the other hand: where p a is the external pressure and q a the flow rate of lubricant exiting at the edge. However, the lubrication region with the hypothetically assumed Γ c edge does not coincide in a significant portion of the cases with the region filled with unbroken lubricant film to 100%. The location of edge (Γ c cav ) of the actual lubrication region is not known beforehand but can be determined on the basis of the properties applicable to it. In case the lubricant film breaks down due to the effect of cavitation and the sub-cavitation pressure is neglected then the cavitational boundary condition introduced by Swift-Stieber [17], commonly called the Swift-Stieber boundary condition in the literature, should be satisfied at the primary and secondary boundary conditions inside the cavitation region and at its boundary: where p c is the cavitational or saturation pressure.
While the satisfaction of the primary and secondary boundary conditions (2)-(3) does not entail any special difficulty, numerous problems are raised by the cavitational boundary condition (4). The difficulty in the resolution process is constituted primarily by the fact that the boundaries of the lubricating film broken down by cavitation are unknown. Consequently, the boundaries of the initially assumed region do not coincide with the real boundaries (Γ c - Figure 2) of the lubricating film. Thus, the relationships written previously for a homogeneous phase continuous lubricating film cannot be applied in the entirety of the presumed region bounded by Γ c . The Elrod-Adams algorithm is widely used for extending the Reynolds equation inside the cavitation zone which includes the θ = ρ/ρc fractional film content [18]. Kumar and Booker [19,20] published the method which is fitted for the FEM procedure. This method separates the density and pressure determination. Using of linear correlation between the density and the viscosity was also suggested by Kumar and Booker [19,20]: Instead of the separation of the variables for the contact and cavitation zone, penalty cavitation method has been proposed by Szávai [16] where the pressure dependent density in the cavitation zone is approximated by a high gradient slope in case of sufficiently low pressure (i.e., under saturation pressure). With these conditions the density can be described as follows [16]: where γ(p) is the penalty function which is γ(p) = c if p < pc otherwise 0, where c is a sufficiently high number.
The pressure dependency of the density exists both in the lubrication region and cavitation zone.
Thus, ρ* range contains the cavitation zone and the lubrication region as well. The volume fraction based on (6) can be obtained as [16]: Applying assumption of Kumar and Booker of the density and viscosity correlation the viscosity can be obtained as [16]:

Kinematic Properties of Lubricant in the Gap
The F(τeq) is a characteristic function of a specific type of a lubricant model and τeq is the equivalent shear stress in = ' ⋅·⋅ ' where σ' is the stress deviator tensor.
In the case of various types of lubricant models, the function F(τeq) can be equal to the forms below [17]: The Elrod-Adams algorithm is widely used for extending the Reynolds equation inside the cavitation zone which includes the θ = ρ/ ρ c fractional film content [18]. Kumar and Booker [19,20] published the method which is fitted for the FEM procedure. This method separates the density and pressure determination. Using of linear correlation between the density and the viscosity was also suggested by Kumar and Booker [19,20]: Instead of the separation of the variables for the contact and cavitation zone, penalty cavitation method has been proposed by Szávai [16] where the pressure dependent density in the cavitation zone is approximated by a high gradient slope in case of sufficiently low pressure (i.e., under saturation pressure). With these conditions the density can be described as follows [16]: where γ(p) is the penalty function which is γ(p) = c if p < p c otherwise 0, where c is a sufficiently high number. The pressure dependency of the density exists both in the lubrication region and cavitation zone. Thus, ρ* range contains the cavitation zone and the lubrication region as well. The volume fraction based on (6) can be obtained as [16]: Applying assumption of Kumar and Booker of the density and viscosity correlation the viscosity can be obtained as [16]:

Kinematic Properties of Lubricant in the Gap
The F(τ eq ) is a characteristic function of a specific type of a lubricant model and τ eq is the equivalent shear stress in τ eq = 1 2 σ · · · σ where σ is the stress deviator tensor. In the case of various types of lubricant models, the function F(τ eq ) can be equal to the forms below [17]: where τ E is the Eyring shear stress and τ L is the limit shear stress. τ L can be taking account as a linear function of pressure of which coefficients τ l0 and χ: Since boundary conditions for the velocity field on the surfaces-as shown in Figure 1-are: where according to [15], and The velocity component in z direction can be decomposed as [15]: where W i is rigid body motion of the contacting body in z direction, while the first member comes from movement of non-plane surface in x and y direction. The filling parameter taken into consideration the derivative of velocity along the gap and the stress in this direction may be written in the following form considering that the value of the filling parameter in the contact region is 1 and the pressure in the cavitation zone is negligibly low, respectively: The stress deviator tensor takes the simplified form below for Newtonian fluids:

Theoretical Description of Thermodynamical Part of the TEHD
The next equation expresses the conservation of energy stating that the change in total energy during unit time is equal to the sum [16]: where Ξ is the dissipation originating in the internal friction of the lubricant: ϑ-the temperature, and c v -specific heat referred to constant volume. D-strain rate σ'-deviator stress tensor If the thermal properties of the continuum are presumed to be constant, then: for an uncompressible liquid (∇· Similarly to the determination of the elasto-hydrodynamic pressure distribution the processes taking place in the lubricant and the contacting bodies as well as in their environment have to be dealt with simultaneously for the determination of the temperature field. The details of the processes taking place in the contacting bodies will not be discussed here, but the assumption employed according to which the temperature fields developing in the bodies can be determined analytically or numerically on the basis of the boundary conditions applicable to the bodies, primarily in respect of temperature, and secondarily in respect of heat transfer.
Because of the nature of the problem the following simplification may be employed similarly to the simplifications seen at the derivation of the Reynolds equation according to [17]: However, this does not mean a significant easement for the solution since the temperature variation along the z axis cannot be eliminated. For the calculation of the temperature field the temperature variation along the gap also has to be determined unless the approach assuming parabolic temperature distribution along the cross-section is employed. In this latter case it is sufficient to determine a typical Under the conditions of elasto-hydrodynamic lubrication-as also outlined before-the flow may be considered to be laminar as the gap size is very small, barely a few µm. Thus-with turbulence neglected-the terms averaged in time may be considered to constitute the parameters usable for describing the fluid motion.
The dissipation in the lubricant is: Taking into consideration that the derivatives along the gap of the u and v fields are significantly larger than the derivative of the w field along the contact plane the dissipation will be the following with negligible error [21]: Thus, the dissipation can be written in the following form for Newtonian fluids:

Thermal Boundary Conditions
At variance with the boundary condition for pressure boundary conditions have to be defined for the temperature field both at the contact surfaces (S 1 , S 2 ) and for the entire cross-section of the oil film at the boundary of the contact zone.
Both surfaces S 1 and S 2 can be divided into two sections (S iϑ , S iQ ) enabling thereby to specify a primary boundary condition for the temperature field, on the one hand, and a secondary boundary condition defining the convection of heat through the surfaces: where ϑ si (x,y) is the temperature of surface S iϑ while q si is the convection of heat through surface S iQ which are either pre-defined or determined during the solution of the associated thermal problem (relationship between the lubricant and the contacting body). Oil inlet and outlet boundaries have free temperature condition (q n = 0). In the case of pure sliding the contact region moves together with one of the surfaces or-worded in another way-this surface is stationary in relation to the contact region. In this case the heat exchange taking place through the surface stationary in relation to the contact zone can be neglected in most cases since it is more intensive by several orders of magnitude through the surface of the body moving in relation to the contact zone. At this time the adiabatic boundary condition proposed by Rohde and Oh [8] can also be applied to the surface moving with contact region: The adiabatic modelling of this surface makes problem solving significantly easier, for example in the case of the quasi-stationary line contact of bodies regarded as infinite semi-spaces as the temperature distribution of the body moving with the contact region cannot be determined analytically or only in a very complex manner beyond the initial transient phase.
While generally accepted boundary conditions (30) and (31) have been developed for contacting surfaces S 1 and S 2 during the evolution of this field of science different boundary conditions may be encountered at the inlet and outlet cross-section of the lubricating film. An applicable approximation could involve the consideration of the derivatives of the temperature of the lubricant films carried on the surfaces taken in parallel with the plane of the surface to be negligibly small and the free surface of the lubricant in contact with air to be adiabatic with good approximation as it was presented in [16] ∂ϑ(x, y, z) A boundary condition similar to this was seen also in previous calculations [17] as well since no other reliable and practically applicable boundary condition could be defined because of the problems associated with backflow and the determination of the thickness of the lubricant film carried on the surfaces and of the initial position of the contact region in connection with this.

Consideration of Cavitation
In the cavitation zone the lubricant flows adhering to the surfaces and broken up into strips. For this reason, the magnitude of the gas filling factor must be taken into account also for the thermodynamical equation. Thus, similarly to (6) and (8)-assuming that condition (4) is satisfied in the hydrodynamic problem-the dissipation according to (26)-(28) will be with a negligible error [16]: The coefficient of thermal conductivity varies as the function of density like the viscosity, thus [16]: but, at the same time, the coefficient of specific heat does not change because of its specific nature: On the basis of the above, the differential equation of thermal conductivity (23)-(24) takes the modified form shown below: If the thermal parameters of the continuum are considered to be constant: Since, according to the law of the conversion of mass: Taking (7) into consideration term div(u) in Equation (37) of thermal conductivity may be written in the following form: or in another form, expanding the substantial derivative: Let us multiply the above equation with θρ L in order to obtain symmetric matrixes in the finite element model: The above equation may be simplified as follows in some basic cases.
If the lubricant density is constant (ρ L = const, dρ L /dt = 0), then: In the case of stationary contact: At stationary contact with constant lubricant density: Returning to original Equation (42) ρ L = ρ c may be considered to be approximately constant in the cavitation zone. Above the lubrication region, however, the value of the filling parameter is θ = 1.
Thus, in the lubrication zone: while in the cavitation zone Let us notice that ρ L may vary only in the lubrication zone where θ = 1 but, at the same time, may vary only in the cavitation zone where, however, ρ L = ρ c with good approximation.
On this basis, the new thermal conductivity equation derived with the use of the filling parameter will be [16]: The following simplified differential equation is obtained with constant density and thermal parameters:

Temperature Variation of Contacting Bodies as the Result of Heat Sources Distributed on Their Surfaces
As presented in Section 3.1, determining the temperatures of the contacting bodies and their surfaces is a high priority task during the solution of the thermo-dynamical problem. As obvious also from boundary condition (29) the temperature of the lubricant adhering to the surfaces in the contact region is identical with the surface temperature and, furthermore, heat defined by (30) affects the bodies through the surfaces: Similarly to the determination of surface deformation, modelling the contacting bodies as semi-infinitive bodies is applicable frequently also for the calculation of the thermal part of thermo-elasto-hydrodynamic problems. If this can be done, then the relationships according to [13] elaborated for the case of heat sources moving in and infinite half-space (at a point in infinite space [13]) in the determination of the temperature distribution of the bodies: while for line loads [16]: where K 0 is a second type, zero order modified Bessel function. Let us introduce the symbol shown below: which will take the following form in the same sense for line contact with (53) taken into consideration [16]: If the above approximation cannot be permitted, then the substructure approach well-known in the finite element practice may be employed in which case the thermal problem of contacting bodies is handled as a substructure of tribology problem. This is less efficient but the combined determination of the temperature field applicable to both the structure and the contact region can also be carried out if the construction of the cause-and-effect matrix necessary for the sub-structural model cannot be solved or only difficultly.

Weak Integral form of the Thermodynamic Equation
The weak integral form of thermodynamic Equation (42) can be set up with the use of a weight function w Q [16]: on the basis of (30) equivalent to: if vector u is introduced as a finite element symbol from w Q = N ϑ and vector u, then: Let us introduce the following symbol [16]: Thus, the discretized thermal equation will be [16]: taking the form below with constant density and thermal parameters:

Coupling of the Temperature Fields Developing in the Contacting Bodies and the Lubricant
Not only the lubricant but also the contacting bodies have to be investigated for the solution of the thermodynamical problem if the boundary conditions according to (30) and (33) are not predefined but have to be determined in the course of solving the thermal problem for coupled (lubricant-contacting body) fields. Two possible routes are available for such coupling. The traditional technique applicable to sub-structures is employed in one of the solutions, retaining the equations set up for both the lubricant and the contacting body at the boundary. This route can be followed easily if the variation solution based on the weak integral of the energy equation is used also for the contacting bodies. According to the other possible route the equation applying to the contacting bodies at the edge is considered to be valid and the variation of the differential equation applying to the lubricant is set up only for the region inside the lubricant. This later possibility is advantageous mainly if an analytical solution is desired to be employed with regard to the temperature of the contacting bodies.
In the finite element method, the unknown fields are attempted to be defined as a linear combination ( j c j ϕ j ) of expediently chosen approximating functions (ϕ j ) and unknown constants (c j ) in the course of solutions based on variation theory. Constants c j are defined so as to satisfy the weak form of the original equation. Following the finite element practice, symbol N will be used for the column vector of the spatial shape functions. While the T is the column vector of the time dependent c j "constant" values.
The approximation of the temperature distribution inside the lubricant is [16]: In this case Legendre type elements were used of which shape functions can be grouped as shown below: • The nodal functions taking the value of 1 at one node of the lattice dividing the region and 0 at the others and varying linearly over the region:Ń • Compiled from these the vector of the approximation functions is obtained in the form below for a 3D element: Let us note that on S i surfaces only the surface node-, edge-, and side shape functions of the elements of N ϑ have got nonzero values. The shape functions have got three groups: the group of the shape function operating inside the lubricant and groups of those related to surfaces S 1 and S 2 .
T can also be grouped similarly: Let us introduce the following symbols: While the internal shape functions assume 0 value at the contact surfaces, their derivatives will not be equal to 0 at the same locations. On this basis (50) may be written in the following the ith surface at z = z Si : Consequently, the surface temperature per (51) in consideration of (69) is: Let us follow the same approximation method based on least squares: Let us find the minimum of (71) as: Thus [16]: The arrangement of this results in: Let the equations below apply: and, furthermore: while: On the basis of these, the equation originating from the thermal boundary conditions applicable to surfaces S 1 , S 2 may be written as: That is: These equations-together with the thermal equation per (62)-have to be solved by iteration or in parallel with the modification that there w Q = N ϑV since the conditions expressed by (78) are given at the S 1 , S 2 edges. If the above boundary condition per (78) holds only for the contact surface marked i then w Q = N ϑVSj . In case Equations (62) and (78) are solved in a single iteration cycle, the systems of equations to be solved remain symmetric. In this case only parameters T V or T V,Si are calculated from Equation (62)-depending on whether the above boundary condition model holds for both contact surfaces or only for one of them-and parameters T sj from Equation (78). If, choosing the parallel solution mode, both equations are desired to be solved in a single step, the thermal equation takes the form below.
Or, if adiabatic boundary condition is employed for surface S j : Studying the equations obtained if can be found that they are not symmetric thus their solution requires more resources but only one step.
At variation with surface deformations, here the temperature distribution applicable to the surface cannot be incorporated into the thermodynamic equation of the fluid through the use of a cause-and-effect matrix but appears as a direct boundary condition and is determined either by iteration or as a coupled problem in a single step.
Together with the discretized Reynolds equation [15] and with the discretized form of gap size constitute [15] and with the relationship applicable to the parameters of displacement caused by pressure [15] the discretized thermal equation a closely coupled system but the thermal and contact problems must be manageable also separately at the same time.
As it can be seen above, modified WLF formula is used to describe the pressure-temperature relationship of the viscosity. For this task, the specific values of coefficients of modified WLF formula in Equation (82) can be seen in Table 1. Table 1. Coefficients of the WLF formula describing the pressure and temperature dependence of lubricant viscosity [16]. There were 15 elements in the gap geometry divided along length of contact area. While only a single element was assumed along the thickness for the calculation of temperature.

WLF Coefficient
The elasto-hydrodynamic part of the problem was solved using of optimized damped Newton-Raphson method [16]. The thermodynamical part of the problem was solved by damped Newton-Raphson [16]. The calculated pressure distribution is shown in Figures 3-6. as well as the temperature distribution and the gap size.
Subsequent to the problem of pure rolling contact, the calculations were carried out with 1, 1.9, 2 sliding ratios. Heat generation is much higher in this case, thus the mode of modelling the temperature distribution along the gap strongly influences the results as demonstrated well in the results published by Wolff and his co-authors. In our solution the temperature distribution along the gap was taken fully into consideration by employing fourth degree approximation.
The change in pressure distribution in comparison to the case of pure rolling contact is clearly apparent although a lower effect can be seen in comparison to the pressure distribution calculated by Wolff et al. with the temperature varying along the gap.
The pressure distribution indicated develops at the gap shape and temperature distribution shown in Figures 3-6 demonstrating well the definite temperature rise near the pressure peak shifting towards the upper surface moving more slowly.
The sliding ratio is:

Discussion
The weak integral form of the Reynolds and energy equation have defined and used as a basis of the finite element method. For approximating the pressure and temperature field, Legendre polynomials have been applied as shape functions. Solving the TEHD problems with a p-FEM based procedure allowed us to replace the smooth mesh with a coarse one. On the p-FEM procedure of thermal part of the TEHD problem necessary improvements have been made. This developed method and its application possibility is introduced in this paper. The developed method works in case of Newtonian

Discussion
The weak integral form of the Reynolds and energy equation have defined and used as a basis of the finite element method. For approximating the pressure and temperature field, Legendre polynomials have been applied as shape functions. Solving the TEHD problems with a p-FEM based procedure allowed us to replace the smooth mesh with a coarse one. On the p-FEM procedure of thermal part of the TEHD problem necessary improvements have been made. This developed method and its application possibility is introduced in this paper. The developed method works in case of Newtonian and non-Newtonian fluid and for both stationery and time dependent case. The thermodynamic model of the fluid-film contact has got non-constant temperature condition inside the film.
As for additional development, it is shown that boundary condition can be redefined as all inlets and outlets can be adiabatic boundaries, therefore, they can be treated uniformly. For the contact surfaces, a substructure model has been defined using analytical solution of the moving heat source where the input temperature data obtained from p-FEM calculation. Instead of an iterative way between the solid and fluid problem, in this paper we present an efficient solution when the thermal model for lubricant and surfaces have been solved as a coupled analysis in one step. In order to handle the cavitation problem a kind of penalty method has been applied on volume fraction parameter that is useful for the p-version finite element method based solution of thermal problem since the continuity equation is valid not only in the contact but in the cavitation zone as well.
The temperature distribution of the fluid film was modelled using higher order FEM approximation. The temperature modelling procedure was coupled to film thickness and pressure calculation. The developed calculation procedure has got a good efficiency for solving TEHL problems. This developed p-FEM procedure can reduce the model size compared to the conventional h-FEM methods and it is also possible to use this process in case of rough surfaces or dynamic load. The calculation performed shows that our p-FEM numerical procedure is valid in case of TEHL computation and the solution is obtained in a stable way.