Abstract In this study, we used the finite volume method to computationally model natural convective flow in packed bed geometry. Using the OpenFOAM® v2112 code, we performed the computational analysis. We successfully meshed the intricate packed bed flow geometry, which consists of several spheres positioned at random. The spheres have sizes of 0.006 and 0.01 m, and the associated Rayleigh numbers are 1.83 × 107 and 8.48 × 107 respectively. We used the packed bed heights of H/d = 5, 10, and 20 in the simulations. By comparing the results of the OpenFOAM® v2112 simulations of the natural convection flow for all self-heating sphere in a packed bed, we demonstrated that the velocity distributions and Nusselt values are in good agreement with the experimental data. Additionally, it was evident from the velocity and temperature distributions in a packed bed core that there was a major temperature rise at nearby low velocity fields and a minor velocity rise in the intermediate and upper elevations. We showed that increasing the height of the pebble-bed core and correspondingly increasing the quantity of spheres inside it makes the flow more difficult and also generates local hot spots. This study is notable for using the finite volume method to evaluate natural convection flow in all self-heating packed beds and for simulating packed bed flow using a significant number of spheres. These two factors contribute to the originality of this work.