Detailed Modeling of the Direct Reduction of Iron Ore in a Shaft Furnace

This paper addresses the modeling of the iron ore direct reduction process, a process likely to reduce CO2 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.


Introduction
The direct reduction (DR) of iron ore, usually followed by electric arc steelmaking, is an alternative route to the standard, blast furnace, basic oxygen route for making steel. Annual DR iron production (86 Mt in 2017) remains small, compared to the production of 1180 Mt of blast furnace pig iron [1]. However, an attractive feature of DR, compared to blast furnace reduction, is its considerably lower CO 2 emissions, which are 40 to 60% lower for the DR-electric arc furnace route, compared to the blast furnace, basic oxygen route [2]. Among DR processes, shaft furnaces represent over 82% of the world's DR iron production, with the two main processes being MIDREX (65%), as shown in Figure 1, and HYL-ENERGIRON (17%) [3].
In a DR shaft furnace, a charge of pelletized or lump iron ore is loaded into the top of the furnace and is allowed to descend, by gravity, through a reducing gas. The reducing gas, comprised of hydrogen and carbon monoxide (syngas), and obtained by the catalytic reforming of natural gas, flows upwards, through the ore bed. Reduction of the iron oxides occurs in the upper section of the furnace, at temperatures up to 950 • C. A transition section is found below the reduction section; this section is of sufficient length to separate the reduction section from the cooling section, allowing an independent control of both sections. The solid product, called direct reduced iron (DRI) or reduced sponge iron, is cooled in the lower part of the furnace, down to approximately 50 • C, prior to being discharged. The modeling of a shaft furnace, simulating the reduction of iron ore by syngas, is a powerful tool for defining optimal operating conditions. Use of such a model can lead to the maximization of conversion or the minimization of energy consumption, among other effects capable of reducing carbon dioxide emissions. As such, numerous iron ore shaft furnace models have been proposed in the literature. Initial studies addressed the reduction of a single pellet by H2, CO, or H2-CO mixtures [4][5][6][7][8][9]. Subsequent studies developed models that simulated the reduction zone of the shaft furnace in one dimension [10,11]. With the aim of correctly describing the lateral gas feed, some studies have introduced two-dimensional models [12][13][14]; however, these models did not consider the presence of methane, which is responsible for important reactions in the process. More recently, several authors introduced other reactions [15] and accounted for the cooling zone [16,17]. Some even developed plant models [18]; however, these works were limited to one-dimensional models.
In this work, we developed further the model of Ranzani Da Costa and Wagner, built to simulate the reduction section of DR shafts, operated with pure hydrogen [13,14,19]. We extended this model to consider CO-H2-CH4 reducing gas, and accounted for transition and cooling sections. The present model, named REDUCTOR, is 2-dimensional in the steady-state regime. The model includes a sophisticated, pellet sub-model. We consider eight heterogeneous and two homogeneous chemical reactions. These features represent a more advanced and detailed model, compared to previous studies. Moreover, the results were validated against two sets of plant data.
The present model, REDUCTOR, differs from the other model we recently reported [18] on the following points. REDUCTOR is a computational fluid dynamics (CFD)-type, two-dimensional model, which describes the shaft furnace alone. The other model is of the systemic type, is onedimensional, and aims to simulate the whole DR plant. The shaft furnace description included in the plant model, though based on similar equations, was intentionally made simpler and faster to run, on process simulation software. Thus, REDUCTOR is more detailed and more precise, but, of course, requires longer computation times. The modeling of a shaft furnace, simulating the reduction of iron ore by syngas, is a powerful tool for defining optimal operating conditions. Use of such a model can lead to the maximization of conversion or the minimization of energy consumption, among other effects capable of reducing carbon dioxide emissions. As such, numerous iron ore shaft furnace models have been proposed in the literature. Initial studies addressed the reduction of a single pellet by H 2 , CO, or H 2 -CO mixtures [4][5][6][7][8][9]. Subsequent studies developed models that simulated the reduction zone of the shaft furnace in one dimension [10,11]. With the aim of correctly describing the lateral gas feed, some studies have introduced two-dimensional models [12][13][14]; however, these models did not consider the presence of methane, which is responsible for important reactions in the process. More recently, several authors introduced other reactions [15] and accounted for the cooling zone [16,17]. Some even developed plant models [18]; however, these works were limited to one-dimensional models.
In this work, we developed further the model of Ranzani Da Costa and Wagner, built to simulate the reduction section of DR shafts, operated with pure hydrogen [13,14,19]. We extended this model to consider CO-H 2 -CH 4 reducing gas, and accounted for transition and cooling sections. The present model, named REDUCTOR, is 2-dimensional in the steady-state regime. The model includes a sophisticated, pellet sub-model. We consider eight heterogeneous and two homogeneous chemical reactions. These features represent a more advanced and detailed model, compared to previous studies. Moreover, the results were validated against two sets of plant data.
The present model, REDUCTOR, differs from the other model we recently reported [18] on the following points. REDUCTOR is a computational fluid dynamics (CFD)-type, two-dimensional model, which describes the shaft furnace alone. The other model is of the systemic type, is one-dimensional, and aims to simulate the whole DR plant. The shaft furnace description included in the plant model, though based on similar equations, was intentionally made simpler and faster to run, on process simulation software. Thus, REDUCTOR is more detailed and more precise, but, of course, requires longer computation times.

Principle
The reduction of hematite ore to iron occurs via two intermediate oxides, namely, magnetite and wüstite (considered as Fe 0.95 O [19]), and by two gaseous reactants, namely, H 2 and CO. The following six reduction reactions were therefore considered: Methane reforming and water gas shift reactions also occur in the gas phase, based on the composition of reduction gas and temperature, through the following reactions: We also considered two other side reactions that could occur in the reactor, especially where an iron layer has formed:

•
Methane decomposition reaction CH 4(g) C (s) + 2H 2(g) (9) • Carbon monoxide disproportionation (inverse Boudouard reaction) The model itself is two-dimensional, axisymmetrical, and steady-state. It is based on the numerical solution of local mass, energy, and momentum balances, using the finite volume method. The geometry in the reduction and transition sections is cylindrical, while conical in the cooling section. This corresponds to the geometry of the shaft furnaces and is necessary to describe correctly the lateral gas feed and outlet cooling gas, as shown in Figure 2. The reactor modeled is a shaft furnace of the MIDREX type.   The solid load is fed from the top of the reactor (z = H) to form a moving bed of solid particles composed of spherical iron ore pellets that descend by gravity. The pellet diameter (d p ) is assumed to be unique and unchanging during the reduction reaction, and the initial pellet composition is known. The gas phase is composed of six species: H 2 , CO, H 2 O, CO 2 , N 2 , and CH 4 . The reducing gas is injected from the sidewall, at a height of z = H Feed,gas which then moves upward, against the solid flow, before finally exiting the furnace at the top. The temperature and composition of this reducing gas are known. A secondary feed gas-the cooling gas-which is introduced from the bottom of the furnace (z = −H inf ), is also considered. This cooling gas exits the furnace, through the wall in the upper part of the conical section. The temperatures of the solid and gas are different and vary, according to their position (r, z) within the furnace. The solid temperature is assumed to be uniform in the interior of the pellets. Thus, this model is based on a faithful description of the physical-chemical and thermal phenomena, from the reactor scale to the crystallite scale, as shown in Figure 2. In the pellet sub-model, the pellet is assumed to be initially comprised of dense grains; these grains later fragment into smaller crystallites at the wüstite stage, in agreement with microscopic observations [19]. Thus, from the reactor to the crystallites, we have a 4-scale model.

Gas Phase
The descending solid pellets, through which the ascending gas flows, can be considered a porous medium, consisting of quasi-stationary solid spheres (the gas velocity is much greater than that of the solid). The Ergun equation (see Appendix A for nomenclature), combined with the continuity equation, thus gives 1 r ∂ ∂r where the terms are in units, mol m −3 s −1 ; K is the permeability coefficient, calculated as Materials 2018, 11, 1865 5 of 16 and the source term S mol,tot corresponds to the net gas production by the non-equimolar reactions. Equation (11) is used to calculate the pressure field, and the gas velocity vector is calculated, using Equation (13): The mass balance for a gaseous species, i, considering axial and radial dispersion, in addition to convection, is written: with the source term, S i , given in Table 1. Table 1. Source terms for the gas species mass balances.
The heat balance for the gas phase-considering convection and conduction-as well as the heat exchanged with the solid and heat brought by the gases evolving from the solid, gives:

Solid Phase
Regarding the grain flow, in the upper cylindrical section, it is considered that pellets descend vertically. In contrast, in the lower section of the conical shape, a radial component of the solid velocity must be introduced. A bibliographical study of granular flows led us to use the model of Mullins [20,21], in which the radial velocity is calculated as proportional to the radial gradient of the axial velocity: where B is taken, as proposed in Reference [20]: The mass balance for a gaseous species j gives: with the source term, S j , given in Table 2.
The heat balance for the solid phase takes into account axial and radial convection, conduction, and heat exchange with the gas phase. The heat of the reactions is attributable to the solid phase, considering that all the reactions occur either inside the pellets (heterogeneous reactions) or at their surfaces (homogeneous reactions, catalyzed by the solid); thus:

Transport Coefficients
The various transport coefficients, D r , D z , λ e f f ,r , and λ e f f ,z , as well as other parameters, like specific heats, are calculated as functions of temperature and composition. Details regarding the relationships were given in [22].

Iron Oxide Reduction
Unlike most of the previous approaches published which are based on the shrinking core model (with one or three fronts, separating the oxides in the pellet), we developed a specific, pellet sub-model. The sub-model was built, according to our experimental findings, to simulate the reduction of a single pellet by H 2 -CO. The reaction rate was used as a function of the local reduction conditions (temperature and gas composition), inside the reactor. We used the law of additive reaction times [23], which considers the different resistances (chemical reaction, diffusion, external transfer) involved in series. Therefore, the time required to attain a certain conversion is approximately the sum of the characteristic times: τ i [14,23]. This sub-model was initially developed for simulating reduction by H 2 only, as detailed previously [14]; we extended this model for reduction by CO. The characteristic times and the reaction rates are listed in Appendix B.

Methane Reforming and Water Gas Shift Reactions
Methane reforming and water gas shift reactions are known to be catalyzed by iron or iron oxides [24,25]; thus, their rates are functions of the composition of the reduction gas, temperature, and mass of the catalyst. The methane reforming rate equation considering the forward and reverse reactions is given by Equation (20): The expression of the reaction rate constant, k 7 , is given in Table 3. Because the reforming of CH 4 was hardly observed on the iron oxide catalysts, as reported in the literature [25], it was considered that such reforming only occurs with iron as a catalyst. We assumed that sufficient iron was formed on the outside of the pellet, when the reduction degree exceeded 50%.  [16,26] Similarly, the rate expression for the water gas shift reaction is given by Equation (21): when occurring on Fe or Fe 0.95 O, and by Equation (22) v when occurring on Fe 2 O 3 or Fe 3 O 4 . Here, besides iron, various iron oxides also catalyze the reaction. The corresponding expressions for k 8 and k 8 are given in Table 3, according the literature [24,25].

Carbonization Reactions
In the DR furnace, carbon can be formed, either from methane decomposition (Equation (9)) or from CO disproportionation (Equation (10)). Both reactions are reversible, and the reverse reactions are functions of the carbon activity. The carbon activity was calculated from Chipman's relationship [27]: where C is the ratio of atomic C to atomic Fe. For sake of simplicity, we did not distinguish between C and Fe 3 C in the solid, with both being considered as C.
The rate equation of the methane decomposition reaction is given by Equation (24) v 9 = k 9 P 0.5 The expression of the reaction rate constant, k 9 , included in Equation (24) was determined, as per the literature [16,26], as listed in Table 3. The rate equation of the carbon monoxide disproportionation reaction is given by Equation (25) v 10 = k 10 P 0.5 and the values of the reaction rate constants, k 10 and k 10 , are also provided in Table 3, from the same references.

Boundary Conditions
The balance equations need a set of associated boundary conditions to be solved. First, the temperature and composition of the solids and gases are assumed to be known (the operating conditions) at their respective inlets (at the top for solids, and at the bottom and sides for gases). In addition, because of axisymmetry and tight walls, one has: − Symmetry axis : zero fluxes ∂T s ∂r = For the gas flow, a known pressure condition is also required at the exits. The top pressure was known but not the pressure of the cooling gas outlet, as shown in Figure 3, at Point 4. The latter was estimated to obtain approximately 90% of the inlet cooling gas expelled from this outlet and, approximately, 10% flowing upwards.

Meshing and Numerical Solution
The system of partial derivative equations was discretized and solved, according to the finite volume method [28]. Meshing of the cylindrical reduction and transition sections is orthogonal, with cells made finer next to the top, as shown in Figure 4, on the left. For the conical section, a nonorthogonal grid was used, as shown in Figure 4, on the right. To easily connect the two sections, the number of radial cells was kept the same. The numerical code was written in the language of FORTRAN 1995.  Figure 3 shows the values of the known boundary conditions for the two simulations conducted, corresponding to two different plants. Plant A is a North American, MIDREX plant, currently in operation, the main operating data of which were provided to us. Plant B was the first MIDREX plant operated in the USA, for which published data are available [10]. The production capacity of plant A is 4.5 times greater that of plant B.

Meshing and Numerical Solution
The system of partial derivative equations was discretized and solved, according to the finite volume method [28]. Meshing of the cylindrical reduction and transition sections is orthogonal, with cells made finer next to the top, as shown in Figure 4, on the left. For the conical section, a non-orthogonal grid was used, as shown in Figure 4, on the right. To easily connect the two sections, the number of radial cells was kept the same. The numerical code was written in the language of FORTRAN 1995.

Meshing and Numerical Solution
The system of partial derivative equations was discretized and solved, according to the finite volume method [28]. Meshing of the cylindrical reduction and transition sections is orthogonal, with cells made finer next to the top, as shown in Figure 4, on the left. For the conical section, a nonorthogonal grid was used, as shown in Figure 4, on the right. To easily connect the two sections, the number of radial cells was kept the same. The numerical code was written in the language of FORTRAN 1995.

Results and Discussion
In this section, the results of the Plant A simulation are first presented and discussed, then a comparison between the calculated and measured data for both plants is given. Results for the values of the different variables, throughout the reactor, are given in separate figures; however, all of these variables must be considered simultaneously for interpretation purposes. Figure 5a shows the pressure and velocity fields inside the bed, throughout the reactor. The color scale refers to the pressure, and the lines refer to the streamlines. The large arrows indicate the locations of the various gas and solid inlets and outlets. These locations are the same (and not repeated everywhere) in the following figures. The pressure decreases almost linearly, from bottom to top. The reducing gas, injected at the sidewall (z = 5.32 m), enters radially and then flows essentially vertically, except in the transition zone. The cooling gas first flows upwards, and then, most of it leaves the furnace radially, at the cooling gas outlet, except for a fraction that rises in the reduction section. Figure 5b,c show the temperature distribution of the gas and solid phases in the reactor. First, it was found that the gas and solid temperatures were very close to each other. This similarity resulted from the high gas-to-solid heat transfer, as was described in a previous study [14]. Downwards from the solid inlet, the solid temperature rapidly increased to reach the gas temperature. Second, the temperatures were not axially or radially uniform, throughout the reactor. The hottest zone was near the reducing gas inlet, with gas introduced at 957 • C. Above this inlet, the temperature decreased, because of methane reforming (as shown later, in Figure 7), an endothermic reaction. Third, the cooling gas not only cooled the solid in the bottom section but also influenced the temperature field in the reduction section, with the gas rising from the cooling zone to the central part of the reduction zone. This maintained a lower temperature alongside the center of the shaft.

Pressure Field, Velocity of Gas and Temperature Field
From these results, radial gradients of temperature were revealed to influence, together with the gas composition profiles, the reduction of the solids and the metallization degree achieved.
comparison between the calculated and measured data for both plants is given. Results for the values of the different variables, throughout the reactor, are given in separate figures; however, all of these variables must be considered simultaneously for interpretation purposes. Figure 5a shows the pressure and velocity fields inside the bed, throughout the reactor. The color scale refers to the pressure, and the lines refer to the streamlines. The large arrows indicate the locations of the various gas and solid inlets and outlets. These locations are the same (and not repeated everywhere) in the following figures. The pressure decreases almost linearly, from bottom to top. The reducing gas, injected at the sidewall (z = 5.32 m), enters radially and then flows essentially vertically, except in the transition zone. The cooling gas first flows upwards, and then, most of it leaves the furnace radially, at the cooling gas outlet, except for a fraction that rises in the reduction section.  Figure 5b,c show the temperature distribution of the gas and solid phases in the reactor. First, it was found that the gas and solid temperatures were very close to each other. This similarity resulted from the high gas-to-solid heat transfer, as was described in a previous study [14]. Downwards from the solid inlet, the solid temperature rapidly increased to reach the gas temperature. Second, the temperatures were not axially or radially uniform, throughout the reactor. The hottest zone was near the reducing gas inlet, with gas introduced at 957 °C. Above this inlet, the temperature decreased, because of methane reforming (as shown later, in Figure 7), an endothermic reaction. Third, the cooling gas not only cooled the solid in the bottom section but also influenced the temperature field in the reduction section, with the gas rising from the cooling zone to the central part of the reduction zone. This maintained a lower temperature alongside the center of the shaft.

Pressure Field, Velocity of Gas and Temperature Field
From these results, radial gradients of temperature were revealed to influence, together with the gas composition profiles, the reduction of the solids and the metallization degree achieved.  Figure 6 plots the evolution of solid mass fractions, throughout the reactor. Figure 6a shows that the hematite was fully converted to magnetite very rapidly in the upper part of the reactor. Subsequently, magnetite was reduced to wüstite, as shown in Figure 6b. Afterwards, wüstite slowly began to reduce to iron, as seen in Figure 6c,d. In the external two-thirds of the reduction section, above the reducing gas inlet-a zone where the gas was rich in H 2 and CO and the temperature high-the conversion to iron was completed, in approximately 7 m. In the central part of the reactor, where the temperature was lower and the gas, lower in H 2 and CO, the conversion was not completed and some wüstite remained in the cooling zone. Though the average metallization degree was approximately 94%, metallization was not uniform, with most pellets being completely reduced, whereas others were not. Figure 6e shows the carbon mass fraction, throughout the reactor. We observed that the carbon was in the same location as Fe, in accordance with the catalytic effect of iron on carbon formation. above the reducing gas inlet-a zone where the gas was rich in H2 and CO and the temperature high-the conversion to iron was completed, in approximately 7 m. In the central part of the reactor, where the temperature was lower and the gas, lower in H2 and CO, the conversion was not completed and some wüstite remained in the cooling zone. Though the average metallization degree was approximately 94%, metallization was not uniform, with most pellets being completely reduced, whereas others were not.  Figure 6e shows the carbon mass fraction, throughout the reactor. We observed that the carbon was in the same location as Fe, in accordance with the catalytic effect of iron on carbon formation.

Gas Mole Fractions
As showcased by Figure 7, the situation here is more complex, due to the numerous reactions occurring. The main features of these reactions are as follows. Near the reducing gas inlet, the reforming of methane occurred, which increased the H2 and CO contents. Above the gas inlet, the H2 and CO contents decreased, while H2O and CO2 were formed, as a result of the reduction reactions. In the central zone, with less reduction, lower amounts of H2O and CO2 were formed, and part of the cooling gas, rich in CH4, was present.

Gas Mole Fractions
As showcased by Figure 7, the situation here is more complex, due to the numerous reactions occurring. The main features of these reactions are as follows. Near the reducing gas inlet, the reforming of methane occurred, which increased the H 2 and CO contents. Above the gas inlet, the H 2 and CO contents decreased, while H 2 O and CO 2 were formed, as a result of the reduction reactions. In the central zone, with less reduction, lower amounts of H 2 O and CO 2 were formed, and part of the cooling gas, rich in CH 4 , was present.  Figure 8 is a summary diagram, based on the above results. The shaft furnace was divided into eight zones and distinguished according to the main chemical and thermal processes occurring. On the left part of the diagram are indicated the molar percentages of H2 and CO, involved in each reaction, and the molar percentage of methane, reformed by H2O or CO2, or decomposed to carbon and H2. This diagram is an illustration of how modeling work can help one to understand the detailed behavior of a reactor. Clearly, these results could not be obtained from other means.  Figure 8 is a summary diagram, based on the above results. The shaft furnace was divided into eight zones and distinguished according to the main chemical and thermal processes occurring. On the left part of the diagram are indicated the molar percentages of H 2 and CO, involved in each reaction, and the molar percentage of methane, reformed by H 2 O or CO 2 , or decomposed to carbon and H 2 . This diagram is an illustration of how modeling work can help one to understand the detailed behavior of a reactor. Clearly, these results could not be obtained from other means.  Figure 8 is a summary diagram, based on the above results. The shaft furnace was divided into eight zones and distinguished according to the main chemical and thermal processes occurring. On the left part of the diagram are indicated the molar percentages of H2 and CO, involved in each reaction, and the molar percentage of methane, reformed by H2O or CO2, or decomposed to carbon and H2. This diagram is an illustration of how modeling work can help one to understand the detailed behavior of a reactor. Clearly, these results could not be obtained from other means.

Validation
Unfortunately, neither interior measurements of solid or gas temperatures, nor compositions, were available for comparison with the calculations. However, from some published data regarding Plant B, and from plant data measurements from Plant A, an overall validation of the model was possible. Table 4 provides a comparison of the simulation results, with the available plant data. It can be seen that the model reproduced the outlet temperatures and compositions quite satisfactorily. From this strong agreement, obtained by simulations of two plants of differing capacities, the model can be considered validated.

Conclusions
This article presented the modeling and simulation of an iron ore, direct reduction shaft furnace. We developed a new mathematical model, with the aim of introducing a more-detailed description of the chemical processes, compared to previous studies. The model presented is two-dimensional, describes three sections in the shaft, and accounts for eight heterogeneous and two homogeneous reactions. The model was validated against plant data from two MIDREX plants of notably different capacities. From the analysis of the calculated 2D maps of temperature and composition of the gas and solid phases, it was possible to gain new insights into the interior behavior of the shaft furnace and identify different zones, according to the chemical and thermal phenomena occurring. One significant result is the presence of a central zone of the shaft of lesser temperature and conversion.
Such a model can be helpful in: Investigating the influence of various parameters and operating conditions (including the reducing gas composition), comparing different furnace configurations, and suggesting improvements [29]. These investigations will be the subject of a future paper.