Thermal residual stress is identified as one of the major reasons of stress concentration in material's microstructures which initiates failure in the microstructure. Considering the details of microstructural features (inclusions shape, size, and distri-bution) can provide better understanding of the thermal residual stress developed in the materials due to the temperature change. In this paper, we have extended the self-consistent clustering analysis (SCA) method for efficient and accurate modeling of thermal residual stress for thermoelastic heterogeneous materials. The governing equations of the thermoelasticity has been implemented through a eigenstrain problem and solved using a Fast Fourier Transform (FFT) based solution scheme and later extended for SCA formulation. The thermoelastic formulation of SCA method has been verified for a single inclusion (homogeneous and inhomoge-neous) problem with the analytical solution of Eshelby and the FFT-based solution for multiple clusters. An example problem with multiple inhomogeneous inclusions has been solved for different temperatures and the residual stress developed has been analyzed. Results show that the thermoelastic formulation of SCA can have order of hundred times speed up compare to the traditional FFT-based solution scheme. The proposed methodology can be implemented in thermomechanical problems and provide efficient multiscale capabilities for prediction of thermal residual stress.