Combined effects of slip velocity and volume fraction of slip spheres on the heat transfer characteristics of multiple slip spheres are numerically investigated within the framework of a free surface cell model combined with a linear slip velocity along the surface of the slip spheres. The governing conservation equations of the mass, momentum, and energy are solved by a segregated approach using a simplified marker and cell algorithm implemented on a staggered grid arrangement in spherical coordinates. The convection and diffusion terms of conservation equations are discretized using quadratic upstream interpolation for convective kinematics and second‐order central differencing schemes, respectively. Prior to obtaining new results, this numerical solver is validated by comparison of present results with the existing literature values. Further new results are obtained for a range of conditions as; Reynolds number, Re: 0.1–200; Prandtl number, Pr: 1–100; volume fraction of slip spheres, Φ: 0.1–0.5 and slip parameter, λ: 0.01–100. The effects of these dimensionless parameters on isotherm contours and local and average Nusselt numbers are thoroughly delineated. Finally, a new empirical correlation for the average Nusselt number of multiple smooth slip spheres is proposed on the basis of present numerical results.