###### Abstract
Traditional numerical methods such as the finite element method (FEM), the boundary element method (BEM), the meshless method, and the extended finite element method (XFEM) are widely used in the study of the thermal response of complex FGMs. For example, [PERSON] et al. [11] analyzed the functional gradient plate with cracks under thermal load based on the FEM and indicated that the spatial distribution of material components has a significant influence on the thermal stress intensity factor (TSIF). [PERSON] and [PERSON] [12] also used the FEM to study the thermal fracture of an FGM thick-walled cylinder with cracks under thermal shock load. [PERSON] and [PERSON] [13] performed a thermo-elastic analysis of FGM hollow cylinders with surface cracks based on the FEM. [PERSON] et al. [14] studied the temperature distribution of the FGM plate under thermo-mechanical coupling and considered the influence of gradient parameters using the meshless Galerkin method. [PERSON] et al. [15] used Galerkin the BEM to simulate the FEG heat conduction process, and the simulation results agree with the analytical solution. [PERSON] et al. [16] investigated the thermo-elasticity of a three-dimensional FGM hollow cylinder based on the BEM and analyzed the temperature and stress changes inside the material. [PERSON] et al. [17] studied the effect of circular holes on the heat conduction of FGM plates through the BEM and compared them with the results of various numerical methods. [PERSON] [18] studied the thermo-elasticity of FGMs based on the combination of the radial integration method and the BEM and applied it to analyze the dynamic fracture of materials. [PERSON] and [PERSON] [19] studied the influence of edge crack length on the TSIF in FGM plates under thermal load based on the element-free Galerkin method. [PERSON] and [PERSON] [20,21] used the meshless method to simulate the steady-state heat conduction and thermal stress distribution in a two-dimensional FGM plate. Based on aerodynamics, [PERSON] et al. [22] studied the effects of the spatial distribution of material components, porosity, and boundary temperature on the transient thermal response of FGMs. [PERSON] et al. [23] used the numerical manifold method (NMM) to study the steady-state heat transfer of FGMs and compared it with the existing analytical solution. They concluded that this method can simulate the steady-state heat conduction problem of FGMs. However, most of these numerical methods are limited to studying the stress intensity factor at the crack tip when dealing with the fracture of FGMs under thermal shock loading. These methods cannot consider the propagation and extension of the crack. The reason is that these numerical methods are mostly continuous, and it is difficult to simulate crack initiation, propagation, and interaction. To this end, we studied FGM thermal-mechanical coupling and thermal cracking behavior using a GPU parallel multiphysics finite-discrete element software, namely MultiFracS, developed by [PERSON] [24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53]. A heterogeneous material evaluation function module is added in the continuous-discontinuous thermo-mechanical coupling model to consider the variation of material parameters with space position. Thus, the heat transfer and thermo-mechanical coupling behavior of FGMs can be studied. It should be pointed out that the difference between the coupled thermo-mechanical model in this paper and in the literature [34] is that the heat conduction model in this paper can consider the hindering effect of cracks on heat transfer and can dynamically update the sharing relationship of triangular element nodes on both sides of the cracks during the crack propagation. However, in the heat conduction model in the literature [34], the adjacent triangular element always shares nodes during heat transfer calculation, Therefore, it is impossible to consider the hindering effect of cracks on heat transfer.