One significant factor influencing geothermal energy exploitation is the variation of the mechanical properties of rock in high temperature environments. Since rock is typically a heterogeneous granular material, thermal fracturing frequently occurs in the rock when the ambient temperature changes, which can greatly influence the geothermal energy exploitation. A numerical method based on the numerical manifold method (NMM) is developed in this study to simulate the thermo-elastic fracturing of rocklike granular materials. The Voronoi tessellation is incorporated into the pre-processor of NMM to represent the grain structure. A contact-based heat transfer model is developed to reflect heat interaction among grains. Based on the model, the transient thermal conduction algorithm for granular materials is established. To simulate the cohesion effects among grains and the fracturing process between grains, a damage-based contact fracture model is developed to improve the contact algorithm of NMM. In the developed numerical method, the heat interaction among grains as well as the heat transfer inside each solid grain are both simulated. Additionally, as damage evolution and fracturing at grain interfaces are also considered, the developed numerical method is applicable to simulate the geothermal-related thermal fracturing process.