A Direct Numerical Simulation (DNS) method has been developed to solve the heat transfer equations for the computation of thermal convection in particulate flows. This numerical method makes use of a finite difference method in combination with the Immersed Boundary (IB) method for treating the particulate phase. A regular Eulerian grid is used to solve the modified momentum and energy equations for the entire now region simultaneously. In the region that is occupied by the solid particles, a second particle-based Lagrangian grid is used, which tracks particles, and a force density function or an energy density function is introduced to represent the momentum interaction or thermal interaction between particle and fluid. The numerical methods developed in this paper have been validated extensively by comparing the present Simulation results with those obtained by others. (C) 2008 Elsevier Ltd. All rights reserved.