As a typical multibody system, rolling bearing consists of several lubricated frictional pairs which transmit motions and forces among components. In this paper, a multibody lubrication dynamics (MLD) model of cylindrical roller bearings (CRBs) is developed by coupling the bearings dynamics sub-model and several lubrication sub-models of each frictional pairs in real-time. The lubrication sub-models are discretized by using finite difference and solved by numerical iteration methods of Gauss-Seidel with line relaxation and discrete convolution-fast Fourier transform (DC-FFT) technique. What’s more, a parallel algorithm is also designed to enhance the calculation efficiency. Its results are validated by the experimental results in literature, and further compared by the results from other dynamics models. It is revealed that the film formation changes not only friction but also geometrical constraints and pressure distribution, which may introduce additional load inside bearing and increase rolling friction under high-speed conditions. As the dynamic oil film stiffness and damper are took into account by MLD, the cage has a lower whirl frequency and forces. The curve-fitting formulas of EHL film thickness are overestimates the thickness value, and thus change the load distribution and friction between roller and raceway.