We have developed a direct numerical simulation approach combined with the immersed boundary (DNS-IB) method for studying heat transfer in particulate flows. In this method, fluid velocity and temperature fields are obtained by solving the modified momentum and heat transfer equations, which are due to the presence of heated particles in the fluid; particles are tracked individually and their velocities and positions are solved based on the equations of linear and angular motions; particle temperature is assumed to be constant. The momentum and heat exchanges between a particle and the surrounding fluid at its surface are resolved using the immersed boundary method with the direct forcing scheme. The DNS-IB method has been used to study the heat transfer of 225 heated spheres in a fluidized bed. By exploring the rich data generated from the DNS-IB simulations, we are able to obtain statistically averaged fluid and particle velocities as well as the overall heat transfer rate in the fluidized bed. Good agreement between the current study and the one by Pan et al. (2002) is found for the hydrodynamic properties of the bed such as pressure gradients within the bed and the relationship between fluidization velocity and bed solid fraction. The particle-averaged Nusselt number is found to increase as the fluidization velocity increases and the bed height rises; particles at the entrance of the bed tend to have the maximum heat transfer rate because of the higher particle-fluid temperature gradients in this region; as the fluid moves upward in the bed, it gets warmer, which reduces particle-fluid temperature gradients and decreases the transfer rate of particles. Published by Elsevier B.V.