A Computational Fluid Dynamics Study on Shearing Mechanisms in Thermal Elastohydrodynamic Line Contacts

: A computational ﬂuid dynamics (CFD) model of the thermal elastohydrodynamically lubricated (EHL) line contact problem has been developed for the purpose of exploring the physical processes that occur inside a thin EHL ﬁlm subjected to shearing motion. The Navier–Stokes equations are solved by using the ﬁnite volume method (FVM) in a commercial CFD software, ANSYS Fluent. A set of user-deﬁned functions (UDF) are used for computing viscosity, density, heat source, temperature of moving surfaces and elastic deformation of the top roller according to well-established equations commonly used in the EHL theory. The cavitation problem is solved by taking into account multiphase mixture ﬂow. The model combinations of Houpert and Ree–Eyring and of Tait and Carreau were used for modeling the non-Newtonian behavior of Squalane and the results were compared. Both rheological models suggest the existence of shear-band and plug-ﬂow at high ﬂuid pressure. Due to the di ﬀ erences in viscosity at GPa-level pressure, the chosen model has substantial inﬂuence on the computed shear stress and temperature distributions in the high-pressure region. This shows the importance of using correct rheology information in the whole range of pressure, temperature, and shear strain rate.


Introduction
Lubricated machine elements designed for transferring high loads and high rotational speeds mainly operate in the elastohydrodynamic lubrication (EHL) regime. In a typical hard EHL contact, parts in relative motion are separated by a highly pressurized thin lubricant film. High pressure inside the film is generated due to a converging gap that is formed by the moving surfaces, operating conditions and lubricant properties. In such films, fluid pressure ranges from about 0.5 to 3 GPa. As a consequence of the high fluid pressure, lubricant viscosity increases several orders of magnitude, and the lubricant can reach a quasi-solid state. The tremendous increase in viscosity has a couple of implications. Firstly, being highly viscous helps the lubricant to stay in the contact zone and to separate the contacting surfaces. The quasi-solid lubricant film acts as a wedge between two moving surfaces and forces them to deform elastically. Secondly, by increasing viscosity more heat is generated in the film. The amount of generated heat becomes even larger when sliding is increased. As a result of the high heat, viscosity and consequently shear stress decreases. All these factors contribute to the friction in such contacts [1][2][3] and understanding its causes is one of the major current issues in the industry.
In experimental studies of EHL two main approaches are commonly encountered [4]. The first approach focuses on measuring EHL quantities by using well-controlled instruments (e.g., ball-on-disc and twin-disc tribometer) that replicate the type of contact in a machine element that is being studied. The second approach is based on shearing a highly pressurized fluid sample at high strain rates by using specially designed rheometers, capable of reaching high pressure levels. Over the years, researchers working in both fields have generated a significant amount of knowledge about EHL. In general, both sides agree on the effects of operating conditions and material properties on EHL film thickness. However, there is still an ongoing debate on rheological model that best describes the physical behavior of EHL films. In many cases, researchers dealing with EHL friction measurements use a Ree-Eyring based model [4] to describe their experimental data. On the other hand, researchers dealing with high pressure rheometry are advocating a Carreau based model [4]. Several articles that heated up the debate were published in 2014 and 2015 by Spikes and Jie [4,5], Bair and his co-authors [6], and it seems that there is still no consensus on this topic.
It should be acknowledged that besides the two main approaches described above, another experimental approach for studying EHL contacts has been developed. In 1989, Cann and Spikes [7] proposed a method based on measuring the temperature inside the contact and extracting the viscosity and shear stress based on the Carslaw-Jaeger model. In the following years, several papers dealing with the proposed method were published [8,9].
When it comes to numerical modeling of EHL, the most common approach is based on solving the Reynolds equation, which is derived from the Navier-Stokes equations by making assumptions that proved to be valid for the studies of thin fluid films. The two most important assumptions are that the thickness of the film is much smaller than its length and that fluid pressure is constant in gap height direction. As being very accurate, the Reynolds equation represents a standard equation in numerical modeling of EHL. However, the development of the computing power in the last decades has enabled a few researchers to obtain a full solution of EHL based on the Navier-Stokes equations.
In 2002, Almqvist and Larsson [10] performed the first thermal CFD simulation of EHL. They modeled Newtonian fluid behavior by using Roelands viscosity equation. They found that the decreasing effect of temperature on viscosity results in better numerical stability and they were able to reach a pressure of approximately 0.7 GPa but with a significant computational cost. In 2004, the same group of authors performed CFD analysis of EHL on the surface roughness scale while using Roelands and Ree-Eyring equations for viscosity modeling [11]. They discovered that the Ree-Eyring model leads to a more stable numerical solution due to having a decreasing effect on viscosity and a suppressing effect on the singularity in the momentum equation. In 2008, Hartinger et al. [12] found that their CFD solution of thermal EHL is in good agreement with the generalized Reynolds solution, for non-Newtonian fluid behavior and pressure of about 0.6 GPa. The authors relied on the Roelands viscosity model modified by Houpert and on the Ree-Eyring rheological model. In 2016, Lohner et al. [13] used commercial multiphysics software to develop a model of the EHL line contact problem based on the generalized Reynolds equation. They have also used the Ree-Eyring model and found that their results are in good correlation with previously published CFD results with pressures of up to approximately 0.6 GPa. In 2017, Hajishafiee et al. [14] performed a CFD analysis of thermal EHL and reached higher pressure levels over 3 GPa. They computed viscosity by using the model combination of Houpert and Ree-Eyring. For achieving numerical stability at high pressure, the authors assumed physically unrealistic elastic behavior of steel parts with reduced Young's modulus up to approximately 690 GPa. They concluded that high computational expense of CFD technique is justifiable in cases when assumptions in the Reynolds equation are not applicable.
What is common for all the introduced CFD studies of EHL is that they all rely on the Roelands/Houpert and Ree-Eyring viscosity and rheology models. However, in many theoretical and experimental studies of EHL the model combination of Tait and Carreau is also widely used. Furthermore, the former CFD studies also reveal that the Reynolds approach might not be sufficient when the assumptions in the Reynolds equation are not applicable. As the aim is to study the physical processes inside a thin EHL film subjected to shearing motion under high load for different rheological behaviors, the CFD approach is used in order to not be misled due to the assumptions of the Reynolds approach. Additionally, differences in results for the model combinations of Houpert and Ree-Eyring and of Tait and Carreau will be presented and as it will be shown, the differences are substantial.

Methods
In the first part of this section, governing equations of fluid dynamics are presented. The second part of the section describes the equations from the field of EHL. In the last part of the section, fluid domain geometry and the used numerical methods are presented.

Governing Equations
Fluids in motion are described by the equations of conservation of mass and momentum. To include thermal effects, the equation of conservation of energy was coupled together with the former two. For a multiphase flow and the mixture model, the governing equations read as follows [15].
The conservation of mass equation [15], widely known as a continuity equation, reads where ρ m is the mixture density and → ν m is the mass-averaged velocity, that read When gravitational and body force are neglected, the conservation of momentum equation [15], commonly known as the Navier-Stokes equation, reads as follows where = τ is the stress tensor and → ν dr,k is the drift velocity for secondary phase k, that read In Equation (5), µ m represents viscosity of the mixture described by the following expression The conservation of energy equation [15] reads where S T is the heat source term that for EHL conditions takes into account heat generated in the lubricant film by the viscous heating (Q shearh ) and by the compression work (−W) [16] as follows

Cavitation Modeling
When modeling thin lubricating films, one must resolve the cavitation problem that occurs at the outlet of the contact zone, at the position where the two contacting surfaces are starting to diverge from each other. As a consequence of cavitation, negative pressure occurs, which must be avoided when computing fluid force (integral of fluid pressure) in the force-balance equation that reads [16] w = x outlet The cavitation problem was solved by using the "full cavitation model" by Singhal et al. [17] that can be used if a fluid is modeled as a mixture of liquid and vapor phases. According to this model, fluid density ρ and vapor volume fraction α are a function of the vapor mass fraction f as follows where indices v and l denote vapor and liquid phase, respectively. Vapor mass fraction f was obtained from the solution of the transport equation coupled together with the continuity and the Navier-Stokes equations, as follows Here, terms R e and R c denote vapor generation and condensation rates, respectively, which are dependent on the vapor pressure [15,17] as follows As recommended by Hartinger et al. [12], a vapor pressure of 5000 Pa was used in the present study.

Lubricant Properties
For modeling EHL conditions, user-defined functions (UDFs) were used for computing density, viscosity, heat source (Equation (9)), temperature profile of the moving surfaces and elastic deformation of the top roller. The UDFs were written in computer language C.

Density Equations
Commonly, lubricant density in isothermal conditions is computed according to the Dowson-Higginson expression [18] that reads where A = 0.6·10 −9 m 2 /N and B = 1.7·10 −9 m 2 /N. For modeling thermal effects, Bos [19] included a linear temperature dependence in the Dowson-Higginson equation, as follows Due to being more physically relevant, Habchi [20] recommends using a density equation that is derived from the Tait equation of state, that reads The disadvantage of the equations that rely on the Tait equation of state is that they require specific fluid data which are not easily found in the literature. The effect of pressure on density was included in Equation (19) by using a modified version of the Tait equation of state for the volume relative to the volume at ambient pressure V 0 , as follows [21] where K 0 is the isothermal bulk modulus at zero pressure that reads The temperature effect on density was included in Equation (19) by using volume at ambient pressure relative to the ambient pressure volume at the reference temperature V R , as follows [21]

Viscosity Equations for Newtonian Fluid Behavior
For the Newtonian fluid behavior and isothermal conditions, viscosity was computed according to the Roelands formula [22] that reads where Z is Roelands pressure-viscosity index that can be computed by using pressure-viscosity coefficient α, as follows Z = α 5.1·10 −9 (ln(µ 0 ) + 9.67) .

Rheological Models
For modeling non-Newtonian fluid behavior two rheological models are used, the Ree-Eyring and the Carreau. The Ree-Eyring model reads [24] η Ree−Eyring p, T, where τ 0 is the Eyring stress that denotes the threshold where the Newtonian fluid behavior ends and the non-Newtonian fluid behavior starts. To avoid division by zero, instead of η 0 in Equation (29) Hartinger [25] suggests using the Houpert equation (Equation (25)). In addition, to avoid numerical error during simulation, Hartinger also recommends using the Ree-Eyring model only if strain rate is larger than an arbitrarily chosen minimum value of . γ = 1 × 10 −8 s −1 , otherwise he suggests using Houpert equation [25]. This can be written in the following fashion According to Björling et al. [21], the Carreau equation for the non-Newtonian viscosity reads where µ is the limiting low-shear viscosity that is described by the Vogel-like formula Here, ϕ represents the dimensionless viscosity scaling parameter that reads [21] It should be mentioned that the model combination of Tait and Carreau does not incorporate limiting shear stress behavior of liquids [21], unlike the model combination of Houpert and Ree-Eyring in which temperature effect mimics this phenomenon. To overcome this drawback of the model combination of Tait and Carreau, Björling et al. [21] suggests truncating shear stress if the limiting value of shear stress is reached, which is defined as follows where Λ parameter is derived from experiments.

Surface Temperature
Surface temperature, evaluated at location x and caused by heat flux q f acting at locationx, was computed according to Carslaw-Jaeger thermal boundary condition that for 1D line contact reads [24] where index s denotes the solid part. However, Equation (35) has the limit of applicability depending on the Peclet number that can be calculated by using formula [24] where L represents a characteristic length scale. For EHL contacts, the characteristic length scale is the Hertzian half width b. Finally, in order to use Equation (35), it is recommended that the Peclet number is larger than 5 [24]. In 2018, Hartinger [25] showed that for the EHL line contact problem, surface temperatures predicted by Carslaw-Jaeger equation closely matched experimental measurements.

The Film Thickness Equation
The elastic deformation of the top roller, evaluated at the location x and caused by the pressure p acting at locationx, is modeled by using the film thickness equation that reads [16] h where E is the reduced elastic modulus and R is the reduced contact radius that read as follows

Mesh Generation and Numerical Method
The EHL line contact is common in many machine elements, such as spur gears or cylindrical roller bearings. The 2D fluid domain geometry shown in Figure 1 was found to be appropriate for solving this kind of problem, since in this type of contact the pressure distribution in the third dimension is uniform. The geometry consists of a roller on a flat surface. Although a bottom flat surface was used instead of a curved surface, the effect of the bottom curvature was taken into account in computing deformation of the top roller through reduced radius of curvature R (Equation (39)). The reduced geometry was chosen for the purpose of direct comparison of the obtained results against the results of Srirattayawong [16] who used similar geometry and the same CFD software. ANSYS SpaceClaim was used for creation of the CAD model and high-quality mesh was generated in ANSYS ICEM CFD. The equations of mass conservation, momentum and energy were solved in ANSYS Fluent that was based on the FVM. FVM approximates the solution of the non-linear partial differential equations on discrete places on mesh geometry which are termed cells. Cell refers to a small surface (2D) or volume (3D) surrounding each mesh point. In order to obtain a solution at each cell center, flux values at the cell edges were used. In this way, the properties of the equations coming from the conservation laws (mass, momentum, energy) were mirrored, where the change of a quantity in the cell center was only due to fluxes across the cell edges. The UDFs, based on the equations presented in the first part of this section, were dynamically loaded with the Fluent solver. The UDFs for density, viscosity and heat source were updated in each iteration while the UDFs for surface temperature and elastic deformation of the top roller were updated at the beginning of each time step. The film thickness equation was implemented in the UDF for dynamic mesh. It was coupled with the force-balance equation (Equation (11)) in such way that the parameter h 0 was corrected at the beginning of each time step until the forces are balanced. As recommended by Srirattayawong [16], motion of the mesh nodes in the interior of the fluid domain was controlled by using spring-based smoothing method.
A pressure-based solver was employed for solving the conservation equations in a sequential and iterative manner [16]. In addition, the chosen algorithm for pressure-velocity coupling was PISO [16]. When it comes to spatial discretization, the Green-Gauss node-based method was chosen for resolving gradients [16]. When using a multiphase mixture flow model, PRESTO pressure interpolation scheme is recommended. For other quantities that are a product of the conservation equations, second-order upwind scheme was used, except for vapor, where only first-order upwind method was available. Finally, for solving time derivatives, second-order implicit method was used.

Results and Discussion
In the first part of the section, the results of the mesh verification test are presented, which prove that the results of the CFD model are independent of the mesh resolution. In the second part of the section, the results for the isothermal conditions are given. In the third part of the section, the thermal results at the pressure of about 1.0 GPa are presented by using Squalane as a lubricant and the model combinations of Houpert and Ree-Eyring and of Tait and Carreau for computing viscosity.

Mesh Verification Test
Solving the EHL problem by using a CFD approach based on the Navier-Stokes equations is computationally very expensive. In addition to a fine mesh, complex physics and a large number of very complex UDFs slow down the computation. In order to speed up the computation, the CFD model was modified for a parallel run on eight cores by using an Intel(R) Xeon(R) E5-2697 v3 CPU. It was found that an even larger number of cores cause communication overhead that can again slow down the computation, therefore eight cores were deemed to be an optimum.
In gap height direction, 10 cells are used [16]. In the horizontal direction, a finer mesh is used in the central region of the fluid domain where higher pressure is generated and all the processes of interest take place. Three mesh resolutions of the central region are tested, as shown in Figure 2 and Table 1.

Results and Discussion
In the first part of the section, the results of the mesh verification test are presented, which prove that the results of the CFD model are independent of the mesh resolution. In the second part of the section, the results for the isothermal conditions are given. In the third part of the section, the thermal results at the pressure of about 1.0 GPa are presented by using Squalane as a lubricant and the model combinations of Houpert and Ree-Eyring and of Tait and Carreau for computing viscosity.

Mesh Verification Test
Solving the EHL problem by using a CFD approach based on the Navier-Stokes equations is computationally very expensive. In addition to a fine mesh, complex physics and a large number of very complex UDFs slow down the computation. In order to speed up the computation, the CFD model was modified for a parallel run on eight cores by using an Intel(R) Xeon(R) E5-2697 v3 CPU. It was found that an even larger number of cores cause communication overhead that can again slow down the computation, therefore eight cores were deemed to be an optimum. In gap height direction, 10 cells are used [16]. In the horizontal direction, a finer mesh is used in the central region of the fluid domain where higher pressure is generated and all the processes of interest take place. Three mesh resolutions of the central region are tested, as shown in Figure 2 and Table 1.   Results presented in Table 1 show that the difference in maximum pressure between the coarsest and the finest mesh resolutions is only 0.02 MPa. The case with the finest mesh had 63% longer CPU time than the case with the coarsest mesh. For the purpose of saving time, the small difference in the results was found acceptable, and thus the coarsest mesh resolution shown in Table 1 was used for obtaining all results presented in this article.

Isothermal Conditions, Newtonain Fluid Behavior, Low Pressure
Input parameters including operating conditions and properties of the considered steel as solid and oil as lubricant are given in Table 2. Dowson-Higginson (Equation (17)) and Roelands (Equation (23)) equations are used for modeling the effect of pressure on density and viscosity, respectively.  Vapor phase parameters for Squalane are assumed to be the same as for oil. 2 Thermal conductivity is taken from the referenced paper by taking an average value from the thermal conductivity vs. pressure graph. 3 Eyring stress is derived from the shear stress vs. strain rate graph from the referenced paper. 4 Roelands pressure-viscosity index is calculated by using Equation (24) and pressure-viscosity coefficient of α = 18.1 GPa −1 [30].
The pressure distribution curve shown in Figure 3a is characterized by the pressure spike near the outlet of the contact zone. Another typical EHL feature is the minimum film thickness at the outlet of the contact zone which can be seen in the film thickness distribution curve (Figure 3b). These results are in good agreement with the results reported by Srirattayawong [16] who reported a maximum pressure of 0.474 GPa and minimum film thickness of 0.1869 µm for the lowest mesh resolution he considered. In the present study, for approximately the same mesh resolution, the maximum pressure is 0.477 GPa and the minimum film thickness is 0.191 µm. The small difference in the obtained results can be assigned to a small difference in mesh or a possible different value of vapor pressure, which was not reported in the referenced study. resolution he considered. In the present study, for approximately the same mesh resolution, the maximum pressure is 0.477 GPa and the minimum film thickness is 0.191 μm. The small difference in the obtained results can be assigned to a small difference in mesh or a possible different value of vapor pressure, which was not reported in the referenced study.
Contour plots for the isothermal conditions are given in Figure 4.
(a) (b)  Contour plots for the isothermal conditions are given in Figure 4.
Contour plots for the isothermal conditions are given in Figure 4. The pressure contours (Figure 4a) show that pressure varies only in the gap length direction while the pressure variation in the gap height direction is not present or it is insignificantly small. The velocity streamlines in Figure 4b show that the moving surfaces drag lubricant particles through the narrow passage between them. As the lubricant particles approach the outlet, due to pressure gradient and moving surfaces, their speed increases and they reach a maximum speed at the minimum film thickness location. Such speed increase is a necessary condition for maintaining continuity of the flow. In Figure 4c, a slight change of oil density in the high-pressure region can be noticed. Figure 4d shows that viscosity is uniformly distributed in gap height direction, which is expected since viscosity is only pressure dependent and pure rolling conditions are imposed. It is also expected to have maximum viscosity at the pressure spike location. Strain rate contours in Figure 4e show that the strain rate is increased in the regions where the localized slope of the velocity profile is the steepest. From shear stress contours shown in Figure 4f, it can be seen that the highest shear stress is in the regions where the lubricant is very viscous and strain rate is high, which is at the pressure spike location. Finally, on the basis of the oil and vapor volume fraction contours given in Figure 4g,h, we can clearly distinguish the regions where liquid and vapor phases are dominating the mixture flow. As expected, the vapor phase is dominant in the cavitation region at the outlet of the contact zone.

Thermal Conditions, Non-Newtonain Fluid Behavior, High Fluid Pressure
To reach higher fluid pressure for less computation time and better numerical stability, ceramics is used as material of the solid parts. Ceramics have more than two times higher Young's modulus of elasticity than steel and thus suffer less elastic deformation under the same applied load. This brings benefit in numerical stability since solution becomes more stable as the solid wall deforms less. Squalane is used as lubricant since its properties for the model combinations of Houpert and Ree-Eyring and of Tait and Carreau can be derived from the available literature (see Table 2). A load of 120 kN is applied on the top roller. The speed of the bottom surface has been set to 3.75 m/s while the top surface rotates with 125 rad/s (1.25 m/s) which gives an entrainment speed of 2.5 m/s and a slide-to-roll ratio (SRR) of 1 (50% of sliding). The Pecklet number for the slower moving surface is 14, and 42 for the faster moving surface which means that the Carslaw-Jaeger thermal boundary condition is appropriate for this case.
The pressure and the film thickness results are presented in Figure 5.  Figure 5a shows that the pressure distribution for the model combination of Houpert and Ree-Eyring is very smooth, while for the model combination of Tait and Carreau it is characterized by a very pronounced pressure spike. This implies that the effect of temperature and strain rate on fluid is much stronger in the case when the model combination of Houpert and Ree-Eyring is used. Note that it is not only the rheological behavior that is different for the two model combinations. The Newtonian viscosity behavior is also different, see section 2.3. Different pressure distributions result in slightly different film thickness distributions, as Figure 5b shows. Namely, the higher fluid pressure in the case when the model combination of Tait and Carreau is used means that the lubricant opposes to the applied load with more force, which results in the less narrow passage between the contacting surfaces.
Tait   Figure 5a shows that the pressure distribution for the model combination of Houpert and Ree-Eyring is very smooth, while for the model combination of Tait and Carreau it is characterized by a very pronounced pressure spike. This implies that the effect of temperature and strain rate on fluid is much stronger in the case when the model combination of Houpert and Ree-Eyring is used. Note that it is not only the rheological behavior that is different for the two model combinations. The Newtonian viscosity behavior is also different, see Section 2.3. Different pressure distributions result in slightly different film thickness distributions, as Figure 5b shows. Namely, the higher fluid pressure in the case when the model combination of Tait and Carreau is used means that the lubricant opposes to the applied load with more force, which results in the less narrow passage between the contacting surfaces.
Contour plots for the two model combinations are presented in Figure 6. From the pressure contours given in Figure 6a, it can be seen that for the same operating conditions, the model combination of Tait and Carreau results in about 0.11 GPa higher maximal pressure. The velocity streamlines for both model combinations (Figure 6b) show that fluid particles in the vicinity of the slower moving top wall are stuck and do not flow easily. This phenomenon is known as a plug-flow and most probably, it will become even more pronounced at higher pressure.
is much stronger in the case when the model combination of Houpert and Ree-Eyring is used. Note that it is not only the rheological behavior that is different for the two model combinations. The Newtonian viscosity behavior is also different, see section 2.3. Different pressure distributions result in slightly different film thickness distributions, as Figure 5b shows. Namely, the higher fluid pressure in the case when the model combination of Tait and Carreau is used means that the lubricant opposes to the applied load with more force, which results in the less narrow passage between the contacting surfaces.
Tait and Carreau Houpert and Ree-Eyring Contour plots for the two model combinations are presented in Figure 6. From the pressure contours given in Figure 6a, it can be seen that for the same operating conditions, the model combination of Tait and Carreau results in about 0.11 GPa higher maximal pressure. The velocity streamlines for both model combinations (Figure 6b) show that fluid particles in the vicinity of the slower moving top wall are stuck and do not flow easily. This phenomenon is known as a plug-flow The strain rate contours presented in Figure 6c indicate the existence of a shear-band for both model combinations. Additionally, it can be noticed that the value of the maximum strain rate is higher for the model combination of Tait and Carreau. The viscosity contours in Figure 6d demonstrate that the shear-band is more pronounced when the model combination of Tait and Carreau is used, which can be assigned to the more viscous fluid, higher temperature and higher strain rate. Figure 6e shows that shear stress is much higher in the Tait and Carreau based case. This is the consequence of the more viscous lubricant and higher strain rate. The higher shear stress results in higher friction force and thus in higher friction. The coefficient of friction was 0.0405 in the Tait and Carreau based case and 0.0176 in the Houpert and Ree-Eyring based case.
When temperature contours (Figure 6f) for both model combinations are compared, it can be seen that the heat source is closer to the slower moving top wall. This is expected behavior since heat needs more time to penetrate into the faster moving bottom wall. The location of the heat source governs the location of the shear-band. As it can be noticed from Figure 6c, the shear-band is in both cases closer to the slower moving top wall. Additionally, Figure 6f shows that the model combination of Tait and Carreau gives significantly higher temperature in the region around the pressure spike location. The big difference in temperature at this particular location is due to the more viscous fluid which results in more viscous heating. Furthermore, warmer fluid at the central part of the contact zone results in warmer fluid at the outlet zone and thus a difference in temperature between the two model combinations can be seen in the cavitation zone. The difference in fluid temperatures causes a difference in wall temperatures. As Figure 7 shows, the walls are much hotter in the Tait and Carreau based case. zone results in warmer fluid at the outlet zone and thus a difference in temperature between the two model combinations can be seen in the cavitation zone. The difference in fluid temperatures causes a difference in wall temperatures. As Figure 7 shows, the walls are much hotter in the Tait and Carreau based case. For the model combination of Tait and Carreau, the shear stress is truncated to the limiting shear stress defined by Equation (34). By dividing shear stress (Figure 6e, left) by limiting shear stress (Equation (34)), the contour plot presented in Figure 8 is obtained, which serves to show regions in the fluid domain where the limiting shear stress is reached. From the contact inlet to the beginning of the cavitation zone, the highest ⁄ ratio is approximately 0.85, which means that the limiting shear stress has not be reached in this area. A ⁄ ratio higher than 1 can be seen only in the cavitation zone where the pressure is close to zero or even negative. Since limiting shear stress has no physical meaning in the cavitation zone, it can be claimed that under given operating conditions, the limiting shear stress of Squalane has not been reached. For the model combination of Tait and Carreau, the shear stress is truncated to the limiting shear stress defined by Equation (34). By dividing shear stress τ (Figure 6e, left) by limiting shear stress τ L (Equation (34)), the contour plot presented in Figure 8 is obtained, which serves to show regions in the fluid domain where the limiting shear stress is reached. From the contact inlet to the beginning of the cavitation zone, the highest τ/τ L ratio is approximately 0.85, which means that the limiting shear stress has not be reached in this area. A τ/τ L ratio higher than 1 can be seen only in the cavitation zone where the pressure is close to zero or even negative. Since limiting shear stress has no physical meaning in the cavitation zone, it can be claimed that under given operating conditions, the limiting shear stress of Squalane has not been reached.
the fluid domain where the limiting shear stress is reached. From the contact inlet to the beginning of the cavitation zone, the highest ⁄ ratio is approximately 0.85, which means that the limiting shear stress has not be reached in this area. A ⁄ ratio higher than 1 can be seen only in the cavitation zone where the pressure is close to zero or even negative. Since limiting shear stress has no physical meaning in the cavitation zone, it can be claimed that under given operating conditions, the limiting shear stress of Squalane has not been reached. The change of the viscosity of Squalane with pressure, under zero strain rate and constant temperature (Newtonian behavior), is presented in Figure 9. The change of the viscosity of Squalane with strain rate, under constant pressure and constant temperature (non-Newtonian behavior), is shown in Figure 10. It should be mentioned that the results presented in these two figures were obtained without truncating shear stress to the limiting value, since Figure 8 shows that the limiting shear stress has not been reached in the CFD simulation based on the model combination of Tait and Carreau. Here, the viscosity is instead dominated by thermal softening. The change of the viscosity of Squalane with pressure, under zero strain rate and constant temperature (Newtonian behavior), is presented in Figure 9. The change of the viscosity of Squalane with strain rate, under constant pressure and constant temperature (non-Newtonian behavior), is shown in Figure 10. It should be mentioned that the results presented in these two figures were obtained without truncating shear stress to the limiting value, since Figure 8 shows that the limiting shear stress has not been reached in the CFD simulation based on the model combination of Tait and Carreau. Here, the viscosity is instead dominated by thermal softening. This is in agreement with the conclusion derived by Habchi et al. [31], who found that when operating conditions are such that the limiting shear stress is not reached, then truncating shear stress to the limiting value is of no importance for predicting friction.
As it can be seen from Figure 9, at zero strain rate and reference temperature of T = 313.15 K, there is only a slight difference in viscosity prediction between the model combinations Houpert and Ree-Eyring (H and R-E) and Tait and Carreau (T and C). However, if temperature is increased, it can be noticed that the difference in viscosity increases. At T = 353.15 K and p = 1.0 GPa the results differ by approximately one order of magnitude. This indicates that the temperature effect is mainly responsible for the difference in viscosity predictions by Tait and Houpert Newtonian viscosity models. This is in agreement with the conclusion derived by Habchi et al. [31], who found that when operating conditions are such that the limiting shear stress is not reached, then truncating shear stress to the limiting value is of no importance for predicting friction.
As it can be seen from Figure 9, at zero strain rate and reference temperature of T = 313.15 K, there is only a slight difference in viscosity prediction between the model combinations Houpert and Ree-Eyring (H and R-E) and Tait and Carreau (T and C). However, if temperature is increased, it can be noticed that the difference in viscosity increases. At T = 353.15 K and p = 1.0 GPa the results differ by approximately one order of magnitude. This indicates that the temperature effect is mainly responsible for the difference in viscosity predictions by Tait and Houpert Newtonian viscosity models.
. When the effect of strain rate on viscosity is taken into consideration (Figure 10), it can be noticed that at the lowest examined value of strain rate and pressure of up to approximately 0.2 GPa, there is only a small variation in viscosity between the two combination of models. However, at pressures above 0.2 GPa, the difference is becoming larger as pressure is increased, for the all three sets of conditions. At γ = 5e7 1/s and p = 1.0 GPa the results differ by approximately one order of magnitude.  When the effect of strain rate on viscosity is taken into consideration (Figure 10), it can be noticed that at the lowest examined value of strain rate and pressure of up to approximately 0.2 GPa, there is only a small variation in viscosity between the two combination of models. However, at pressures above 0.2 GPa, the difference is becoming larger as pressure is increased, for the all three sets of conditions. At . γ = 5e7 1/s and p = 1.0 GPa the results differ by approximately one order of magnitude. When the effect of strain rate on viscosity is taken into consideration (Figure 10), it can be noticed that at the lowest examined value of strain rate and pressure of up to approximately 0.2 GPa, there is only a small variation in viscosity between the two combination of models. However, at pressures above 0.2 GPa, the difference is becoming larger as pressure is increased, for the all three sets of conditions. At γ = 5e7 1/s and p = 1.0 GPa the results differ by approximately one order of magnitude.

Conclusions
The presented results in this article demonstrate that the CFD technique can be successfully used in numerical modeling of the EHL line contact problem. The CFD approach has the advantage of not relying on assumptions as the Reynolds approach. One of the biggest disadvantages of the CFD approach is a high computational cost, which is directly caused by the low under-relaxation factors required to control the motion of the solid parts and pressure generation inside a lubricant.
The obtained solution for the isothermal conditions at low pressure of about 0.5 GPa was found to be in a very good agreement with the recent literature, which speaks in favor of the validity of the conducted numerical procedure.
The obtained results for the thermal conditions at high pressure of about 1 GPa show that the solution depends highly on the considered rheological model. This indicates that the choice of a rheological model and/or lubricant properties, at higher pressure levels, can cause significant deviation between numerical simulations and experimental measurements. This is in agreement with the recent study by Hartinger [25], who concluded that at the pressure of about 0.6 GPa, his CFD model of the EHL line contact problem overestimates shear-thinning effect on viscosity when compared to experimental measurements of friction and surface temperatures.
From this study it is not possible to say if one of the two considered model combinations is better than the other, but shows that it is of significant importance to correctly determine the effective viscosity at high pressure and high shear strain rates. Otherwise, the occurrence of shear-bands and the friction cannot be predicted correctly.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature α
Pressure-viscosity coefficient Pa −1 α k Volume fraction of phase k β Temperature-viscosity coefficient