Optimization of the Iron Ore Direct Reduction Process through Multiscale Process Modeling

Iron ore direct reduction is an attractive alternative steelmaking process in the context of greenhouse gas mitigation. To simulate the process and explore possible optimization, we developed a systemic, multiscale process model. The reduction of the iron ore pellets is described using a specific grain model, reflecting the transformations from hematite to iron. The shaft furnace is modeled as a set of interconnected one-dimensional zones into which the principal chemical reactions (3-step reduction, methane reforming, Boudouard and water gas shift) are accounted for with their kinetics. The previous models are finally integrated in a global, plant-scale, model using the Aspen Plus software. The reformer, scrubber, and heat exchanger are included. Results at the shaft furnace scale enlighten the role of the different zones according to the physico-chemical phenomena occurring. At the plant scale, we demonstrate the capabilities of the model to investigate new operating conditions leading to lower CO2 emissions.


Introduction
Steel is one of the key materials of today's industrial world. Moreover, its production is characterized by high energy consumption along with important carbon dioxide emissions. World steel production amounts to 6% of anthropogenic CO 2 emissions [1]. Beside the current reference steelmaking route (integrated plant with sinter and coke making, blast furnace, and basic oxygen furnace), less CO 2 emitting alternative routes are attracting greater interest. Direct reduction (DR) in a shaft furnace followed by electric arc steelmaking thus results in 40 to 60% lower CO 2 emissions compared with the reference route [2]. Direct reduction converts solid iron ore pellets into so-called direct reduced iron (DRI), using a mixture of CO and H 2 produced by reforming natural gas. Due to independence from coke and coal imports, sizeable units to meet demand, reduced investment costs, and reduced construction time [3], direct reduction units are becoming more numerous, particularly in countries where natural gas is cheap and abundant. These advantages have spurred increased research activity, namely in modeling and simulation, in order to correctly evaluate these processes and ultimately optimize them. The present work contributes to this effort with the key objective to provide an accurate simulation of the whole DR process, enabling us to propose new avenues for CO 2 emission mitigation. A distinctive feature of our approach is the multiscale nature of modeling, from, say, 1 µm (grain size) to 1 hm (plant size). We thus developed a single pellet model, a shaft model, and a plant model, using different and complementary modeling strategies. These are presented in the next sections, together with references to previous works.

Single Pellet Model
Iron ore pellets (roughly spherical, typically 7-15 mm diameter) for DR are industrially produced from natural hematite grains (irregular, typically 20 µ m). In a shaft furnace, these pellets are chemically reduced at 600-950 °C by both CO and H2 in three steps, hematite Fe2O3 ⟶ magnetite Fe3O4 ⟶ wustite Fe0.95O ⟶ iron Fe, thus involving six gas-solid reactions at the grain scale. Numerous models for gas-solid reactions are available in the literature, from the simplest unreacted shrinking core model to pore models and grain models (Figure 1) [4]. The grain model matches well the grainy structure of the DR pellets; however, it was not used as such in the current study because it requires a numerical solution of the diffusion-reaction equations, which is incompatible for computation time reasons because of its integration in a multi-particle reactor model (next section). Instead, we used an alternate approach, the additive reaction times law, which gives an approximate analytical solution to calculate the reaction rates and can account for mixed regimes [5]. According to experimental observations, we also introduced an evolution of grain and pore structure with the reactions as follows: the hematite grains become covered with small pores after conversion to magnetite and wustite. At the wustite stage, the pores enlarge and the wustite grains break down into smaller particles termed "crystallites". These subsequently grow (typically from 1 to 10 µ m) and join to form the molten-like structure of sponge iron (another name for DRI). These changes of course influence the kinetics of the transformation, mostly via the diffusion resistances. Details of the single pellet model used and of its results are given in [6,7].

Previous Works
The shaft furnace (Figure 2) is the core of the DR process. Iron ore pellets are charged at the top, descend due to gravity, and encounter an upward counter-flow of gas. The reducing gas (CO and H2, plus CH4, CO2, and H2O, at about 950 °C) is injected peripherally at mid-height and exits at the top. In the lower section of the furnace, of conical shape, cold natural gas is injected to cool the iron pellets produced. The upper section (reducing zone and intermediate zone) is cylindrical (typically height 15 m; diameter 5 m). Two processes, MIDREX and HYL-ENERGIRON, share most of the DR market.

Previous Works
The shaft furnace (Figure 2) is the core of the DR process. Iron ore pellets are charged at the top, descend due to gravity, and encounter an upward counter-flow of gas. The reducing gas (CO and H 2 , plus CH 4 , CO 2 , and H 2 O, at about 950 • C) is injected peripherally at mid-height and exits at the top. In the lower section of the furnace, of conical shape, cold natural gas is injected to cool the iron pellets produced. The upper section (reducing zone and intermediate zone) is cylindrical (typically height 15 m; diameter 5 m). Two processes, MIDREX and HYL-ENERGIRON, share most of the DR market. Their shaft furnaces exhibit some differences (mostly in gas composition and pressure, size, and internal equipment details) but their basic characteristics are similar.
Materials 2018, 11, x FOR PEER REVIEW 3 of 18 Their shaft furnaces exhibit some differences (mostly in gas composition and pressure, size, and internal equipment details) but their basic characteristics are similar. Reflecting the great effort put in the simulation and development of the DR processes, numerous mathematical and numerical models of the shaft furnace were published, differing by the assumptions, the physico-chemical phenomena accounted for, the numerical scheme, the sections of the furnace considered, and the use of the model, etc. The most recent ones [6][7][8][9][10][11][12][13][14][15] are classified in Tables 1 and 2 and below. The basic description is that of an axisymmetric, porous moving bed reactor with counter-current flow between an ascending gas and a descending solid. The main differences are as follows: • The number of reduction steps. This is related to the intermediary iron oxides taken into account. Works with one reduction step consider the direct transformation of hematite to iron, those with two steps consider further the presence of wustite, whereas those with three steps consider also magnetite.

•
The nature of the inlet gas mixture. Most authors considered an inlet gas consisting of CO and H2 only. Two models [13,14] considered the actual mixture, with CO2, H2O, N2, and CH4 as additional gas components.

•
The number of dimensions included in the geometrical description. The standard is onedimensional models. However, three works [6,14,15] considered two dimensions: along the height and the radius of the furnace.

•
The description of the pellet transformation: shrinking core or grain models.

•
The presence of supplementary reactions. The reduction reactions (in one, two or three steps) are always included. Additionally, methane decomposition and steam methane reforming were also taken into account in [10,[12][13][14], together with the water gas shift and Boudouard reactions in [14].

•
The type of heat transfer. All works included heat convection between solid and gas, as well as the heat of reaction. Heat transfer by conduction and radiation, through a contribution to effective conductivity, was also sometimes considered [7,8,11,14,15].

•
The presence of the cooling zone was only taken into account in three papers [10,13,14]. • Some models were validated against plant data, others not.
Results that are common to these works are as follows:  Reflecting the great effort put in the simulation and development of the DR processes, numerous mathematical and numerical models of the shaft furnace were published, differing by the assumptions, the physico-chemical phenomena accounted for, the numerical scheme, the sections of the furnace considered, and the use of the model, etc. The most recent ones [6][7][8][9][10][11][12][13][14][15] are classified in Tables 1  and 2 and below. The basic description is that of an axisymmetric, porous moving bed reactor with counter-current flow between an ascending gas and a descending solid. The main differences are as follows:

•
The number of reduction steps. This is related to the intermediary iron oxides taken into account. Works with one reduction step consider the direct transformation of hematite to iron, those with two steps consider further the presence of wustite, whereas those with three steps consider also magnetite.

•
The nature of the inlet gas mixture. Most authors considered an inlet gas consisting of CO and H 2 only. Two models [13,14] considered the actual mixture, with CO 2 , H 2 O, N 2 , and CH 4 as additional gas components.

•
The number of dimensions included in the geometrical description. The standard is one-dimensional models. However, three works [6,14,15] considered two dimensions: along the height and the radius of the furnace.

•
The description of the pellet transformation: shrinking core or grain models.

•
The presence of supplementary reactions. The reduction reactions (in one, two or three steps) are always included. Additionally, methane decomposition and steam methane reforming were also taken into account in [10,[12][13][14], together with the water gas shift and Boudouard reactions in [14].

•
The type of heat transfer. All works included heat convection between solid and gas, as well as the heat of reaction. Heat transfer by conduction and radiation, through a contribution to effective conductivity, was also sometimes considered [7,8,11,14,15].

•
The presence of the cooling zone was only taken into account in three papers [10,13,14]. • Some models were validated against plant data, others not.
Results that are common to these works are as follows: • The molar content of H 2 and CO decreases from the reduction zone bottom to its top whereas that of CO 2 and H 2 O increases. Inversely, the content of iron oxides decreases from top to bottom with mainly iron exiting from the shaft bottom.

•
The solid and gas temperatures are equal along the shaft except in a thin layer at the top near the solid inlet where a great temperature difference exists, the pellets are charged cold. Moreover, both temperatures increase from the shaft top to its bottom. • Reaction enthalpy is of great importance, namely, the rather endothermic nature of H 2 reduction reactions vs. the exothermic nature of CO reduction reactions, as well as for the models that include it, the endothermic nature of steam methane reforming.
Key points that can be deduced when comparing differing results indicate that: • the inclusion of three reduction steps and the grain model better represent the pellet transformation; • all the components of the gas mixture (excepting N 2 but including CH 4 ) play a role in the transformations; • two dimensions depict more accurately the reactor internal behavior and the output results; 2D results revealed the presence of two zones: one peripheral where the bustle gas is the reducing gas, and one central where the gas stems from the cooling and transition zones; • taking supplementary reactions into consideration along with the cooling and transition zones better represents gas phase transformations and can account for carbon deposition.
The present shaft model has been built considering these results and on the basis of the most detailed (2D, 3 zones, 10 reactions-named REDUCTOR) shaft model [6,14]. However, the goal here is to simulate the whole DR plant, thus it was not practicable to incorporate such a comprehensive, 2D, CFD-type model in the plant simulation software. We selected Aspen Plus, a commercial simulation tool widely used in the chemical industry, as the software. Processes like the steelmaking blast furnace and the reverse chemical looping process are examples of processes somewhat close to the DR process that have been modeled with Aspen Plus. Thus, we had to build an Aspen Plus model of the DR shaft derived from REDUCTOR. Nouri et al. [9] reaction, convection no no yes Shams and Moazeni [10] yes yes yes Palacios et al. [11] reaction, convection, diffusion no no no Bayu and Alamsari [12] Alamsari et al. [13] reaction, convection yes no yes Ranzani da Costa et al. [7] reaction, convection, diffusion no no no Hamadeh et al. [6,14] yes yes yes Ghadi et al. [15] no no yes

Aspen Plus Shaft Model
The results given in [14] show that, in addition to the initial 3-stage division, the reduction zone can be sub-divided into two zones. In the first one, Zone 1, the widest, peripheral zone, the oxides are almost completely reduced at the bottom of the zone. The inlet gas in this zone is the major fraction of the bustle reducing gas plus a part of the gas coming upward from the transition zone. The inlet solid is the major fraction of the oxide pellets charged. The second zone, Zone 2, is a narrow central zone in which only partial metallization occurs. Its inlet gas comes from the transition zone, whereas its inlet solid is a fraction of the oxide pellets. Figure 3 depicts this partitioning and also shows its translation into an Aspen Plus based configuration.

Aspen Plus Shaft Model
The results given in [14] show that, in addition to the initial 3-stage division, the reduction zone can be sub-divided into two zones. In the first one, Zone 1, the widest, peripheral zone, the oxides are almost completely reduced at the bottom of the zone. The inlet gas in this zone is the major fraction of the bustle reducing gas plus a part of the gas coming upward from the transition zone. The inlet solid is the major fraction of the oxide pellets charged. The second zone, Zone 2, is a narrow central zone in which only partial metallization occurs. Its inlet gas comes from the transition zone, whereas its inlet solid is a fraction of the oxide pellets. Figure 3 depicts this partitioning and also shows its translation into an Aspen Plus based configuration. Nevertheless, due to the countercurrent of the solid and gas flows, this description implies using several loops for the numerical solution; e.g., the gas coming from the transition zone influences the solid in the central zone, which in turn influences the transition zone gas. Likewise, the solid from the transition zone influences the gas from the cooling zone, which in turn influences the solid. These loops make it difficult and time-consuming to simulate the process. An effort was thus made to simplify the model. The flow rate of the solid in the central zone being considerably smaller than that in the peripheral zone, the former was not considered as an input to the transition zone, thus cancelling the first loop. Moreover, the transition zone had every little impact on the temperature of the descending solid, thus the link between the descending transition solid and the cooling zone was removed. In place, a fictitious solid stream equal to that leaving the peripheral zone was used as an input. The resulting temperature was thus imposed on the solid leaving the transition zone. This solid was considered as the ultimate output of the peripheral zone. Lastly, an additional cooling zone was added for the central zone. Attention was made to avoid any loops between the two sections. Considering this, Figure 4 highlights the restructuring of the Aspen Plus configuration to avoid the aforementioned loops. Herein, the cross signs emphasize deleted streams, which were responsible for loops, whereas the dotted lines indicate new streams. As a result of this restructuring, a new cooling zone (COLD2) was added and the resulting temperature from the cold section was returned to the exit of the transition zone. Of course, these changes, made for the sake of simplification, do not alter the overall mass and heat flowrates and balances. Nevertheless, due to the countercurrent of the solid and gas flows, this description implies using several loops for the numerical solution; e.g., the gas coming from the transition zone influences the solid in the central zone, which in turn influences the transition zone gas. Likewise, the solid from the transition zone influences the gas from the cooling zone, which in turn influences the solid. These loops make it difficult and time-consuming to simulate the process. An effort was thus made to simplify the model. The flow rate of the solid in the central zone being considerably smaller than that in the peripheral zone, the former was not considered as an input to the transition zone, thus cancelling the first loop. Moreover, the transition zone had every little impact on the temperature of the descending solid, thus the link between the descending transition solid and the cooling zone was removed. In place, a fictitious solid stream equal to that leaving the peripheral zone was used as an input. The resulting temperature was thus imposed on the solid leaving the transition zone. This solid was considered as the ultimate output of the peripheral zone. Lastly, an additional cooling zone was added for the central zone. Attention was made to avoid any loops between the two sections. Considering this, Figure 4 highlights the restructuring of the Aspen Plus configuration to avoid the aforementioned loops. Herein, the cross signs emphasize deleted streams, which were responsible for loops, whereas the dotted lines indicate new streams. As a result of this restructuring, a new cooling zone (COLD2) was added and the resulting temperature from the cold section was returned to the exit of the transition zone. Of course, these changes, made for the sake of simplification, do not alter the overall mass and heat flowrates and balances. This representation however poses the problem of finding the adequate split for the streams between the various sections. We identified three critical splits. The first is at the reduction gas entry level, the bustle gas being split between the transition zone and Zone 1. The second is at the top solid inlet, the entering solid being split between Zone 1 and Zone 2. The third split is at the cooling gas stream that goes up to the transition zone. The first split was set to a predefined value. The second split was calculated considering that the ratios of the inlet solid flow to the input gas flow for each zone are equal. The third split was set to the constant value of 0.13 as given in [6].
This representation, set and done, is not yet sufficient to directly calculate the final conversion and temperatures. Especially for the transition and reduction zones, Aspen Plus does not offer any built-in model that allows their calculation in a way that could correctly depict the reaction kinetics and the variations of the different variables along the height of the shaft. An external calculator, written in Fortran and based on REDUCTOR's equations, was therefore used to compute these values, which were later rendered to the corresponding Aspen Plus blocks. Conversely, for the cooling section, a simple Aspen Plus heat exchanger could be used. The principle of the calculation, namely in the reduction and transition zones, and the list of reactions are given in Appendix A.

Results from the Aspen Plus Shaft Model
Two case studies were considered in this work, based on real industrial data. The first case corresponds to the Contrecoeur plant, located near Montreal, Canada, and currently operated by ArcelorMittal. It was constructed in the late 1970's and is a MIDREX series 750 module. It is characterized by a rather cold input solid and a high content of C in the input gas. The second case relates to the Gilmore plant, built near Portland, Oregon, USA. Now decommissioned, it was the first operating MIDREX plant. It is one of the most referenced plants in literature, a MiniMod module with a CDRI production of 26.4 tons/h. The input operating conditions for both plants are given in Appendix B. These are of quite different capacities, with different reducing gas compositions and different pellet diameters, and thus represent a good test for validating a model. Table 3 compares the results obtained in the model for the outlet gas and solid streams with those provided in literature for the studied cases. The difference for each parameter was calculated as well as the total error. As can be seen, the absolute error for the Contrecoeur case is equal to 6%, whereas that for Gilmore is 4%. The biggest differences pertain to methane and nitrogen composition as well temperature. Other differences include hydrogen and water compositions in the Gilmore case, and the carbon dioxide composition in the Contrecoeur case. These differences can be related to the formulas chosen from the literature for the gas phase reaction rates. For example, methane decomposition and Boudouard reactions seem to be somewhat underestimated in Contrecoeur and overestimated in Gilmore. These values could only be further consolidated through the realization of up-to-date experiments to correctly characterize these reactions. Nonetheless, the results seem to be globally satisfactory, especially if they are compared with other model results [10,14]. The related This representation however poses the problem of finding the adequate split for the streams between the various sections. We identified three critical splits. The first is at the reduction gas entry level, the bustle gas being split between the transition zone and Zone 1. The second is at the top solid inlet, the entering solid being split between Zone 1 and Zone 2. The third split is at the cooling gas stream that goes up to the transition zone. The first split was set to a predefined value. The second split was calculated considering that the ratios of the inlet solid flow to the input gas flow for each zone are equal. The third split was set to the constant value of 0.13 as given in [6].
This representation, set and done, is not yet sufficient to directly calculate the final conversion and temperatures. Especially for the transition and reduction zones, Aspen Plus does not offer any built-in model that allows their calculation in a way that could correctly depict the reaction kinetics and the variations of the different variables along the height of the shaft. An external calculator, written in Fortran and based on REDUCTOR's equations, was therefore used to compute these values, which were later rendered to the corresponding Aspen Plus blocks. Conversely, for the cooling section, a simple Aspen Plus heat exchanger could be used. The principle of the calculation, namely in the reduction and transition zones, and the list of reactions are given in Appendix A.

Results from the Aspen Plus Shaft Model
Two case studies were considered in this work, based on real industrial data. The first case corresponds to the Contrecoeur plant, located near Montreal, Canada, and currently operated by ArcelorMittal. It was constructed in the late 1970's and is a MIDREX series 750 module. It is characterized by a rather cold input solid and a high content of C in the input gas. The second case relates to the Gilmore plant, built near Portland, Oregon, USA. Now decommissioned, it was the first operating MIDREX plant. It is one of the most referenced plants in literature, a MiniMod module with a CDRI production of 26.4 tons/h. The input operating conditions for both plants are given in Appendix B. These are of quite different capacities, with different reducing gas compositions and different pellet diameters, and thus represent a good test for validating a model. Table 3 compares the results obtained in the model for the outlet gas and solid streams with those provided in literature for the studied cases. The difference for each parameter was calculated as well as the total error. As can be seen, the absolute error for the Contrecoeur case is equal to 6%, whereas that for Gilmore is 4%. The biggest differences pertain to methane and nitrogen composition as well temperature. Other differences include hydrogen and water compositions in the Gilmore case, and the carbon dioxide composition in the Contrecoeur case. These differences can be related to the formulas chosen from the literature for the gas phase reaction rates. For example, methane decomposition and Boudouard reactions seem to be somewhat underestimated in Contrecoeur and overestimated in Gilmore. These values could only be further consolidated through the realization of up-to-date experiments to correctly characterize these reactions. Nonetheless, the results seem to be globally satisfactory, especially if they are compared with other model results [10,14]. The related differences are comparable although the present model is of a different type and simpler than the two CFD-type models.  Figure 5 shows the evolution of the solid component flow rates with shaft height in the first reduction zone for the Contrecoeur (a) and Gilmore (b) cases. Height zero corresponds to that of the bustle gas inlet. As can be seen, hematite disappears very rapidly near the solid inlet in both cases. Magnetite on the other hand has a different behavior; it disappears after 2.6 m in the Contrecoeur case, against 4.6 m in the Gilmore case, with more of this oxide formed in the latter case. Wustite on the other hand disappears rapidly in the Gilmore case, its presence window being only 2 m, against about 7 m in the Contrecoeur case. Finally, as expected, both cases give a 100% conversion of iron oxide to pure iron in this Zone 1. Figure 5 also shows the evolution of the carbon deposited in the pellets, which continuously increases from solid entrance before dropping near gas entrance. Exit solid has some carbon content in the Contrecoeur case, whereas the carbon content drops to zero in the Gilmore case. This is related to the inversion of the Boudouard equilibrium (C + CO 2 2CO), which is favorable to carbon deposition at lower temperatures and lower CO 2 content. differences are comparable although the present model is of a different type and simpler than the two CFD-type models.  Figure 5 shows the evolution of the solid component flow rates with shaft height in the first reduction zone for the Contrecoeur (a) and Gilmore (b) cases. Height zero corresponds to that of the bustle gas inlet. As can be seen, hematite disappears very rapidly near the solid inlet in both cases. Magnetite on the other hand has a different behavior; it disappears after 2.6 m in the Contrecoeur case, against 4.6 m in the Gilmore case, with more of this oxide formed in the latter case. Wustite on the other hand disappears rapidly in the Gilmore case, its presence window being only 2 m, against about 7 m in the Contrecoeur case. Finally, as expected, both cases give a 100% conversion of iron oxide to pure iron in this Zone 1. Figure 5 also shows the evolution of the carbon deposited in the pellets, which continuously increases from solid entrance before dropping near gas entrance. Exit solid has some carbon content in the Contrecoeur case, whereas the carbon content drops to zero in the Gilmore case. This is related to the inversion of the Boudouard equilibrium (C + CO2 ⇌ 2CO), which is favorable to carbon deposition at lower temperatures and lower CO2 content. These profiles can be related to the gas phase reactions, which are presented for both cases in Figure 6. It can be seen that methane, water, and carbon dioxide flow rates decrease above the gas entrance, whereas those of hydrogen and carbon monoxide increase. This is the opposite in the upper half of the shaft, except for methane, which sees its flow rate reach equilibrium. This inversion can be explained by the preponderance of the iron oxide reduction reactions, as well as the inverse These profiles can be related to the gas phase reactions, which are presented for both cases in Figure 6. It can be seen that methane, water, and carbon dioxide flow rates decrease above the gas entrance, whereas those of hydrogen and carbon monoxide increase. This is the opposite in the upper half of the shaft, except for methane, which sees its flow rate reach equilibrium. This inversion can be explained by the preponderance of the iron oxide reduction reactions, as well as the inverse Boudouard reaction (2CO C + CO 2 ) as evidenced by the carbon production in Figure 5. The decrease in methane can be related mainly to the steam reforming (CH 4 + H 2 O CO + 3H 2 ) with carbon deposition from methane (CH 4 C + 2H 2 ) having a rather negligible impact. This little impact is emphasized by the decline in carbon flow rate near the gas inlet, which is due to the direct Boudouard reaction. in methane can be related mainly to the steam reforming (CH4 + H2O ⇌ CO + 3H2) with carbon deposition from methane (CH4 ⇌ C + 2H2) having a rather negligible impact. This little impact is emphasized by the decline in carbon flow rate near the gas inlet, which is due to the direct Boudouard reaction. Among the other calculated results (temperature, reaction rates, and values for the other zones), we selected for presentation and comparison the evolution of the iron compounds flow rates with shaft height in Zone 2 ( Figure 7). The situation differs from that of Zone 1. Whereas hematite is readily reduced into magnetite as previously reported, magnetite remains present over most of the height, being slowly reduced in the Contrecoeur case and almost not changing in the Gilmore case until the gas inlet. Conversion to wustite and iron only occurs at the zone bottom, with 60% wustite remaining in the solid at the exit. This situation results from a lower temperature and less CO and H2 for the reduction in Zone 2 compared with Zone 1. The other reactions, involving methane, carbon monoxide and carbon dioxide, hardly take place in Zone 2. These results emphasize the impact this second zone has on the low exit gas temperature and on the average metallization degree. The importance of the transition zone should also be noted. The main reaction in this zone is methane decomposition (CH 4 → C + 2H 2 ). Carbon is herein deposited on produced iron, leading to the final carbon content of the DRI. This carbon comes in addition to the carbon possibly deposited via the inverse Boudouard reaction in Zone 1 (Figure 5). At the rather low temperature of the transition zone, no iron oxide reduction by solid carbon occurs. The hydrogen produced by the methane decomposition reaction is sent to the second zone, contributing to the final metallization degree. Moreover, the gas-solid temperature difference is reversed in this zone, the descending solid being hotter than the ascending gas. Among the other calculated results (temperature, reaction rates, and values for the other zones), we selected for presentation and comparison the evolution of the iron compounds flow rates with shaft height in Zone 2 ( Figure 7). The situation differs from that of Zone 1. Whereas hematite is readily reduced into magnetite as previously reported, magnetite remains present over most of the height, being slowly reduced in the Contrecoeur case and almost not changing in the Gilmore case until the gas inlet. Conversion to wustite and iron only occurs at the zone bottom, with 60% wustite remaining in the solid at the exit. This situation results from a lower temperature and less CO and H 2 for the reduction in Zone 2 compared with Zone 1. The other reactions, involving methane, carbon monoxide and carbon dioxide, hardly take place in Zone 2. These results emphasize the impact this second zone has on the low exit gas temperature and on the average metallization degree. in methane can be related mainly to the steam reforming (CH4 + H2O ⇌ CO + 3H2) with carbon deposition from methane (CH4 ⇌ C + 2H2) having a rather negligible impact. This little impact is emphasized by the decline in carbon flow rate near the gas inlet, which is due to the direct Boudouard reaction. Among the other calculated results (temperature, reaction rates, and values for the other zones), we selected for presentation and comparison the evolution of the iron compounds flow rates with shaft height in Zone 2 ( Figure 7). The situation differs from that of Zone 1. Whereas hematite is readily reduced into magnetite as previously reported, magnetite remains present over most of the height, being slowly reduced in the Contrecoeur case and almost not changing in the Gilmore case until the gas inlet. Conversion to wustite and iron only occurs at the zone bottom, with 60% wustite remaining in the solid at the exit. This situation results from a lower temperature and less CO and H2 for the reduction in Zone 2 compared with Zone 1. The other reactions, involving methane, carbon monoxide and carbon dioxide, hardly take place in Zone 2. These results emphasize the impact this second zone has on the low exit gas temperature and on the average metallization degree. The importance of the transition zone should also be noted. The main reaction in this zone is methane decomposition (CH 4 → C + 2H 2 ). Carbon is herein deposited on produced iron, leading to the final carbon content of the DRI. This carbon comes in addition to the carbon possibly deposited via the inverse Boudouard reaction in Zone 1 (Figure 5). At the rather low temperature of the transition zone, no iron oxide reduction by solid carbon occurs. The hydrogen produced by the methane decomposition reaction is sent to the second zone, contributing to the final metallization degree. Moreover, the gas-solid temperature difference is reversed in this zone, the descending solid being hotter than the ascending gas. The importance of the transition zone should also be noted. The main reaction in this zone is methane decomposition (CH 4 → C + 2H 2 ). Carbon is herein deposited on produced iron, leading to the final carbon content of the DRI. This carbon comes in addition to the carbon possibly deposited via the inverse Boudouard reaction in Zone 1 (Figure 5). At the rather low temperature of the transition zone, no iron oxide reduction by solid carbon occurs. The hydrogen produced by the methane decomposition reaction is sent to the second zone, contributing to the final metallization degree. Moreover, the gas-solid temperature difference is reversed in this zone, the descending solid being hotter than the ascending gas.

Plant Model
The next scale is that of the whole DR plant, a schematic diagram of which is given in Figure 8 (case of a MIDREX-type plant). The main units around the shaft furnace are: the reformer, the scrubber, and the heat exchanger. Downstream from the shaft furnace the top gas is scrubbed and divided into two streams. The first stream loops in the system. It is mixed with natural gas, reheated, and sent to the catalytic tubes of the reformer to produce the reducing gas. The second one (plus possibly some extra natural gas) is burnt with air to provide the reformer with heat. The hot flue gas is used in an exchanger to heat both gas streams entering the reformer. The MIDREX reformer is a kind of 'dry' reformer where the primary reaction is CH 4 + CO 2 −→ 2H 2 + 2CO. In other DR processes, methane is rather reformed by H 2 O (HYL), the reformer is smaller or missing (HYL-ZR) or the reducing gas is produced from other fuels (ENERGIRON). This variety was a further incentive to develop the present DR plant model using Aspen Plus, since we expect to use it for process optimization, including modifying parts of the configuration. Contrary to the shaft furnace, the plant scale was scarcely investigated through modeling. Only two authors studied the interaction between the reformer and the shaft furnace, using mathematical models for each unit [6,16].

Plant Model
The next scale is that of the whole DR plant, a schematic diagram of which is given in Figure 8 (case of a MIDREX-type plant). The main units around the shaft furnace are: the reformer, the scrubber, and the heat exchanger. Downstream from the shaft furnace the top gas is scrubbed and divided into two streams. The first stream loops in the system. It is mixed with natural gas, reheated, and sent to the catalytic tubes of the reformer to produce the reducing gas. The second one (plus possibly some extra natural gas) is burnt with air to provide the reformer with heat. The hot flue gas is used in an exchanger to heat both gas streams entering the reformer. The MIDREX reformer is a kind of 'dry' reformer where the primary reaction is CH4 + CO2 ⟶ 2H2 + 2CO. In other DR processes, methane is rather reformed by H2O (HYL), the reformer is smaller or missing (HYL-ZR) or the reducing gas is produced from other fuels (ENERGIRON). This variety was a further incentive to develop the present DR plant model using Aspen Plus, since we expect to use it for process optimization, including modifying parts of the configuration. Contrary to the shaft furnace, the plant scale was scarcely investigated through modeling. Only two authors studied the interaction between the reformer and the shaft furnace, using mathematical models for each unit [6,16].  Figure 9 illustrates the corresponding Aspen Plus model. Starting from the shaft furnace top gas exit the scrubber is the first step, in which part of the water vapor is separated from the rest of the stream. The adopted Aspen sub-model is a splitter, in which the separated H2O fraction is defined ( 2 ). The remaining off-gas going is then split in two streams: one ( ) going to the reformer and the other to the burner. On the other hand, reformer natural gas is also split in two fractions, one sent prior to the reformer ( 4 ) and the other after it. The reformer was modeled using the built-in Aspen Plus 'Gibbs reactor'. This reactor does not consider kinetics, but calculates the product composition minimizing its Gibbs energy; i.e., at equilibrium. This is a simplification often used in the literature [18,19]. The key design variable is Tref, the reforming temperature. The gas stream exiting the reformer is later mixed with natural gas and air prior to being sent back to the shaft furnace. A prior step was added in the form of a Gibbs reactor where O2 injection and burning is realized. Herein, the reactor temperature ( , 2 ) is the key design variable. The combustion system on the other hand has as inlets the remaining off-gas, natural gas, and combustion air. These streams are heated to the combustion temperature, burnt in the burner producing hot flue gases later cooled prior to discharge. The burner was also modeled as a stoichiometric reactor, employing gas  Figure 9 illustrates the corresponding Aspen Plus model. Starting from the shaft furnace top gas exit the scrubber is the first step, in which part of the water vapor is separated from the rest of the stream. The adopted Aspen sub-model is a splitter, in which the separated H 2 O fraction is defined (split H 2 O ). The remaining off-gas going is then split in two streams: one (split gas ) going to the reformer and the other to the burner. On the other hand, reformer natural gas is also split in two fractions, one sent prior to the reformer (split CH 4 ) and the other after it. The reformer was modeled using the built-in Aspen Plus 'Gibbs reactor'. This reactor does not consider kinetics, but calculates the product composition minimizing its Gibbs energy; i.e., at equilibrium. This is a simplification often used in the literature [18,19]. The key design variable is T ref , the reforming temperature. The gas stream exiting the reformer is later mixed with natural gas and air prior to being sent back to the shaft furnace. A prior step was added in the form of a Gibbs reactor where O 2 injection and burning is realized. Herein, the reactor temperature (T r,O 2 ) is the key design variable. The combustion system on the other hand has as inlets the remaining off-gas, natural gas, and combustion air. These streams are heated to the combustion temperature, burnt in the burner producing hot flue gases later cooled prior to discharge. The burner was also modeled as a stoichiometric reactor, employing gas combustion reactions, with its temperature (T bur ) also as a key design variable. The burner's goal is to provide energy for the reforming system. This step is then followed by flue gas cooling in order to recover energy, namely for pre-heating purposes. The key flue gas cooling variable is (T f lue ).

Aspen Plus Plant Model
Materials 2018, 11, x FOR PEER REVIEW 11 of 18 to provide energy for the reforming system. This step is then followed by flue gas cooling in order to recover energy, namely for pre-heating purposes. The key flue gas cooling variable is ( ).

Figure 9.
Whole plant Aspen Plus model, with gas recirculation and key design variables.
The principle used for modeling the plant in Aspen Plus was to make use of Aspen Plus Design Specs. A Design Spec is a condition that the Aspen Plus solver will make satisfied by acting on a 'controlled variable'. In this case, on the one hand, reformer variables were controlled so that the gas exiting the reformer is identical to the bustle gas, leading to a closed looped process. On the other hand, burner variables were controlled in order to guarantee an adequate heat balance in the process. The design specs refer to the equality of flow rates for hydrogen, water vapor, carbon monoxide, carbon dioxide, methane, and nitrogen between the bustle-gas and the exit gas respectively. The controlled variables are: , and , . Table 4 lists the various design specifications as well as the corresponding controlled variables, and the resulting values. The obtained values were compared to available data in the case of Contrecoeur; no such data were available for Gilmore. As it can be seen, the values do not differ in the first case except for 4 . Indeed, this is not really a discrepancy given that the same amount of methane is reformed. According to plant data, 92% of the methane is sent to the reformer where 46% of it is reformed, the yield being limited by kinetics. Conversely, the reformer modeled as a Gibbs reactor gives an almost complete reforming and thus only 43% of the methane is needed in the reformer, the rest being by-passed. Another approach would have been to adopt a 1D model of the reformer with kinetics, like [6,16]. Downstream of the reformer, the gas temperature is adjusted using a small air injection. The principle used for modeling the plant in Aspen Plus was to make use of Aspen Plus Design Specs. A Design Spec is a condition that the Aspen Plus solver will make satisfied by acting on a 'controlled variable'. In this case, on the one hand, reformer variables were controlled so that the gas exiting the reformer is identical to the bustle gas, leading to a closed looped process. On the other hand, burner variables were controlled in order to guarantee an adequate heat balance in the process. The design specs refer to the equality of flow rates for hydrogen, water vapor, carbon monoxide, carbon dioxide, methane, and nitrogen between the bustle-gas and the exit gas respectively. The controlled variables are: split gas , split H 2 O , T re f , T r,O 2 , n CH 4 ,re f , n air , split CH 4 , n CH 4 ,bur , T f lue and n air,bur . Table 4 lists the various design specifications as well as the corresponding controlled variables, and the resulting values. The obtained values were compared to available data in the case of Contrecoeur; no such data were available for Gilmore. As it can be seen, the values do not differ in the first case except for split CH 4 . Indeed, this is not really a discrepancy given that the same amount of methane is reformed. According to plant data, 92% of the methane is sent to the reformer where 46% of it is reformed, the yield being limited by kinetics. Conversely, the reformer modeled as a Gibbs reactor gives an almost complete reforming and thus only 43% of the methane is needed in the reformer, the rest being by-passed. Another approach would have been to adopt a 1D model of the reformer with kinetics, like [6,16]. Downstream of the reformer, the gas temperature is adjusted using a small air injection. Results in the Gilmore case resemble the Contrecoeur case for T re f and n CH 4 .bur . The first is related to equilibrium whereas the second showcases the lack of need for excess fuel. Both processes however differ for the remaining variables. In the Gilmore case, split gas has a higher value leading to greater off-gas recycling. This may indicate the need for greater CO 2 reforming as well as the need to conserve unreacted CO-H 2 reductants present in the off-gas. The smaller split H 2 O value can be understood in light of a greater H 2 O reforming in the connected reformer. The smaller amount of natural gas sent to the reformer is related to the smaller gas flow rate entering the shaft furnace. The equivalent input air flow rate can be related to the greater amount of nitrogen present in the gas entering the furnace. Higher split CH 4 values can be understood in light of the equilibrium that is reached in the reformer but also a smaller methane fraction in the gas stream entering the shaft furnace. The greater amount of input process air n air,re f required in the burner is related to the greater N 2 flow rate in the bustle gas.

Process Visualization
The adopted process model and simulation tool further enable us to create a mimic board of the plant. Figure 10 highlights the flow diagrams of this process in the Gilmore case. As can be seen, all process units are shown, the shaft reactor, scrubber, reformer, oxygen injection, burner, and air preheater. Moreover, all streams to and from these units are highlighted along with their composition. The figure puts a great emphasis on the gas loop, by providing temperature and composition changes along the operations; e.g., the reduction in H 2 , CO and CH 4 content in the shaft reactor and their recovery in the reforming system. Moreover, it can be seen that the burner provides the reformer system with heat by burning part of the scrubbed off-gas, with little need for input methane. Certain key parameters are shown in yellow, like the mass flow rate of DRI, its metallization and carbon content, together with the ratio of equivalent carbon present in the flue gas to the quantity of DRI produced (C/DRI, here equal to 0.119) as well as the same ratio but for equivalent CO 2 (CO 2 /DRI). In fine, this visualization makes it easier to compare simulations to plant data or to compare simulation runs between each other.

Process Optimization
Having modeled the shaft reactor as well as the gas loop system, various process optimization schemes can be considered. Our aim in this paper was to address the reduction in CO2 emissions. Keeping the same plant configuration (MIDREX-type), the chosen means was the modification of the inlet reducing gas composition. According to the literature, the two important ratios are H 2 CO ⁄ and (H 2 + CO) (H 2 O + CO 2 ) ⁄ ; i.e., the hydrogen-to-carbon monoxide ratio and the reducing power of the reducing gas. The initial values for these parameters were respectively equal to 1.5 and 12 for Contrecoeur and 1.75 and 8.73 for Gilmore. Optimal ratios were recorded close to 1 and 12 for both parameters respectively [9,20]. Considering this, optimization runs are provided hereafter for Gilmore. The initial value for the carbon-to-DRI ratio was 0.119, a figure also in line with the literature, where the range of 0.105-0.120 was reported [2]. ⁄ ratio over 12 did not further reduce the carbon-to-DRI ratio.

Process Optimization
Having modeled the shaft reactor as well as the gas loop system, various process optimization schemes can be considered. Our aim in this paper was to address the reduction in CO 2 emissions. Keeping the same plant configuration (MIDREX-type), the chosen means was the modification of the inlet reducing gas composition. According to the literature, the two important ratios are H 2 /CO and (H 2 + CO)/(H 2 O + CO 2 ); i.e., the hydrogen-to-carbon monoxide ratio and the reducing power of the reducing gas. The initial values for these parameters were respectively equal to 1.5 and 12 for Contrecoeur and 1.75 and 8.73 for Gilmore. Optimal ratios were recorded close to 1 and 12 for both parameters respectively [9,20]. Considering this, optimization runs are provided hereafter for Gilmore. The initial value for the carbon-to-DRI ratio was 0.119, a figure also in line with the literature, where the range of 0.105-0.120 was reported [2]. Table 5 highlights the optimization runs carried out by decreasing H 2 /CO from 1.75 to 1.05 and increasing (H 2 + CO)/(H 2 O + CO 2 ) from 8 to 16. As can be seen, the carbon-to-DRI ratio decreases by as much as 12% between the original case and the optimal case. This case corresponded to a H 2 /CO ratio of 1.23 and (H 2 + CO)/(H 2 O + CO 2 ) of 12. Attempts to further decrease H 2 /CO led to model divergence, the reformer not being able to produce the desired bustle gas composition. Increasing the (H 2 + CO)/(H 2 O + CO 2 ) ratio over 12 did not further reduce the carbon-to-DRI ratio. These results illustrate the potential of the tool. The key differences between the optimized and the original operations pertain to greater split H 2 O (0.72 vs. 0.44) and split gas ratios (0.75 vs. 0.64). More H 2 O is thus withdrawn in the scrubber, whereas greater gas recycling occurs. This can be translated in a greater importance for CO 2 reforming and a greater conservation of the reducing gases. This conservation further leads to a smaller requirement of reformer natural gas (64 mol/s vs. 71 mol/s). Also, a smaller energy demand occurs in the reformer (22 MW vs. 24 MW). It should be further noted that these results were obtained without little (but positive) changes in metallization degree (93.3% for the optimal case vs. 93% for the original case) and carbon content (2% vs. 1.8%).
Calculations were also performed for the Contrecoeur case. The best carbon-to-DRI ratio was 0.099 kg/kg (a 17.5% reduction from the original case 0.12 kg/kg), obtained using H 2 /CO and (H 2 + CO)/(H 2 O + CO 2 ) ratios of 1.39 and 12.5, respectively. These values are in line but somewhat different from those of the Gilmore case, due to the content in other gas components of the bustle gas.

Discussion
The obtained results were already given a first interpretation in the previous sections. Hereafter, we focus on the context and significance of the present work. The developed process model differs from previously published works by its multi-scale nature, from grains (µm) to the shaft furnace (hm), and ultimately to integrated plant simulation (Aspen Plus).
The core of the DR process is the reduction shaft furnace. The corresponding Aspen Plus model was made sufficiently sophisticated to reproduce the main results from the more complex, CFD-type models of this reactor, but sufficiently light enough to be integrated in a whole plant flowsheet model. Thus, based on the recent results of [14], which proposed a virtual division of the shaft into zones distinguished according to the physico-chemical and thermal phenomena occurring, the shaft furnace was modeled as a set of interconnected zones. The ones in the reduction and intermediate sections were discretized in horizontal slices. This description, intermediate between 1D and 2D, enables us to retain the advantages of describing the axial variations of the calculated variables and of a short computation time (typically 3 h on a standard PC). The results from this model agree well with both available industrial data and results from previous models. Moreover, the two cases simulated corresponding to plants of different operating conditions and of quite different throughput, this demonstrates the adaptability and the robustness of the model.
The overall Aspen Plus model of a DR plant includes models for the principal units: the shaft furnace, the natural gas reformer, the scrubber, and the main heat exchanger. To the best of our knowledge, this is the first published work on a systems model of the whole DR plant using process simulation software. The decisive interest of such an approach is that the recycling gas can be accounted for. In MIDREX-type DR plants, most of the top gas exiting the shaft furnace is recycled, being first sent to the reformer then fed as the reducing gas at the gas inlet of the shaft furnace. Using the overall model, it becomes possible to study the interactions between the respective performances of the reformer and of the shaft furnace.
As a test seeking to mitigate the CO 2 emissions from the plant, the reducing gas composition at the shaft gas inlet was varied. A series of simulations show that CO 2 emissions and natural gas consumption could be reduced tuning the ratios H 2 /CO and (H 2 + CO)/(H 2 O + CO 2 ) at 1.23 and 12, respectively. The C-to-DRI ratio can be lowered from 0.119 to 0.105 kg/kg, a 12% reduction.
However, this optimization by hand is a first piece of work. The next stage would be to use a mathematical optimizer for the same purpose. The overall model could also be used to test different plant configurations. A variety of DR plants (MIDREX-type of different sizes, HYL-type of different sizes, HYL-ZR, ENERGIRON with different reducing gas sources) were designed and built. This means that further configurations could be designed to optimize criteria like the CO 2 emission, but also energy consumption, operating costs, etc. The present model can give useful answers.

Conclusions
In conclusion, this paper showcases a systemic, multiscale process model developed for the simulation, visualization, and further optimization of a gas-based DR iron reduction plant. A multi-scale approach was adopted, going from grain-size to model the iron pellet and its transformation, to reactor-size where reaction kinetics and heat and mass exchanges were simulated, and ultimately plant-scale where the gas recycling loop was modeled to obtain a closed process, and the subsequent carbon emissions were calculated for two cases. Simulation results were found in good agreement with industrial data. Moreover, hand optimization tests provided a 12% decrease in carbon emissions. Future works will address computer-based optimization along with the testing of novel process configurations. The calculation loop can be broken down as follows. All gas and solid streams are initialized as equal to the input gas and solid streams respectively. Calculation then starts from the first slice up to the last. Reaction and heat and mass exchanges occur between the ascending gas (i) and the descending solid (i + 1) entering the slice. This exchange impacts the exiting ascending gas (i + 1) and descending solid (i). The difference between the newly calculated values and the previous values is then recorded for each stage. Once the calculation is over, the overall relative error is calculated as the sum of the difference between the old and new value divided by the sum of the total values for each of the following variables: molar flow rate for gas and solid components as well as the temperatures for gas and solid phases. The iteration then continues until either this error is lower than a certain tolerance (10 −4 ) or the number of iterations is superior to a maximum number of iterations (7000). Once the calculation is over, results are returned to the Aspen Simulation. These include: the reaction extents, the exit solid and gas temperatures, as well as the temperature and composition profiles for both gas and solid streams along the shaft. Table A1 provides the list of all the reactions considered in the reduction and transition zones. The first set of equations relate to iron oxide reduction: Fe 2 O 3 → Fe 3 O 4 → Fe 0.95 O → Fe , with CO and H 2 as reduction agents. Their reaction rates are controlled by various resistances (inter-granular, intra-granular, inter-crystallite, intra-crystallite, and chemical) summed together courtesy of the additive reaction times law that considers that the different matter transportation steps occur in series. The rest of reactions are the methane steam reforming, water gas shift, methane decomposition and Boudouard. The kinetic rates are calculated based on stoichiometric factors with a possibility for reverse reactions. The expressions for the kinetic constants are based on literature [17,[21][22][23], with the role of iron as a catalyst for methane reactions. Reaction heat is calculated as the sum of each equation's heat of reaction and is attributed to the solid phase. Moreover, solid and gas phases transport heat by convection and conduction, and exchange heat through convective transfer.