A three-dimensional immersed boundary method (IBM) is applied for the solution of the thermal interactions between spherical particles in a viscous Newtonian fluid. At first, the free convection of an isolated isothermal sphere immersed in a viscous fluid is analyzed as a function of the Grashof number. A new correlation for the heat transfer rate from a single sphere is obtained, which is valid in the ranges 0.5 ≤ Pr ≤ 200 and 0 ≤ Gr ≤ 500. Second, the free convection heat transfer rate from pairs of spheres (bispheres) and from small spherical clusters immersed in air (Pr = 0.72) is investigated using this numerical technique. For bispheres, their orientation and the thermal plume interactions within a range of interparticle distances may cause the enhancement of the heat transfer rate above the values observed for two isolated spheres. For the simple triangular particle clusters, where the particles are in contact, it was observed that the average heat transfer rate per sphere decreases with the increased number of spheres in the cluster.