This paper addresses the modeling of the iron ore direct reduction process, a process likely to reduce CO₂ emissions from the steel industry. The shaft furnace is divided into three sections (reduction, transition, and cooling), and the model is two-dimensional (cylindrical geometry for the upper sections and conical geometry for the lower one), to correctly describe the lateral gas feed and cooling gas outlet. This model relies on a detailed description of the main physical–chemical and thermal phenomena, using a multi-scale approach. The moving bed is assumed to be comprised of pellets of grains and crystallites. We also take into account eight heterogeneous and two homogeneous chemical reactions. The local mass, energy, and momentum balances are numerically solved, using the finite volume method. This model was successfully validated by simulating the shaft furnaces of two direct reduction plants of different capacities. The calculated results reveal the detailed interior behavior of the shaft furnace operation. Eight different zones can be distinguished, according to their predominant thermal and reaction characteristics. An important finding is the presence of a central zone of lesser temperature and conversion.