Abstract: Rolling element bearings as one of the most important industrial components are widely used in all kinds of mechanical products. The main component of its material for the low carbon steel and the crystal phase is mainly composed of martensite and austenite and other cubic lattice crystals with obvious anisotropic characteristics. Since today’s bearings are often designed to work under higher speeds and higher loads, even if the bearings are correctly mounted and adequately lubricated, they will fail with fatigue cracks inside the bearings under long term harsh operating conditions. Under these conditions, the effects of inhomogeneity and anisotropy on fatigue life become more important. A coupling algorithm for fluid-solid-temperature multi-physical fields was established, from which the effect of polycrystalline anisotropic material on the thermal elastohydrodynamic lubrication (EHL) was investigated. The Reynolds equation and the stress balance equation were used as the coupling boundary, and the stress-strain relation should be satisfied inside the solid. In order to represent the heterogeneity of polycrystalline anisotropic material on a microscopic scale, in numerical simulation it was necessary to model the stress concentration zone on the sub-surface of the contact zone, which required a very dense mesh to describe the stiffness differences between local grain boundaries. The oil film pressure and solid deformation were solved simultaneously by the multigrid (MG) method, and the temperature fields in solids and oil were solved by the line-by-line scanning method. The results showed that the MG method applied in this study was well suited for solving lubricated contact problem considering polycrystalline anisotropic material, which meant that the powerful solver could easily cope with millions of unknowns in the system of equations. For isothermal EHL condition, the pressure and film thickness by the coupling algorithm were consistent with the results obtained by the traditional EHL solver, so that the calculation accuracy was validated. In contrast to homogeneous isotropic material, anisotropic material grained with different crystal orientation angles could lead to local stress concentrations at grain boundaries. After coupling the temperature field, the stress field underwent an obvious counterclockwise deflection by the action of the oil film traction force. The anisotropy of the solid grains caused violent fluctuations in the pressure and temperature rise profiles, and produced a concentration of the stress at the grain boundaries, while the influence on the film thickness profile was very small. The establishment of this algorithm helped to understand the relation between the microstructure of the material and the macroscopic lubrication contact. It was further demonstrated that the MG method could achieve the required performance for accurate simulation of the sub-surface stress field, which could be applied to the sub-surface analysis and optimization of bearing steel material as well as the computational diagnosis of lubrication contact processes. It was of great significance for studying the influence of the microstructure of bearing steel on fatigue limits and the tribological behavior of complex interface systems under service conditions.