This study addresses stability challenges in oil shale reservoirs during the in situ conversion process by developing a thermo-hydro-mechanical damage (THMD) coupling model. The THMD model integrates thermo-poroelasticity theory with a localized gradient damage approach, accounting for thermal expansion and pore pressure effects on stress evolution and avoiding mesh dependency issues present in conventional local damage models. To capture tensile–compressive asymmetry in geotechnical materials, an equivalent strain based on strain energy density is introduced, which regularizes the tensile component of the elastic strain energy density. Additionally, the model simulates the multi-layer wellbore structure and the dynamic heating and extraction processes, recreating the in situ environment. Validation through a comparison of numerical solutions with both experimental and analytical results confirms the accuracy and reliability of the proposed model. Wellbore stability analysis reveals that damage tends to propagate in the horizontal direction due to the disparity between horizontal and vertical in situ stresses, and the damaged area at a heating temperature of 600 °C is nearly three times that at a heating temperature of 400 °C. In addition, a cement sheath thickness of approximately 50 mm is recommended to optimize heat transfer efficiency and wellbore integrity to improve economic returns. Our study shows that high extraction pressure (−4 MPa) nearly doubles the reservoir’s damage area and increases subsidence from −3.6 cm to −6.5 cm within six months. These results demonstrate the model’s ability to guide improved extraction efficiency and mitigate environmental risks, offering valuable insights for optimizing in situ conversion strategies.