Hybrid Trefftz Voronoi finite elements on polygon meshes are developed to capture thermal response of heterogeneous materials with randomly dispersed inclusions/voids. In this approach, the temperature intra-element field within the polygon is approximated by the homogeneous solutions to the heat conduction governing equation which are also called the T-complete functions. Whereas an auxiliary temperature frame field is independently defined on the polygonal boundary. The continuity at the matrix–inclusion interface is guaranteed by exterior and interior T-complete functions and that across element boundaries is enforced by incorporating the known hybrid functional at the element level. Thus, only the hybrid functional should be established for the matrix region. In conjunction with the divergence theorem, the domain integral involved in the functional vanishes, which leads to the element stiffness equation including boundary integrals only. It is reported that there exists a linear relationship between the maximum number of T-complete functions and the number of Gauss points sampled on each element side. Numerical results demonstrate that the proposed methodology is remarkably more efficient than the conventional finite element method (ABAQUS) and consequently becomes popular in micromechanical and microthermal modeling of composite and porous materials.