Steady-state smooth surface hydrodynamic lubrications of a pocketed pad bearing, an angularly grooved thrust bearing, and a plain journal bearing are simulated with the mass-conservation model proposed by Payvar and Salant. Three different finite difference schemes, i.e., the harmonic mean scheme, arithmetic mean scheme, and middle point scheme, of the interfacial diffusion coefficients for the Poiseuille terms are investigated by using a uniform and nonuniform set of meshes. The research suggests that for the problems with continuous film thickness and pressure distributions, the results obtained with these numerical schemes generally well agree with those found in the literatures. However, if the film thickness is discontinuous while the pressure is continuous, there may be an obvious deviation. Compared with both the analytical solution and other two schemes, the harmonic mean scheme may overestimate or underestimate the pressure. In order to overcome this problem artificial nodes should be inserted along the wall of the bearings where discontinuous film thickness appears. Moreover, the computation efficiency of the three solvers, i.e., the direct solver, the line-by-line the tridiagonal matrix algorithm (TDMA) solver, and the global successive over-relaxation (SOR) solver, are investigated. The results indicate that the direct solver has the best computational efficiency for a small-scale lubrication problem (around 40 thousand nodes). TDMA solver is more robust and requires the least storage, but the SOR solver may work faster than TDMA solver for thrust bearing lubrication problems. Numerical simulations of a group of grooved thrust bearings were conducted for the cases of different outer and inner radii, groove depth and width, velocity, viscosity, and reference film thickness. A curve fitting formula has been obtained from the numerical results to express the correlation of load, maximum pressure, and friction of an angularly grooved thrust bearing in lubrication.