Thermal interface material (TIM) is pivotal for the heat dissipation between layers of high-density electronic packaging. The most widely used TIMs are particle-filled composite materials, in which highly conductive particulate fillers are added into the polymer matrix to promote heat conduction. The numerical simulation of heat transfer in the composites is essential for the design of TIMs; however, the widely used finite element method (FEM) requires large memory and presents limited computational time for the composites with dense particles. In this work, a numerical homogenization algorithm based on fast Fourier transform was adopted to estimate the thermal conductivity of composites with randomly dispersed particles in 3D space. The unit cell problem is solved by means of a polarization-based iterative scheme, which can accelerate the convergence procedure regardless of the contrast between various components. The algorithm shows good precision and requires dramatically reduced computation time and cost compared with FEM. Moreover, the effect of the particle volume fraction, interface thermal resistance between particles (R-PP), interface thermal resistance between particle and matrix (R-PM), and particle size have been estimated. It turns out that the effective conductivity of the particulate composites increases sharply at a critical filler volume fraction, after which it is sensitive to the variation of filler loading. We can observe that the effective thermal conductivity of the composites with low filler volume fraction is sensitive to R-PM, whereas the it is governed by R-PP for the composites with high filler content. The algorithm presents excellent efficiency and accuracy, showing potential for the future design of highly thermally conductive TIMs.