Hybrid boundary node method (hybrid BNM) is applied for thermal analysis of 3D composites. A new formulation is derived for inclusion-based composites, which reduces degrees of freedom (DOFs) compared with conventional multi-domain solver. A version of fast multipole method (FMM) is coupled with new formulation for large scale analysis. Numerical examples are presented to analyze thermal behavior of composites with many inclusions.