Nonmetallic inclusions such as sulfides and oxides are byproducts of the steel manufacturing process. For more than half a century, researchers have observed microstructural alterations around the inclusions commonly referred to as “butterfly wings.” This paper proposes a model to describe butterfly wing formation around nonmetallic inclusions. A 2D finite element model is developed to obtain the stress distribution in a domain subject to Hertzian loading with an embedded nonmetallic inclusion. It was found that mean stress due to surface traction has a significant effect on butterfly formation. Continuum damage mechanics (CDM) was used to investigate fatigue damage and replicate the observed butterfly wing formations. It is postulated that cyclic damage accumulation can be the reason for the microstructural changes in butterflies. A new damage evolution equation, which accounts for the effect of mean stresses, was introduced to capture the microstructural changes in the material. The proposed damage evolution law matches experimentally observed butterfly orientation, shape, and size successfully. The model is used to obtain S-N results for butterfly formation at different Hertzian load levels. The results corroborate well with the experimental data available in the open literature. The model is used to predict debonding at the inclusion/matrix interface and the most vulnerable regions for crack initiation on butterfly sides. The proposed model is capable of predicting the regions of interest in corroboration with experimental observations.