Water-lubricated rubber bearings (WLRBs) are widely used in the propulsion systems of ships and military crafts. Based on the mass conservation boundary condition, the elastohydrodynamic lubrication model of WLRBs is established after considering the elastic deformation of the rubber liner and solved by the finite difference method. An improved algorithm is proposed to track the mass conservation boundary. For the first time, the algorithm allocates independent finite difference grids for the film pressure and fluid fraction to overcome the instability issues and to ensure the flow conservation of each grid cell. The accuracy and stability of the algorithm are verified by the experimental data in the literature, comparison of film pressures calculated by two different finite difference schemes, and flowrate calculations of two different bearings. The lubrication characteristics under the mass conversation boundary condition are compared with classical boundary conditions. The results indicate that the film pressure distribution, velocity distribution, and size and shape of the cavitation zone vary greatly under different boundary conditions in the clearance divergence region. The bearing capacity, attitude angle, and friction force under the mass conservation boundary condition are the largest; those under the double Reynolds boundary condition are smaller; and those under the Reynolds boundary condition are the smallest.