To establish a solid foundation for the computational fluid dynamics (CFD) of grease for the rolling bearing, this study incorporated the viscoelasticity into the grease CFD for the first time. The White-Metzner viscoelasticity model was implemented by user-defined functionalities in a commercial CFD program. The proposed method was examined under two different situations, i.e. the cone penetration test and the stringiness test. Two greases with different consistencies were studied to evaluate the applicability of the model. Significant improvements were achieved particularly for the grease with the lower consistency value compared to the existing viscosity-only simulation with the Carreau-Yasuda model. Furthermore, this study shed light on the mechanism of the characteristic grease flow by analysing the development of the stress in grease.