Steady-state heat transfer problems in heterogeneous solid are simulated by developing an adaptive extended isogeometric analysis (XIGA) method based on locally refined non-uniforms rational B-splines (LR NURBS). In the XIGA, the LR NURBS, which have a simple local refinement algorithm and good description ability for complex geometries, are employed to represent the geometry and discretize the field variables; and some special enrichment functions are introduced into the approximation of temperature field, thus the computational mesh is independent of the material interfaces, which are described with the level set method. Similar to the approximation of temperature field, a temperature gradient recovery technique for heterogeneous media is proposed, and based on the Zienkiewicz-Zhu recovery technique a posteriori error estimator is defined to automatically identify the locally refined regions. The convergence and performance properties of the developed method are verified by using three numerical examples. The numerical results show that (1) The convergence speed of the adaptive local refinement is faster than that of the uniform global refinement; (2) The convergence rate of the high-order basis functions is faster than that of the low-order basis functions; and (3) The existing inclusions change the local distributions of the temperature, and the extreme values of the temperature gradients take place around the inclusion interfaces.