Particle-resolved direct numerical simulations of non-isothermal gas-solid flows have been performed and analyzed from microscopic to macroscopic scales. The numerical configuration consists in an assembly of random motionless spherical particles exchanging heat with the surrounding moving fluid throughout the solid surface. Numerical simulations have been carried out using a Lagrangian VOF approach based on fictitious domain framework and penalty methods. The entire numerical approach (numerical solution and post-processing) has first been validated on a single particle through academic test cases of heat transfer by pure diffusion and by forced convection for which analytical solution or empirical correlations are available from the literature. Then, it has been used for simulating gas-solid heat exchanges in dense regimes, fully resolving fluid velocity and temperature evolving within random arrays of fixed particles. Three Reynolds numbers and four solid volume fractions, for unity Prandtl number, have been investigated. Two Nusselt numbers based, respectively, on the fluid temperature and on the bulk (cup-mixing) temperature have been computed and analyzed. Numerical results revealed differences between the two Nusselt numbers for a selected operating point. This outcome shows the inadequacy of the Nusselt number based on the bulk temperature to accurately reproduce the heat transfer rate when an Eulerian-Eulerian approach is used. Finally, a connection between the ratio of the two Nusselt numbers and the fluctuating fluid velocity-temperature correlation in the mean flow direction is pointed out. Based on such a Nusselt number ratio, a model is proposed for it.