Butterfly wings are closely related to the premature failure of rolling element bearings. In this study, butterfly formation is investigated using the developed semi-analytical three-dimensional (3D) contact model incorporating inclusion and material property degradation. The 3D elastic field introduced by inhomogeneous inclusion is solved by using numerical approaches, which include the equivalent inclusion method (EIM) and the conjugate gradient method (CGM). The accumulation of fatigue damage surrounding inclusions is described using continuum damage mechanics. The coupling between the development of the damaged zone and the stress field is considered. The effects of the inclusion properties on the contact status and butterfly formation are discussed in detail. The model provides a potential method for quantifying material defects and fatigue behavior in terms of the deterioration of material properties.