<p>This study employs an accurate and stable thermo-mechanical smoothed particle hydrodynamics with phase field damage (TM-SPH-PFD) numerical model to investigate thermal shock fracture behavior in brittle solids. The Robin temperature boundary condition is imposed on the material boundary using the finite difference method. By integrating damage parameters and modifying governing equations, the proposed model is capable of simulating and predicting the progressive degradation of both mechanical integrity and heat transfer efficiency caused by thermal shock cracks. To ensure numerical stability, a multi-rate time integration scheme is adopted for coupling the temperature and mechanical fields. Firstly, the correctness of applying the Robin boundary condition is verified through analytical solution comparison. Subsequently, parametric studies are conducted on ceramics and rocks, examining the effects of temperature differentials, specimen dimensions, and convective heat transfer coefficients on crack patterns and temperature evolution. Finally, the failure process of a dry hot rock wellbore under combined thermal shock and mechanical stress is simulated, demonstrating the potential of the TM-SPH-PFD model for engineering applications.</p>