The Effective thermal conductivity (keff) is one of the key thermal properties for packed beds in the presence of a stagnant fluid. In this study, a thermal discrete element model (DEM) based on the original Cheng-Yu-Zulli analytical model for mono-sized packed beds has been improved and implemented especially for mixed beds of different particle sizes or materials. In order to perform the DEM simulation for packed beds, a thermal contact theory considering three heat transfer mechanisms (solid contact conduction, solid-fluid-solid conduction and radiation) was derived and applied in the network of Voronoi cells obtained by radical Voronoi tessellation of the relevant beds. The numerical model was validated through a comparison with experimental results already reported in literature and a good prediction for the effective thermal conductivity was obtained for both mono-sized and multi-sized packed beds in a wide range of solid-to-fluid conductivity ratio. The model also showed a good performance to study the heat flow distribution as well as the coupled thermo-mechanical behavior of packed beds.