Numerical modeling of the thermo-elastic fracture using a hybrid phase-field method is presented. The model’s capability to simulate complex thermo-elastic crack propagation paths under the temperature and mechanical loads’ influence is demonstrated through a few example problems. The evolution of complex crack patterns in a ceramic plate subjected to quenching are simulated. Further, the behavior of functionally graded specimens with a spatial variation of material properties when subjected to thermal shock is investigated. A staggered solution scheme is used to solve the coupled system of multi-field partial differential equations using an open-source finite element software FEniCS. © 2021 Taylor & Francis Group, LLC.