A computational fluid 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 film subjected to shearing motion. The Navier–Stokes equations are solved by using the finite volume method (FVM) in a commercial CFD software, ANSYS Fluent. A set of user-defined 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 flow. 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-flow at high fluid pressure. Due to the differences in viscosity at GPa-level pressure, the chosen model has substantial influence 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.