Large bearings employed in wind turbine applications have half-contact widths that are usually greater than 1 mm. Previous numerical models developed to investigate rolling contact fatigue (RCF) require significant computational effort to study large rolling contacts. This work presents a new computationally efficient approach to investigate RCF life scatter and spall formation in large bearings. The modeling approach incorporates damage mechanics constitutive relations in the finite element (FE) model to capture fatigue damage. It utilizes Voronoi tessellation to account for variability occurring due to the randomness in the material microstructure. However, to make the model computationally efficient, a Delaunay triangle mesh was used in the FE model to compute stresses during a rolling contact pass. The stresses were then mapped onto the Voronoi domain to evaluate the fatigue damage that leads to the formation of surface spall. The Delaunay triangle mesh was dynamically refined around the damaged elements to capture the stress concentration accurately. The new approach was validated against previous numerical model for small rolling contacts. The scatter in the RCF lives and the progression of fatigue spalling for large bearings obtained from the model show good agreement with experimental results available in the open literature. The ratio of L10 lives for different sized bearings computed from the model correlates well with the formula derived from the basic life rating for radial roller bearing as per ISO 281. The model was then extended to study the effect of initial internal voids on RCF life. It was found that for the same initial void density, the L10 life decreases with the increase in the bearing size.