Evaluation of the impact of cluster formation in a direct reduction shaft furnace through numerical simulation

The Direct Reduction (DR) process has been growing worldwide, and there are strong context suggestions that it will grow even more. One of these factors is the environmental pressure that occurs worldwide, and there are already projects to migrate Blast Furnace route steel plants to the Direct Reduction (DR) route, due to its smaller carbon footprint. Considering the importance of this process and the challenges of carrying out experimental tests on a pilot scale, an adequate way to evaluate the process and its impacts is through numerical simulations. There are different techniques applied to models that describe the counter-current reactor in the DR process, but none of them account for the clustering phenomenon. Clustering occurs because of the sintering of the metallic iron on the surface of the pellets in such a way that they attach to each other, forming clusters that hinder the gas flow through the shaft. The present study attempted to adapt a numerical model of a DR process to account for the effect of the cluster formation. Some clustering index equations from literature and some developed as part of this study were used and tested in the model, as a function of temperature, by varying the solid volume fraction in the control unit. The equation that resulted in the adjusted output closest to the current empirical value was implemented in the model and proved to be successful.


Metallurgy and materials
Metalurgia e materiais DRI (direct reduced iron) production in 2019 reached 108.1 Mt, representing an increase of 7.3% compared to the year 2018. It was the fourth consecutive year of increase. More than 75% of this DRI production comes from a reduction process in shaft furnaces (Midrex, 2020). Some countries such as India, Iran, Alge-ria, USA and Russia have stood out in this scenario due to new projects and plant reconditioning. In addition, the need for technologies with lower greenhouse gas emissions is a relevant driver for new investments in the future. It is estimated that a DR route with a 100% H 2 atmosphere may reduce emissions by approximately 80% compared to a conventional Blast Furnace route (Costa et al., 2013, Arens et al., 2017, Hille, 2018, Ahman et al., 2018, Dorndorf, 2018, Hybrit, 2018. However, the process has some challenges in its operation to guarantee dynamic process control and product quality. One of the challenges is the sintering of metallic iron formed at high temperatures and the associated reducing atmosphere causing the formation of clusters (Zhang et al, 2012, Yi et al, 2013, Battle, 2014, Alencar, 2016. This condition comes from when some raw material with a high susceptibility for sticking enters the shaft furnace. So, an increase in the bed pressure is noted, which tends to increase gradually.
The reduction process that takes place inside a direct reduction reactor is complex, as it involves a set of thermochemical reactions in a heterogeneous phase system. In this context, a huge number of models have already been developed trying to represent the DR process. These models range from singleparticle approaches to multiparticle and components models with topochemical and grain models (Venkateswaran and Brimacombe, 1977, Yu and Gillis, 1981, Takenaka and Kimura, 1986, Negri, 1991and 1995, Parisi and Laborde, 2004, Piotrowski et al., 2005, , Pineau et al., 2006, Thurnhofer, 2006, Valipour et al., 2006, Ajbar et al., 2011, Nouri, 2011, Costa et al., 2013, Shams and Moazeni, 2015, Kazemi, 2017, Castro, 2018, Rocha et al., 2019. Since the DR process can be summarized as a counter-current bed composed of solid and gas into a mass and heat transfer, the most acceptable and useful models correspond to numerical models that comprise the multi interaction and multiple phase theory. The insertion of elements and components that best describe the phenomena associated with the direct reduction process is a continuous target to achieve better operational predict-ability and quality values. A recent example that can be cited is the study of the particle shape factor charged in a Corex reactor (You et al., 2019). But, there are still no models that consider the phenomenon of clustering formation impacting the bed's permeability and its stability.
In this sense, the present study used data from other studies (Pereira, 2012, Alencar, 2016, Griscom et al., 2000 concerning the clustering index in different temperatures to provide enough data to establish parameters and correlations that could be added in a model based on transport equations applied to reactive multiphase and multicomponent systems in order to estimate the impact of cluster formation inside the bed on the operational parameters.

Modeling concept
The standard test to measure the susceptibility of clustering formation is the ISO11256. This test is based on an isothermal reduction at 850 °C in a reduc-ing atmosphere composed of CO and H 2 . At the end of the test, the reduced mass is weighed and put into a tumbling cycle. At the end of each revolution, the mass in the form of a cluster is weighed and subjected to the drum again (ISO, 2007). Eq.1 shows how the international standard ISO11256 standard establishes the clustering index.
The mathematical modeling of the DR process in the shaft furnace can be represented by the transport equation of momentum, energy and chemical species (Melaaen, 1992, Castro et al., 2018, Patankar, 2018, Rocha et al., 2019. In this approach, two phases can be discriminate, gas and solid. The gas phase is represented by the reductant gas injected in the reduction zone, the gas generated in the reduction reaction and reform, the cooling gas, and the carburization reaction. The solid phase is the metallic burden, that is composed of pellets and lump. The model can also represent the size distribution, divided into granulometric ranges for the pellets and lump. In this implementation, kinetic constants were assumed to be the same for all materials. Thus, the reduction rate is different due to physical characteristics, such as size, shape factor and porosity. The principle of solid arrangement in the control volume is assumed. Thus, this approach considers that the gas phase fills the voids. Therefore, as in the rule of continuity, then Equation 6 can be applied. CI1 (%) = 0.38 x T(°C) -308.67 (Alencar, 2015) CI2 (%) = 0.445 x T (°C) -376.75 (Alencar, 2015) CI3 (%) = 0.435 x T (°C) -347.25 (Alencar, 2015) CI4 (%) = 0.0004e0.0122 x T(°C) (Pereira, 2012) ε g + ε s = 1 Where: m r is the total mass after reduction and cm i is the mass (g) of the cluster portion after each revolution "i".
Previous studies (Pereira, 2012, Zhang et al., 2012, Yi et al., 2013, Alencar, 2015 and 2016) evaluated the metallic iron morphology present at different pellet clusters considering different reduction temperatures. According to these studies, as the reduction temperature increases, the clustering index also rises. Based on these data, it was possible to formulate Equations 2, 3, 4 and 5 that correlate with the clustering index (CI) as a function of temperature. Three of these equations consider the effect of temperature as linear and one considers the effect as exponential. All of them were evaluated in this study.

a. General conditions
(3) For each phase considered in the control volume, the conservation rule is applied. Thus the momentum, energy and chemical species transport equation can be calculated. To represent the industrial process, the corresponding boundary conditions are applied in these equations. The fluxes and volumetric composition to the reduction and cooling gas inlet are assumed to be known. For the gas phase, the reduction and cooling zone are considered. Due to pressure difference, the gas does not pass from one zone to another. The solids assume movement against the upward gas flow and it is allowed to pass through these two zones. In the wall region, the gas is non-slip, and for the interaction with the solid particles, a friction coefficient is adopted. To represent the cooling of the refractories, a heat exchange coefficient is applied to the walls. For the two phases, the flows are assumed to be completely developed at the outlet.
Equations 8 to 11 represent the phase momentum, mass conservation, energy, and chemical species, respectively.
Regarding phase interactions, the model considers the momentum, energy and mass changes. The momentum interactions are modeled by semi-empirical equations (Eq. 12 and 13) with adjusted constants to the shaft furnace process.
Equation 14 is applied to modeling the heat transfer between the solid phase components and the gas phase.

Momentum interaction:
In the above equations, the index i and j represent the phases and velocity components, and n represents the chemical reactions. The terms μ, ε, ρ are the dynamic viscosity volumetric fractions and phase densities, respectively. U and u are vectors and velocity components. P is the pressure and F are the exerted strength by other phases.
Cp, k and ∆H are heat capacity, thermal conductivity and heat due to chemical reactions, respectively. R is the reaction rates and M is the molecular weight of the species.
The indexes g and s indicate gas and solid respectively, m is the granulometric range of the solid phase component, f is the volume fraction of the granulometric range of the component in the solid phase, F is the interaction strength between the phases, d is the mean diameter of the granulometric range of the component in the solid phase and ϕ is the mean shape factor in the granulometric range of the component in the solid phase.

b. Process governing equations and boundary conditions c. Phase interactions d. Convection and radiation energy transfer
The term ε g in the above equation corresponds to the gas fraction and the ε s is the solid fraction in the volume. To calculate the solid fraction in the volume, the mixture rule is applied according to Equation 7.
The granulometric ranges for pellets and lumps in the solid phase are represented by the index "m" in the above equation.
The above equation is accountable for convective and radiation heat transfer among the phases and T is the phase temperature.

e. Phase Properties
Based on its components, the mixing rule is used to calculate the solid and gaseous phase properties. From the ideal gas law, we can determine the gas-phase density.
In the above equation, the subscript j means the gas species. The viscosities of these components are determined as follows (Bird et al., 1960, Reid et al., 1988.
The Eucken´s polyatomic gas approximation is used to determine the thermal conductivities for the components (Wilke, 1950, Neufeld et al., 1972. The wilke method (Wilke, 1950) is used to determine the viscosity and thermal code of the phase.
Where it is possible to determine the binary diffusivity of the gas species according to the equation below.
The terms of the above equation are determined as follows: The terms that appear in the equation are presented below: The above terms are: Boltzmann constant and σ and ∈ which are related to the collision gas types (Bird et al., 1960, Reid et al., 1988.
The equation below defines the gas temperature as a function of gas composition and enthalpy.
Considering the intra-bed radiation and the limit layer convection, the thermal conductivity is determined (Bird et al., 1960, Reid et al., 1988, Akiyama et al., 1992. The model uses the finite volume method (FVM) to discretize the equation of mass, energy, momentum and species taking into account a general coordinate system (Piotrowski et al., 2005, Pineau et al., 2006. These equations are integrated over a controlled volume; thus, a set of algebraic equations is created. The power-law scheme is used to determine the algebraic equation coefficients (Pineau et al., 2006). To improve the coupling of the pressure and velocity fields, the SIMPLE algorithm technique is used iteratively. To efficiently solve the algebraic equation system, the tridiagonal matrix algorithm based on the line by line solution is used. To improve the convergence, the alternated direction implicit is applied in the model solver. The solution is obtained iteratively to a stationary condition. The solution accuracy is reached by convergence criteria of 10-6 for the energy and momentum equation. The convergence target to the mass balance equation considers a difference of 0.01.
The term e s is the emissivity of solid components and the above constants are determined as follows:

f. Numerical solution
Based on the solid properties, the component's heat capacity are calculated as a function of temperature.
The empirical equations (Parisi and Laborde, 2004) for the solid fraction com-ponent ε m are actualized according to the new diameter in Equation 38: Figure 1 schematically shows the configuration of the control volumes in different aspects. Scenario A shows a reference situation with the pellets' average diameter (d m ). Scenario B demonstrates the real aspect of the pellets that, under reduction and high temperatures, faces softening, deformations, and the sticking phenomena, leading to the gas flow disturbance. On the other hand, scenario C is the representation of the ap-proach given in the present study, where the smaller diameter of the pellets generates a greater packing factor and, consequently, greater pressure drop, resulting in the same practical effect observed in scenario B.
The evaluation of the clustering equations with different temperatures was validated, aiming similar results that were achieved in two industrial cases. In both scenarios, the reactor was a shaft furnace type and the impact of the clustering index was converted into relative productivity. After comparing and defining the best prediction equation for CI and its impact on solid volume fraction, new temperature scenarios were tested to evaluate the sensibility related to it.
Before starting the tests involving the new elements in the model, such as the equation that correlates solids fraction as a function of the clustering index, an initial validation of the model was carried out based on a shaft furnace's industrial data from a mill located in MENA. Tables 1 and 2 show, respec-tively, the main operational data of this industrial reactor and the predominant gas composition in the reduction and cooling zones.   Taking into account the tests and results mentioned above where it is possible to obtain a correlation between temperature and clustering index, the model applied the equations to consider the presence of pellet clusters inside the DR reactor.

Validation of the reference model with industrial data
All equations formulated were applied in the model considering the increase of the solid fraction in the control volume and consequently decreasing the gas perme-ability to the region. In order to increase the solid volume fraction, the mean diameter d m of each solid component is recalculated as a function of CI (Eq. 37), where CI is the Clustering Index in decimal format. It can be seen in Table 1 that the operational practice of the selected reactor is in line with the typical modus operandi of DR plants, where the raw material pellet enters the reactor at room temperature, and the bustle gas temperature is around 1000 °C (Atsushi, 2010). In addition, Table 2 shows that the reduction gas composition has an H 2 /CO ratio of approximately 1.85, which is a value that guarantees good conditions for reduction kinetics and thermal balance in the system (Capriotti, 2012). Finally, the cooling gas is very rich in CH 4 to favor the carburization and self-reforming process inside the reactor.

Pellets clustering modeling
The pellet burden used in this validation step also represented the market share consumed by this mill studied. Three pellets from different suppliers were simulated according to chemical and granulometric quality highlighted in Tables 3 and 4. According to Table 3, all pellets in the mix have similar iron content. However, the acid gangue content is lower in pellet C, which also has the highest values of CaO and MgO. Meanwhile pellets A and B are very similar in general, but pellet B has higher binary basicity and higher phosphorus content. In terms of size distribution, Table 4 presents in detail the particle size fractions of each pellet tested in the simulation. These data are one of the model inputs and are used to calculate the average particle size and consequently the bed porosity. Through Table  4, it is possible to say that pellet A has the distribution of larger pellets, while pellet C has the distribution of smaller pellets. The most used pellet in the reactor is A (50%), followed by C (40%) and B (10%). Table 5 summarizes the simulation comparing the results obtained from the model based on transport equations, multiphase and multicomponent systems with the industrial results provided by the studied mill. It appears that the key process parameters, such as metallization, carburation and production in the model were comparable to the industrial values collected. There are no values in the remaining parameters similar to the industrial averages, but were within the confidence interval of the industrial values measured over 30 days. Figure 2 shows the temperature profile of solids and gas. The results obtained in the simulation are aligned with what is expected for an typical industrial operation. First, the temperatures of solids and gases are close to each other and the highest temperature observed inside the reactor in both cases is slightly lower than the bustle gas temperature because  after entering the shaft, the bustle gas expands its volume and subsequently experiences a slight pressure drop. In addition, the regions of higher temperatures in the reduction zone are close to the walls, whereas the cooling zone presents lower temperatures due to the carburization reaction with CH 4 , which is endothermic. In all tests performed, each of these equations generated a CI that replaced the Eq.
15 and led to different outputs of the solid fraction and consequently permeability. The CIs that result from each equation for different temperature ranges are shown in Table 6. All cases/equations had a version called "Adjusted" that represents the return to the base case's original condition for Dp (pressure drop). This understanding reflects the operational practices that always aim to work within an ideal range of Dp in order to keep operational safety and stability. Thus, all four evaluation cases use equations that correct the solid fraction volume, which decreases the permeability and in-creases the Dp of the reactor. Therefore, in the adjusted cases, the Dp is reestablished, correcting gas flow and/or temperature, and thus these adjusted cases lead to lower productivity, as shown in Table 7.
Such numerical simulation results in Table 7 were contrasted with the practical experience based on two real situations involving two DR shaft furnace reactors from different locations.
In one case, a pellet consumes a CI of 10%, which was replaced by another pellet of the same supplier, but with a CI of 39%. In the other case, the same pellet with a 10% Clustering Index was changed by one with a 13% CI. The impacts in both scenarios were similar. It was necessary to reduce the temperature and flow rate of the bustle gas over a period and, as a result, the furnace pro-  Table 7 -DRI and process parameters for each case evaluated.

CI approach results
ductivity decreased. Even though they are different furnaces, dealing with the data in both cases led to a correlation, where each additional 1% CI in the bur-den would result in a relative decrease of production around 0.83%. Table 8 presents the data from each numerical simulation based on the four equations formulated. Figure 3 shows the decrease of production-related to each 1% of CI, aiming to highlight the best value compared to the empirical rule. It can be seen from Table 8 and Figure 3 that the Equations CI1 and CI4 are the ones that come closest to the industrial value of 0.83%. However, when observing Table 6, it is noted that extrapolating the Equation CI4 to temperature values just above 1000 °C would result in CI values above 100%, which is some-thing impossible. Therefore, the most suitable equation considered to represent the impact of the clustering inside the reactor was the Equation CI1. Figure 4 depicts the influence of the bustle gas temperature on CI. Each slice of the furnace represented by the letters in Figure 4 refers to a simulation made at a specific bustle gas temperature in its steady state. As can be seen, the fraction of solids (clusters) increases with increasing bustle gas temperature, where Equation Cl1 was applied to conduct this investigation. This aspect is in line with the industrial reports of buildup problems (Costa et al., 2013, Voelker, 2019 because clusters usually concentrate. In more severe cases, these clusters can accumulate inside the furnace and drastically obstruct the permeability, requiring a maintenance cold stop. This is why reactors cannot be found operating at bustle gas temperatures much higher than 1000°C.

References
The experience of seeking literary references about the temperature impact on the formation of the clusters and associating these equations with the solid's fraction volume and their respective bed permeability was novel and successfully tested using an adjusted mathematic approach. Based on this study, new evaluations of potential coating agents may be performed more assertively through labo-ratory tests and the developed numerical model, saving costs and time with more complex and industrial tests. In addition, from the data obtained with this study, it can be stated that: • The adaptation of a pre-existing model based on the finite volumes method focusing on the RD process was suitable to assess the impact of the cluster formation inside a DR shaft furnace reactor.
• From the equations evaluated in this study, the Equation CI1 (CI1 (%) = 0.38 x T (°C) -308.67) was the one that showed the best adherence to the empirical result of real case studies.
• This equations for clustering modelling used in different temperatures presented the same industrial trend, starting with the reactor walls and going towards the center.