The Immersed Boundary (IB) method proposed by Uhlmann for Direct Numerical Simulation (DNS) of fluid flow through dense fluid-particle systems is extended to systems with interphase heat transport. A fixed Eulerian grid is employed to solve the momentum and energy equations by traditional computational fluid dynamics methods. Our numerical method treats the particulate phase by introducing momentum and heat source terms at the boundary of the solid particle, which represent the momentum and thermal interactions between fluid and particle. Forced convection heat transfer was simulated for a single sphere and an in-line array of 3 spheres to assess the accuracy of the present method. Non-isothermal flows past stationary random arrays of spheres are investigated to assess the capability of our simulation method for dense particulate systems. All results are in satisfactory agreement with reported experimental and numerical results. (c) 2013 Elsevier Ltd. All rights reserved.