Process Modeling and Optimization of a Submerged Arc Furnace for Phosphorus Production

Abstract

This article presents a process model of a phosphorus-producing, submerged arc furnace. The model successfully incorporates accurate, multifield thermodynamic, kinetic, and industrial data with computational flow dynamic calculations and thus further unifies the sciences of kinetics and equilibrium thermodynamics. The model is structurally three-dimensional and uses boundary conditions, initial values, and material specifications provided by industrial measurements, laboratory experiments, and a combination of empirical and thermodynamic data. It accounts for fully developed gas flows of gaseous product generated from within the packed bed; the energy associated with chemical reactions, heating, and melting, as well as thermal conductivity and the particle–particle radiation within the burden. The model proves the existence of a narrow, gas–solid reduction zone where the bulk of phosphorus is produced. It shows that fast reaction rates in this narrow reaction zone in combination with long residence times diminish the influence changing reaction rates have on the process. It indicates that most heat exchanged between the new pellets entering the furnace and the gaseous product produced in the reduction zone takes place in the top 0.5 m of the furnace bed. The gaseous product and flow information shows low and recirculating gaseous flow velocity areas that cause dust accumulation.

Introduction

A phosphorus-producing, submerged arc furnace (SAF) is a complex system in both physical and chemical aspects. It involves multiphase, high-temperature reduction reactions, energy conversion, and distribution from electric power through arcs and conduction. Although the underlying chemical reactions and, to a lesser extent, the kinetics driving the main reactions in the submerged arc furnace are well documented, it is the complex interaction between these phenomena and the intricate burden characteristics that cause the specific power consumption (the energy consumed for every ton of P₄ produced) to be a largely uncontrollable, dependent variable.

The most complex electric powered furnace is the SAF, and over the years, it has proven difficult to model. Larsen et al.[1] developed a numerical two-dimensional (2D) model for an alternating current (AC) arc in a silicon metal SAF using Fluent. The arc is the main energy source in the furnace (≈90 pct). The conservation equations for mass, momentum, and energy together with time-dependent Maxwell’s equations were solved. The model assumes symmetric furnace conditions and takes only one phase in consideration for a three-phase SAF. Therefore, the significant interactions among phases cannot be reflected by the model. Andresen and Tuset[2] simulated fluid flow, heat transfer, and heterogeneous chemical reactions above 2123 K (1850 °C) in the arc region of a silicon furnace with a 2D Fluent model. The AC arc in the gas-filled cavity is simplified with a direct current (DC) arc. Saearsdottir et al.[3] investigated arc behavior in silicon and ferrosilicon furnaces using a 2D model. The behavior of the solid bed region and overall furnace thermal performance was not included. Sridhar and Lahiri[4] developed a 2D single electrode model for the current and temperature distribution in a SAF for ferromanganese production. The current distribution is calculated by solving Maxwell equations for the magnetic field. Joule heating and heat conduction are modeled in the slag and solid zones, and the arc and gas flows are not included. Heat sources as a result of chemical reactions are included. The influence of electrode depth, electrical conductivity of the charge, and the slag on power consumption was studied. Ranganathan and Godiwalla[5] used the same approach for a ferrochromium SAF to simulate the effects of bed porosity, charge preheat, and charge chemistry on the temperature profile. Vanderstaay et al.[6] modeled the submerged arc electric furnace smelting using the manganese-rich slag method process by using a computational thermodynamics package. It was assumed that higher manganese and iron oxides are reduced to MnO and FeO before entering the zone where molten slags and alloys form and equilibrate. The model predictions were compared with data from Thermit Alloys (P) Limited (Maharasthra, India), an Indian ferroalloy smelter. It then was used to examine the effects of changing the amount of carbon reductant and temperature on several performance indicators. No attempt has been made so far to develop a three-dimensional (3D) model of a complete furnace to simulate chemical reactions in the packed bed combined with gas generation and flow. This prompted the authors to develop the model set out in the article.

The Phosphorus Producing Process

The main process reaction as defined by the Wohler process[7] proceeds according to Eq. [1]. From the Moeller feed consisting of pelletized (fluoro)apatite (Ca₁₀(PO₄)₆F₂), coke, and silica (in the form of gravel), it produces a calcium-silicate slag, calcium fluoride, carbon monoxide and the desired product, phosphorus gas as follows:

$$ {\text{Ca}}_{10} \left( {{\text{PO}}_{4} } \right)_{6} {\text{F}}_{2} + 15{\text{C}} + 9{\text{SiO}}_{2} \to 3{\text{P}}_{{2({\text{g}})}} + 9\left[ {\left( {{\text{CaO}} {\cdot} {\text{SiO}}_{2} } \right)} \right] + {\text{CaF}}_{2} + 15{\text{CO}}_{{({\text{g}})}} $$

(1)

A 2D depiction of the SAF is shown in Figure 1.[8] Gravity delivers the feed to the SAF through 10 evenly distributed feed chutes ensuring a constant packed bed volume. The gaseous product leaves the furnace through two symmetrically spaced outlet vents situated above the thin alloy layer of ferrophosphorus tap hole in the roof of the furnace. The ferrophosphorus is tapped off, usually once per day. The specific power consumption is the most symptomatic indication of the profitability of the process and the optimization thereof is the first priority at any phosphorus producing company. Molten slag, however, continuously is tapped alternately through two water-cooled tap holes located 400 mm above the furnace floor. Because of the large production volume of slag, a seemingly low content of P₂O₅ in the slag will result in substantial losses of unreduced, potential product. The second priority is to control the process to keep the P₂O₅ in the slag as low as possible.

Fig. 1

A 2D depiction of a SAF.[8] The constant packed bed volume, three feed chutes, electrodes, gaseous product space and slag, and ferrophosphorus alloy layers are shown. The product gas leaves the furnace through two outlet ducts shown in Figure 3

The feed material forms the major electrical and flow resistance of the smelting furnace circuit. As the feed materials descend toward the hot zone in the furnace, they start to soften and melt, significantly lowering the electrical resistance. A conductive path is provided between the electrodes where the Joule heating (I²R) is released to attain the high temperatures and energy levels necessary to effect the endothermic reactions. These reactions take place in a dome-shaped area or a reduction zone around the tips of the electrodes. The shape, size, and reach of this reduction zone through the packed bed are fundamental to understanding and controlling the SAF. The continuous tapping of the slag, the insignificant heel associated with the process, and the fluctuating state of the furnace[9] make the slag–metal reactions within this particular furnace not negligible but certainly less influential than found by Eksteen[10] for ferrochrome and ferromanganese furnaces. Work by Dresen et al.[8] and Van der Pas[11] revealed that the solid–gas reactions inside the packed bed and not the liquid–gas reactions inside the slag are responsible for most production of the gaseous product. It is also inside the packed-bed region where the rate-limiting steps of the process are located.[12] It is for these reasons that the focus of the work in this project is on the packed-bed zone and on the gas–solid reactions.

This study also addresses a third priority of understanding and the eventual minimization of the varying energy consumptions created through changing mixtures of primary feed materials as a result of ever-increasing turbulent ore markets. By generating a virtual window to the inside of the furnace, operators can get an idea of the changing conditions inside the furnace as a result of these changing feed materials.

General Model Construction

A literature survey indicated that no attempt has been made so far to develop a 3D model of a complete furnace to simulate chemical reactions in the packed bed combined with gas generation and flow.

Model Characteristics

The presented model of a SAF has the following unique characteristics:

• (a) Being 3D

• (b) Using a furnace modeling structure that is exactly the same as the actual furnace

• (c) Using boundary conditions, initial values, and material specifications provided by industrial measurements, laboratory experiments, and a combination of empirical and thermodynamic data

• (d) Accounting for fully developed gas flows of gaseous product not introduced to but generated from within the packed bed through the Wohler reaction

• (e) Accounting for the energy associated with chemical reactions, heating and melting, and phase transformations in the packed bed

• (f) Accounting for the thermal conductivity and the particle–particle radiation within the packed bed

• (g) Implementing the integration of accurate, multifield thermodynamic data and kinetic reaction rate data with computational fluid dynamic calculations

The model does not simulate Maxwell equations and Joule heating, and the liquid phase of molten slag and alloy phase are modeled as stagnant liquid phases.

Operation Data of the Industrial Furnace

To develop a reliable model for the SAFs, 6 months of actual process data from a production plant were collected, filtered, and reconciled. A detailed description of these procedures is reported elsewhere.[13] The 8-hour averages then were arranged in order of ascending Moeller feed (ton/8 h), using the five ranges shown in Table I. Of these ranges, five main scenarios representing the averaged data of that particular range were developed. It is on these five scenarios and derivatives thereof that the model development is centered. Table I is a depiction of the averaged values that each of the five scenarios represent. The energy inputs are averaged during an 8-hour period. The corresponding pellet, silica gravel, and coke flow rates also are averaged values for the same period, whereas the gaseous outlet temperature value is averaged for both outlet pipes. The P₂O_{5 (slag)} value is the P₂O₅ content of the slag from a single 8-hour sample.

Table I Averaged Feed Flow Rates, Power Consumption, Gaseous Outlet Temperatures and P₂O_{5 (slag)} Content in the Slag of the Five Scenarios Derived from Actual Process Data

From Table I, the three flow rate and energy input values are used as (some of) the input parameters to the model, whereas the average gaseous product outlet temperatures and the P₂O₅ content of the slag are used for validation purposes. For all five scenarios in Table I, a fixed coke-to-P₂O₅ addition (C-fix) of 0.48 (standard deviation = 0.012) was used. The averaged chemical analyses of the Moeller feed (apatite pellets, coke, and silica gravel) of all scenarios are shown in Table II.

Table II The Relative Concentration of the Moeller Feed Material to the Furnace and Subsequently to the Furnace Model. This Is the Format Used as Input to Subsequent Factsage[14] Work Performed to Assist with the Model Development

In Table I, scenario 3 represents the average, steady-state conditions at the production plant in question and is referred to as the base case scenario, which was used for initial model development.

Furnace Dimensions, Structure and Computational Grid

The furnace model was constructed with Fluent version 6.1.18,[15] and the mesh was generated in Gambit version 2.04.[16] A hexagonal mesh scheme (Cooper type) was used for the furnace and both a tetrahedral/hybrid mesh scheme (T-grid type) and hexagonal mesh scheme were used for the outlet ducts. The total number of cells for the base case scenario is 414,000, and the typical solving times were between 12 and 20 hours. A cut-through figure of the 3D constructed furnace along with selected parts of the grid is shown in Figure 2.

Fig. 2

A 3D depiction of the SAF’s dimensional grid as constructed in Gambit.[16] The inlet zones, the electrodes and the slag-packed bad interface are displayed. Parts of the furnace structure were deliberately left out in order to depict the inside of the furnace

The entire body of the furnace lining was constructed as a conducting solid using the individual manufacturer specifications of each type of lining (Table III).

Table III Physical Properties of the Concrete Top, the Chamotte, and the Carbon Lining and the Three Electrodes of the Furnace

Figure 3 shows carbon bricks in the bottom and lower part of the furnace wall, chamotte bricks in the upper part of the wall, and concrete for the furnace lid on top as well as the feed chutes. The chamotte bricks are made from sintered clay and consist mainly of SiO₂ and Al₂O₃. The electrodes are manufactured from high-quality carbon. The outlet duct has a zero thickness and virtually was extended to establish a fully developed flow and at the same time minimize the effect of reversed flow. Below the packed bed, the slag and the ferrophosphorus phases are simplified as nonflowing, liquid layers but still retain all material properties associated with the respective phases.

Fig. 3

A 3D depiction of the SAF as constructed in Gambit.[16] Four of the five boundary conditions introduced to the model (Table VII) are shown with the dotted lines

Gu and Irons[17] employed a channel arc model that showed how the arc from an electrode could be approximated by a cylinder. For this reason, three cylindrical inlet zones created between the electrode tips and the slag surface and the diameter of the electrodes were used as the diameter of the inlet zones. The volumes of these zones vary depending on the operating height of the electrodes and provide the model with individual energy input values unique to each electrode. For the base case scenario, each inlet volume is 1.31 m³ in volume. The packed bed inside the furnace is represented with the Fluent porous media model as well as with four user-developed submodels. These submodels are discussed later. The formation of freeze lining was not taken into consideration.

Reactant and Product Material Properties

The density and porosity values of the apatite pellets and coke particles were measured by using a mercury intrusion porosimetry technique. These values are shown in Table IV.

Table IV The Density and Porosity Values of Apatite Pellets, Coke, and Silica Gravel

Given a weight ratio in the Moeller of 68.9 pct apatite pellets, 10.7 pct coke, and 20.4 pct silica gravel, the Moeller bulk density (ρ _bulk) equals 1843 kg/m³. The thermal conductivity and radiation aspects of the Moeller are introduced in combination with the user-developed model construction. The slag density was calculated using the slag composition from the base case scenario and the individual component density data from HSC Chemistry version 5.0.[18] The slag thermal conductivity was provided by The Slag Atlas.[19] The slag heat capacity was determined experimentally at the industrial plant laboratory and validated through HSC Chemistry. The gaseous product density is calculated by Fluent with the ideal gas law. The gaseous product thermal conductivity above the packed bed (but still inside the furnace) was determined by using data from Yaws.[20] The gaseous product heat capacity was calculated for CO and P₄ gas by using HSC Chemistry. The ferrophosphorus density and heat capacity were obtained from industrial data. These values are summarized in Table V.

Table V The Density, Thermal Conductivity, and Heat Capacity Values of Slag, Ferrophosphorus, and Gaseous Product

Input Conditions

The two energy inputs to the furnace model are summarized in Table VI. Figure 4 shows a graphic depiction of the energy sources.

Table VI A Summary of the Two Energy Input Conditions to the Furnace Model

Fig. 4

A 2D depiction of the two energy input conditions as introduced to the model (Table VI)

Gravity delivers the new Moeller feed continuously at approximately 573 K (300 °C) to the top of the packed bed in the furnace. To maintain this temperature at the top of the bed, a negative energy input is employed over the top 0.25 m of the packed bed domain as follows:

$$ Q_{\text{newfeed}} = \left[ {\int_{{313\,{\text{K}}}}^{\text{Tcell}} {Cp_{\text{pb}} } \times dT - \int_{{313\,{\text{K}}}}^{{573{\text{T}}}} {Cp_{\text{pb}} } \times dT} \right] \times W \times \rho_{\text{bulk}} \times \left[ {{\frac{1}{{z_{\text{cell}} }}}} \right] $$

(2)

where Cp _pb is the heat capacity of the packed bed and W is the rate of descent in the downward direction of the actual packed bed at the production plant m/s. z _cell is the height m of the control volume within the grid and ρ _bulk is the bulk density of the Moeller bed.

Boundary Conditions

The five boundary conditions defined for the furnace model are summarized in Table VII.

Table VII A Summary of Five Boundary Conditions Defined for the Furnace Model. Four of These Boundary Conditions (1 to 4) Are Shown in Figure 3

Exposed Furnace Surface

A large part of the stainless steel surface of the furnace is exposed to the ambient air (the exposed surface in Figure 3) and has a measured surface temperature of ≈373 K (100 °C) during steady-state operation. The average radiative heat loss over this oxidized steel surface is ≈0.8 pct of the total energy consumption. In the model, radiative heat losses were accounted for by introducing a mixed regime across the relevant surface that included radiation and convective heat losses. The parameters used were radiation ambient temperature = 623 K (350 °C), emissivity = 0.8 for oxidized steel,[21] convection ambient temperature = 573 K (300 °C), and convection heat transfer coefficient = 10 W/m² × K.

Heating the Moeller from 313 K to 573 K (40 °C to 300 °C)

The part of the feed chutes not protruding into the furnace is more than 20 m in length as it extends upward. The average temperature of the Moeller as it is fed to the furnace at the top of these external feed chutes is 313 K (40 °C). By the time the Moeller enters the furnace, the temperature has increased to approximately 573 K (300 °C). However, these feed chutes do not form part of the domain. To reflect reality, the model domain was treated as a quasiexternal heat transfer problem to ensure that the Moeller temperature, as it enters the furnace, can be initialized safely at 573 K (300 °C). This energy is provided by the energy input 1 and is calculated as follows:

$$ Q_{\text{feedchutes}} = \dot{m}_{\text{moeller}} \times \int_{{313\,{\text{K}}}}^{{573\,{\text{K}}}} {Cp_{\text{pb}} \times dT} $$

(3)

It is accounted for by introducing the appropriate negative heat flux over the same area as the electrode cooling water flux. In the base case scenario, this amounts to −2.56 MW, which brings the total negative flux introduced to the top of the electrodes to (−5 MW) + (−2.56 MW) = −7.56 MW.

Mass Flow

No mass flow boundary condition needed to be defined because the mass is generated from within the packed bed in the reaction zone of the furnace.

Physical Models

Turbulence Model

In this study, the standard k − ε model was used to model the turbulence. The k − ε model is a semiempirical model based on model transport equations for the turbulence kinetic energy k and its dissipation rate ε and assumes fully turbulent flow. It is robust, economical, and reasonably accurate. Multiple simulations showed that the k − ε model (with default settings provided by Fluent) solved the turbulence effectively, whereas more advanced models like the renormalization-group k − ε model, the Spalart–Allmaras model, and the Reynolds stress model did not provide any significant improvements.

Radiation Model

In this study, the P1 radiation model is used to model the radiation above the packed bed. Radiation contributes greatly to the heat transfer in this gas zone. Inside the packed bed, a submodel was developed specifically to simulate the radiative effects. The P1 radiation model is based on the expansion of the radiation intensity I into an orthogonal series of spherical harmonics and is the simplest case of the more general P_n-approximations model. Multiple simulations showed that the P1 model (with default settings provided by Fluent) solved the radiation effectively. Among the other models, the discrete transfer radiation model (DTRM) does not include scattering, whereas the discrete ordinate model is too time consuming and does not improve the solution significantly.

User-Developed Model Construction

The furnace model does not have a fixed mass input such as normal computational fluid dynamics inlet boundaries. The mass input is the gaseous product (a mixture of P₂ and CO) originating from the chemical reaction that is generated within the packed bed. The gaseous product is generated as a result of the interaction between the chemical reaction energy, heating and melting energy of the packed bed material, as well as unique particle-to-particle radiation and thermal conductivity phenomena. These complex relationships were facilitated through the development of four user-developed models.

User-Developed Model 1: Mass Source Generation

Defining a Scalar Value for P₂O₅ Concentration

A scalar was defined in the model that represented the concentration of unreduced P₂O₅ in the pellets. An initialization value of 230 kg of P₂O₅ per cubic meter corresponds to 29.1 wt pct P₂O₅ in the pellets and a packed-bed porosity of 37.7 pct as shown in the base case scenario. This initialization value is a function of the type of apatite ore, the feed composition, bulk density, particle size distribution, and individual particle sizes.

Formation of a New Transport Equation for the Scalar (P₂O ^S₅ )

The next step is to simulate downward packed bed movement. In this way, P₂O ^S₅ can be transported down through the furnace bed. To make this transfer possible, Fluent can solve the transport equation for a user-defined scalar in the same way that it solves the transport equation for a scalar such as mass fractions of species. Mathematically, this process can be described as the formation of an additional advection term within the Fluent domain. This advection term in the differential transport equation has the following general form:

$$ \nabla \times \tilde{\psi }\phi $$

(4)

where ϕ represents P₂O ^S₅ . In the default advection term, $\tilde{\psi }$ is the product of a scalar density (ρ) and a velocity vector $\left( {\tilde{\nu }} \right)$:

$$ \tilde{\psi }_{\text{default}} = \rho \times \tilde{\nu } $$

(5)

where the velocity vector can be specified in three dimensions U, V, or W. Because of the packed bed moving in a downward direction (radial movement is assumed negligible), W has been defined as the downward velocity of this advection term; in other words it is the rate of descent of the actual packed bed at the production plant. For the base case scenario, the downward velocity W = −1.538 × 10⁻⁴ m/s, thereby modeling solid flow as plug flow.

P₂O₅ Consumption and Gaseous Product Generation

Throughout the iterative solution process of the model, the unreduced P₂O₅ left in the pellets $\left( {{\text{P}}_{ 2} {\text{O}}_{{ 5\, ( {\text{new)}}}}^{\text{S}} } \right)$ is used to calculate the rate of change of P₂O₅ concentration in the packed bed as the reaction takes place. For these calculations, the following Eq. [6] was used:

$$ \begin{gathered} 2{\text{Ca}}_{3} \left( {{\text{PO}}_{4} } \right)_{2} + 10{\text{C}} + 6{\text{SiO}}_{2} \to {\text{P}}_{{4\left( {\text{g}} \right)}} \hfill \\ + 6\left( {{\text{CaO}} \times {\text{SiO}}_{2} } \right) + 10{\text{CO}}_{{\left( {\text{g}} \right)}} \hfill \\ \end{gathered} $$

(6)

The onset temperature of ≈1423 K (1150 °C) for the P₂O₅ reduction was determined experimentally, and the reaction also is assumed to take place in the following specific steps[9]:

$$ \begin{gathered} {\text{Step}}\; 1 :\,{\text{Liberation}} \hfill \\ {\text{Ca}}_{ 3} \left( {{\text{PO}}_{ 4} } \right)_{ 2} + 3 {\text{SiO}}_{ 2} \to 3\left( {{\text{CaO}} \times {\text{SiO}}_{ 2} } \right) + {\text{P}}_{ 2} {\text{O}}_{{ 5 ( {\text{g)}}}} \hfill \\ \end{gathered} $$

(7)

$$ \begin{gathered} {\text{Step}}\; 2 :\,{\text{Diffusion}} \hfill \\ 2{\text{P}}_{2} {\text{O}}_{{ 5 ( {\text{g)}}}} + 10{\text{C}} \to {\text{P}}_{{4({\text{g}})}} + 10{\text{CO}}_{{ ( {\text{g)}}}} \hfill \\ \end{gathered} $$

(8)

Equation [6] is a first order solid–solid and solid–gas reaction, which can be modeled according to a shrinking core model.[12] In the shrinking core model, the P₂O₅ gas in Eq. [7] is liberated from a reaction surface moving from the outside to the inside of the pellet and is controlled by the removal rate of gases from the reaction surface. The reaction constants, k _rr 1/s, for a variety of apatite feed ores are calculated from empirical equations obtained experimentally through the kinetic investigation of an actual feed sample provided by the production plant.[11] Equation [9] shows the experimentally determined reaction constant equation for X-type ore. The use of actual ore highlights the engineering nature of the study. Equation [9] is written as follows:

$$ k_{\text{rr}} = 0.003 \times {\text{T}} - 0.4 $$

(9)

The rate of change of the P₂O₅ as a result of Eq. [7] now is calculated through the following Eq. [10] where t = 1 second, and k is the reaction constant of Eq. [9]:

$$ \Updelta {\text{P}}_{2} {\text{O}}_{5}^{\text{S}} = - {\text{k}}_{\text{rr}} \,{\cdot}\, {\text{P}}_{2} {\text{O}}_{{5\,({\text{old}})}}^{\text{S}} {\cdot} {\text{e}}^{{ - {\text{kt}}}} \left[ {{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {60}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${60}$}}} \right] $$

(10)

The remaining amount of P₂O₅ $\left( {{\text{P}}_{ 2} {\text{O}}_{{5\left( {\text{new}} \right)}}^{\text{S}} } \right)$ then is given by Eq. [11] as follows:

$$ {\text{P}}_{2} {\text{O}}_{{5({\text{new}})}}^{\text{S}} = {\text{P}}_{2} {\text{O}}_{{5({\text{old}})}}^{\text{S}} + \Updelta {\text{P}}_{2} {\text{O}}_{5}^{\text{S}} $$

(11)

Equation [12] initializes a new P₂O₅ value for the next iteration in the solving algorithm of the model as follows:

$$ {\text{P}}_{2} {\text{O}}_{{5\left( {\text{old}} \right)}}^{\text{S}} = {\text{P}}_{2} {\text{O}}_{{5\left( {\text{new}} \right)}}^{\text{S}} $$

(12)

Equations [10] through [12] are embedded in a coded loop structure within the model. $\Updelta {\text{P}}_{2} {\text{O}}_{5}^{\text{S}}$ is the amount of P₂O₅ that reacts according to Eq. [6]. This amount of P₂O₅ decrease in the packed bed is used to determine the volumetric mass generation rate of the gaseous product. From Eq. [6], it is determined that each kilogram of P₂O₅ produces 0.44 kg of P₄ and 1 kg of CO gas. The amount of P₄ and CO generated then is introduced into the model, and the amount of P₄ and CO generated is used in the determination of the energy distribution within the packed bed.

Additional Mass Sources

One additional aspect to consider is the additional coke fed to the furnace. This feeding is done (1) to prevent the carbon bricks from reacting, (2) to compensate for the loss of coke during slag tapping, and (3) to compensate for the reduction of other oxides like Fe₂O₃ and MnO. Additional reactions involving the extra coke develop an excess amount of gas leaving the furnace that is unconnected to Eq. [6]. This extra gas is introduced artificially to the furnace model at the three inlet volumes (zones) underneath the electrodes as shown in Figure 3.

Energy Distribution Ratio

It is important to note that the distribution ratio of the input energy in the furnace between the energy required for chemical reactions and the energy for heating and melting of the feed is not specified when solving the model. These values are generated entirely as dependent variables and are essential for a more realistic prediction of the furnace behavior.

User-Developed Model 2: Chemical Reaction Energy

For the development of the energy consumption part of the model, chemical reaction and equilibrium software Factsage version 5.4.1[14] was used in combination with the relative concentration of the Moeller feed material to the furnace (Table II). These individual component flow rates were recalculated based on an input of 1 kg of P₂O₅ per second and were solved with Factsage. The total energy needed as a result of gaseous product formation and phase transformation H _R is calculated as a function of temperature and is depicted in Figure 5. In this way, not only the energy associated with the primary reduction reaction (Eq. [6]) but also with all secondary reactions (which include alkali reduction reactions) is accounted for.

Fig. 5

The energy required for phase transformation, the energy required for subsequent gaseous product and slag formation, as well as the theoretical amount of product (given an initial amount of 1 kg P₂O₅) as calculated by Factsage

The 12 underlined values indicated in Figure 5 are used to construct a trend line from which the energy required for gaseous product formation alone at a specific temperature can be obtained. These values then are multiplied with the mass loss rate of the P₂O₅ value within each individual cell contained in the packed-bed reactor domain to obtain the total energy sink to be employed over the same domain. The calculation is shown in the following Eq. [13]:

$$ Q = - \left( {{\frac{{H_{\text{R}} \times \Updelta {\text{P}}_{ 2} {\text{O}}_{ 5}^{\text{S}} }}{{V_{\text{cell}} }}}} \right) {\frac{W}{{m^{3} }}} $$

(13)

User-Developed Model 3: Heating and Melting Energy

By using data from Table II combined with Factsage, heat capacity values for the packed bed (Cp _pb) were calculated for temperatures between 373 K and 3373 K. The raw Cp data with corresponding trend lines are presented in Figure 6.

Fig. 6

A graphical depiction of the heat capacity values for the packed bed as a function of temperature. Included are the trend lines with which the Cp values are represented within the coded structure of the model. The data from Table II were used to generate the values. The four regions marked from 1 to 4 are areas of large Cp value changes

In the same figure, four regions of sharp changes in the Cp value are identified in Factsage calculation. These changes are explained as follows:

• (a) Between 1317 K (1044 °C) and 1318 K (1045 °C), Ca₁₀(PO₄)₆F₂ → 3Ca₃(PO₄)_{2(whitlockite)} + CaF₂.

• (b) Between 1372 K (1099 °C) and 1373 K (1100 °C), 3Ca₃(PO₄)_{2(whitlockite)} → Ca₃(PO₄)₂.

• (c) From 1423 K (1150 °C), the main reaction starts taking place, and the combined Cp of the product and the reactant still increases.

• (d) From ≈1540 K (1267 °C), the main reaction proceeds, but now the combined Cp of the product and the reactant starts decreasing.

With temperature-dependent Cp values available, an energy sink as a result of heating and melting is employed over each individual cell in the packed-bed domain through Eq. [14] within the coded structure of the model. The energy required for heating and melting of the feed therefore is compensated for within every model iteration. Equation [14] is written as follows:

$$ Q_{\text{HM}} = - \left( {Cp_{\text{pb}} \times \rho_{\text{bulk}} \times \Updelta T_{{{\text{cell}}\left( {\text{W}} \right)}} \times W} \right)\quad {\frac{W}{{m^{3} }}} $$

(14)

User-Developed Model 4: Particle–Particle Radiation and Effective Thermal Conductivity Model

Because the porous media model of Fluent calculates the thermal conductivity of the packed bed as an average of the thermal conductivity based on the volume fraction of the solid (furnace feed) and of fluid (gaseous product), the radiation heat transfer between the particles in the packed bed is not accounted for. To improve the standard porous media model for high temperatures, the two conductivity values for the solid (furnace feed) and the fluid (gaseous product) were replaced by a single, temperature-dependent effective thermal conductivity value representing the actual furnace conditions more accurately. The aim was to have an effective thermal conductivity value to be used in Fluent that incorporated both conductive as well as particle–particle radiative aspects based on actual process conditions. To estimate such a thermal conductivity value, a representative volume of the packed bed was constructed by Adema.[22] This volume was created to be a stripped-down version of the real packed bed that contained pellets (and their corresponding density, thermal conductivity, and specific heat capacity values as determined by the industrial plant) as the solid phase and gas consisting of CO and P₄ as the gas phase. This volume is shown in Figure 7 and has dimensions of 17.1 × 5.7 × 5.7 cm.

Fig. 7

A view of the solved Fluent packed bed model with which the particle–particle radiation and effective thermal conductivity model was constructed.[22] The unit for the scale on the left is temperature

The stacking of the spheres (pellets) inside the volume has a cubic formation, and each of the spheres is 2 cm in diameter. The spheres have an overlap of 0.1 cm to simulate compaction and partial melting of the bed but mainly to improve meshing. The packed-bed structure has a porosity of 37.7 pct. A T-grid mesh with 470,000 cells was constructed. On the packed-bed volume, a hot and a cold face with set temperatures were defined. Temperature differences between the two faces were between 288 K (15 °C) and 333 K (60 °C). The remaining boundaries were defined as symmetry faces. No actual gas flow was modeled, and the standard porous model was disabled during modeling. The radiation was modeled using the DTRM, which assumes grey radiation and diffuse surfaces and does not include scattering. As the pore spaces are small, the absorption coefficient of the gas was assumed negligible, and the solid surface emissivity was assumed to be equal to 1. The effective thermal conductivity is determined using Eq. [15], which is based on Fourier’s law of conduction, as follows:

$$ k_{\text{eff}} = q\left[ {{\frac{L}{{A\left( {T_{\text{hot}} - T_{\text{cold}} } \right)}}}} \right] $$

(15)

Fluent is used to determine the heat flow q through the cold and hot faces. The k _eff was determined at different temperatures and with different temperature gradients. The result is a temperature-dependent value for k _eff ranging from 1.15 W/m × K at 373 K (100 °C) to 9.21 W/m × K at 2273 K (2000 °C) represented by Eq. [16] as follows:

$$ k_{\text{eff}} = 1e^{ - 6} \times T^{2} + 0.002 \times T + 0.41 $$

(16)

Equation [16] was used to generate the effective thermal conductivity values that incorporate both conductive as well as particle–particle radiative effects. Similar results for grids of 300,000, 500,000, and 1,000,000 cells and volume dimensions of 5.7 × 5.7 × 5.7 cm, 17.1 × 5.7 × 5.7 cm, and 28.5 × 5.7 × 5.7 cm confirmed grid and volume independence. Although a particle–particle radiation model, therefore, is not directly available in the full-scale model, the effect of such radiation is taken into account by means of the aforementioned effective thermal radiation value.

Results and Model Validation

In this chapter, the modeling results for the base case scenario are presented along with the validation methodology. It concludes by performing a parametric study by using the validated model. It reveals important new findings in regard to the process.

Base Case Model Results

Reduction Zone and P₂O₅ Consumption

Because the bulk of the reduction of phosphorus takes place in the solid–gas region (packed bed) of the furnace,[9,11] it is important to understand the characteristics and behavior of the main reduction zone within the feed where temperatures range from about 1423 K (1150 °C) to about 1773 K (1500 °C). The former temperature refers to the onset of P₂O₅ liberation, whereas the latter refers to the temperature in which the calcium silicate formation increases to a point where the pellet looses its structure and melts, thereby slowing down the diffusion process. The reduction zone as calculated by the model is indicated in Figure 8. It is within the boundaries of this thin zone where the main reduction occurs and the gaseous product is generated. Figure 8 also shows the decrease in the concentration of P₂O_5(pellets) from 29.1 wt pct (and scalar value of 230) to 2.8 wt pct P₂O_5(slag).

Fig. 8

A 2D depiction of the P₂O₅ concentration inside the furnace as it is consumed. Clearly shown is a narrow, gas–solid reduction zone where the bulk of the phosphorus is produced. The unit on the left is kg of P₂O₅ per cubic meter

This finding proves the existence of a narrow, gas–solid reduction zone where the bulk of the phosphorus is produced. This is an important insight. The upper boundary of this reduction zone mainly is determined by the chemical composition of the pellet, bed conductivity, and bed permeability (i.e., the major variables controlling conduction and convection phenomena within the bed). The lower boundary of this reduction zone (and, therefore, the thickness thereof) mainly is determined by the structural qualities of the sintered pellet, the porosity of the pellet, and to a lesser extent, the reaction rate. The liberated P₂O_5(g) must diffuse through the pellets pores and react with the carbon on the coke surface. This is the rate-limiting step, not the reaction rate. If within the reduction zone, the pellets were to collapse mechanically, then this diffusion step and the important convection phenomena transporting the energy would be impeded, and the lower boundary would be reached.

Temperature and Pressure Distributions

Figure 9 shows the temperatures inside the furnace along with the location of the reduction zone. Although it is not shown here because of the 1800 K (1527 °C) upper visual display limit of the graphics, the average temperature in the inlet zone is ±2700 K (2427 °C). The colder areas indicated with the two dotted ellipses are where cold, fresh feed at a temperature of 573 K (300 °C) is introduced artificially to the model (Eq. [2])—a situation that reflects reality. The location of the temperature probe in the outlet duct also is shown.

Fig. 9

A 2D depiction of the temperatures inside the furnace along with the reduction zone. The average temperature in the inlet zone is 2700 K (2427 °C), but because of the 1800 K (1527 °C) upper visual display limit of the graphics, it is not shown. The unit for the scale on the left is Kelvin

Figure 10 shows a 45-Pa pressure drop inside the packed bed as calculated by the porous media model (the standard condition of 1 atm was assumed).

Fig. 10

A 2D depiction of the pressure drop inside the packed bed. The unit for the scale on the left is relative pressure in Pascal

Gas Formation and Flow Characterization

Figure 11 shows the velocity contours of the gaseous product from the outlet ducts. The furnace under investigation produced ≈3.26 kg/s of gas in the base case scenario, and the model flow rate of 3.15 kg/s (1.58 kg/s per pipe) corresponds to within 3.4 pct of the actual gaseous flow rate. No constraints were put on this value when solving the model. It is essential for validation purposes and is generated entirely as a dependent variable.

Fig. 11

A 2D depiction of the velocity contours of the gaseous product for both outlets. The unit for the scale on the left is m/s

Figure 12 shows the direction of the velocity contours and shows the average calculated gas velocity to be 4.2 m/s. This finding provided insight into the reason for the recurring build up of gaseous outlet pressures; low, recirculating gaseous flow velocity areas in the top parts of the ducts caused dust accumulation that resulted in the increased pressures.

Fig. 12

A 2D depiction of the velocity vectors of the gaseous product for one of the outlets. The unit for the scale on the left is m/s

Energy Distribution Within the Furnace

According to Robiette and Allen[23] and Bailey,[24] the energy in an industrial phosphorus furnace is distributed between heating and melting of the material (≈40 pct) and chemical reactions (≈45 pct), whereas cooling losses (cooling water), electrical losses (Joule heating), and radiative heat losses account for the remaining energy (≈15 pct). For the base case scenario, the total energy input was 45.4 MW, and after convergence, the total energy output was 46.1 MW. The difference is less than 1.6 pct. The energy was distributed as follows:

• (a) 41 pct (18.5 MW) was used for chemical reactions

• (b) 40 pct (18.2 MW) was used for heating and melting of the material in the packed bed

• (c) 19 pct (8.7 MW) was attributed to cooling losses that include the latent heat of the gaseous product leaving the furnace

These values correspond well with the values in the literature. Figure 13 shows the location in the packed bed where the 18.2 MW is extracted as a result of heating and melting. No heating and melting energy are extracted in the reduction zone (Figure 8), only chemical reaction energy. The reduction zone is indicated with two broken arrows. Above the reduction zone, most energy is being used for heating the Moeller feed while some slag formation takes place. Below the reduction zone, most of the energy is used for melting as the pellets collapse to form slag together with the silica gravel. However, some energy still is used to heat up the coke.

Fig. 13

A 2D depiction of the energy requirements as a result of heat and the melting of the burden. The unit for the scale on the left is Watt

Base Case Model Validation

Several strategic temperature measurements at the production plant were used to validate the model.

The Radial Furnace Lining Temperature

Eight probes were installed radially inside the carbon lining of the furnace at a height corresponding to the reduction zone. The average temperature value of one of these probes measured during a period of two months was 585 K (312 °C) (±61 K (±61 °C)). This value agrees with the model-predicted temperature of 660 K (387 °C) in the corresponding position in the model within 2 standard deviations.[13]

The Vertical Furnace Lining Temperature B

Five probes were installed vertically in the furnace. The average temperature differences between the top and bottom, vertically installed, temperature probes calculated during the same two months was 70 K (±35 K) (±35 °C). This value agrees with the model temperature difference of 30 K (30 °C) within 2 standard deviations.[13]

The Gaseous Outlet Temperatures

Two online, gaseous outlet temperature values from the furnace were monitored continuously. The actual position of these probes in the gaseous outlet pipe of the furnace is shown Figure 9. The predicted temperature of 832 K (559 °C) from the model differs somewhat from the 728 K (455 °C) actually measured, but it still falls within the limits of the ±126 K (±126 °C) standard deviation value of the actual temperature measurement.

The Slag Temperature

An average slag temperature at the production plant was measured as 1750 K (1477 °C) (±50 K (±50 °C)) by a pyrometer (considering that the slag temperature at the plant still was measured only once a day in 2003, the representativeness of the value can vary). The temperature of the model is 1813 K (1540 °C) and agrees within 2 standard deviations.[13]

The P₂O_5(slag) Analysis from the Laboratory

The actual P₂O_5(slag) for the base case scenario was 1.6 wt pct compared with the 2.8 wt pct of the model prediction. Given (1) the small amounts of P₂O₅ in the slag, (2) the 16 pct standard deviation of P₂O_5(slag) in the actual samples,[13] and the fact that successive interpolations can lead to small errors in the surface integration report,[14] this prediction is good.

The Specific Power Consumption (SPC)

The most important parameter at the production plant is the SPC. It provides a no-nonsense assessment of the efficiency of the process. For the base case scenario, the actual SPC is 14.3 MWh/ton P₄, and the model SPC is 15.3 MWh/ton P₄, which amounts to an absolute difference between the two values of 6.9 pct. No constraints were put on the variables required for the calculation of the SPC while solving the model. These solved model SPC values form part of the overall validation process (Table VIII).

Table VIII Validation Data Comparing Actual and Model Values

Changing Production Rates

The remaining four scenarios identified in Table I were modeled. Four subsequent model output parameters along with the corresponding actual production values are shown in Figures 14 and 15. These parameters are as follows:

Fig. 14

The model vs the actual values of the P₂O_5(slag) and SPC values for scenarios 1 through 5 as listed in Table I

Fig. 15

The model vs the actual values of the gaseous product flow rates and the gaseous outlet temperatures for scenarios 1 through 5 as listed in Table I

• (a) SPC

• (b) P₂O_5(slag)

• (c) Gaseous product outlet temperature

• (d) Gaseous product flow rate

Figure 14 shows that small differences exist for the SPC values between the actual value and the model value. However, the relative differences remain fairly constant at high flow rates. A reason for the model SPC values being consistently higher than the actual values becomes clear when combined with the P₂O_5(slag) results. With the way SPC is defined, a decreased P₄ flow rate in the model will decrease the SPC.

With the additional P₂O₅ lost to the slag (2.8 pct − 1.6 pct = 1.2 pct for scenario 3), a corresponding amount of P₄ is not being produced in the model. This issue is reflected in the model’s slight underprediction of the gaseous product flow rate in Figure 15.

Figure 15 shows a diverging trend between the actual and model gaseous product temperatures. Here, the model shows better predictive accuracy at higher flow rates (i.e., when the furnace is running optimally). The model is steady state and therefore will mimic better a process closer to steady state. Low actual flow rates are the result of instabilities characterized by, for example, incidents of bridge formation and Moeller burden collapses. These instabilities can cause the average gaseous temperature during an 8-hour period to drop. For scenario 1, the 154 K (154 °C) difference (854 K − 700 K (581 °C − 427 °C)) shown in Eq. [17] represents an energy discrepancy of 0.39 MW − 2.1 pct of the energy used for heating and melting as follows:

$$ 0.39\,{\text{MW}} = 2.59{\frac{\text{kg}}{\text{s}}} \times 1000\,{\frac{\text{J}}{{{\text{kg}} \times {\text{K}}}}} \times {\frac{{154\,{\text{K}}}}{1,000,000}} $$

(17)

These predicted values are acceptable when considering the ±126 K (126 °C) standard deviation associated with the actual temperature measurements.

The model has been tested across a wide range of actual process data. The next step is to test the sensitivity of key parameters.

Parametric Study

Scheepers[13] performed a series of sensitivity analyses on nine important parameters by solving 37 scenarios. These parameters were furnace flow rate (residence time), fixed carbon-to-P₂O₅ ratio in the Moeller bed (C-fix), effective thermal conductivity, reaction kinetics, packed-bed porosity, P₂O₅ pct in the furnace feed, energy transported away from the furnace through the electrode cooling water, and electrode operating depths. In this article, the sensitivity of the following parameters are discussed: (1) reaction kinetics and (2) energy transported away from the furnace through the electrode cooling water

Five model output parameters along with the corresponding actual production values are shown in Figures 16 to 21. These parameters are as follows: (1) the average gaseous temperature at the outlet, (2) the average slag temperature, (3) the average packed bed temperature, (4) SPC, and (5) P₂O_5(slag).

Fig. 16

The influence of varying ore reaction rates on the slag-, packed bed-, and gaseous temperature values, in which scenario 6 has the highest reaction rate and scenario 8 has the lowest reaction rate

Influence of Reaction Kinetics

The reaction rate constants of type X-ore is shown in Equation 9. The following additional scenarios were created and solved:

• (a) Scenario 6; k _rr1 = 0.0003 T − 0.39 (the highest)

• (b) Scenario 7; k _rr3 = 0.0003 T − 0.40

• (c) Scenario 8; k _rr4 = 0.0001 T − 0.125 (the lowest)

In Figure 16, the gaseous product, the slag, and the average packed-bed temperatures for all four scenarios indicate little variation. The reason for this finding is the short time required for reduction[8,11,25] when compared with the 8 to 12 hours residence time of the Moeller in the furnace.

Even with the reaction rate increased (scenario 6) or decreased (scenario 8), the time available for reduction is more than enough considering the long time it takes for the burden to reach the slag zone. It also should be remembered that the starting temperature for Eq. [6] of 1423 K (1150 °C) (the temperature in which the equation’s equilibrium constant ≥ 1) is mostly independent of reaction kinetics, which means that the point in which the reactions for the scenarios commence always will start at the same heights within the packed bed. However, the thickness of the reduction zone will change.

Therefore, if the reaction rate is slowed excessively—as is the case with scenario 8 (Figure 17)—the bottom end of the reduction zone touches the slag interface more, and the P₂O_5(slag) starts increasing. However, at high reaction rates, melt formation on a microscale can block the porous channels of a pellet, thus preventing the P₂O₅ gas molecule to diffuse to the carbon and thus decrease the reduction rate. This situation is not considered in the model.

Fig. 17

The influence of varying ore reaction rate values on P₂O_5(slag) and SPC values. It is only at very low reaction rates, as in scenario 8, that the increased thickness of the reduction zone causes the bottom end of this zone to touch the slag interface

Influence of Electrode Water Cooling (Boundary Condition 3)

Plant operation revealed that during normal operational conditions, an average of 5 MW is lost through the three electrodes by means of the cooling water. This is in part caused by a large portion of the energy released by the arc being dissipated in the electrode.[25] In the model, the energy is introduced to the user-developed inlet (arc) zones, and no current actually flows through the electrodes. For this reason, the thermal conductivity value of the electrodes provided by the manufacturer in Table III was increased from 23 to 4500 W/m × K, using augmented thermal conduction to compensate the Joule heating of the electrodes. In this way, the required 5 MW could be extracted from the furnace model and thus prevent the subcooling of the top part of the electrode grid. It is recognized that such artificial augmentation could have an impact on modeling results. To assess its influence, the negative energy sink of −5 MW associated with scenario 1 was decreased to (1) scenario 9 (negative energy sink = −4.5 MW) and (2) scenario 10 (negative energy sink = −4 MW)

In Figure 18, the gaseous product, the slag, the average packed-bed temperatures, and the exposed electrode surface temperatures just above the packed bed are depicted.

Fig. 18

The influence of changing cooling flux to the top part of the electrode on the slag-, packed bed-, and gaseous predicted temperature values

With the electrode cooling at −4 MW, the exposed electrode surface temperature increased by ≈100 K (100 °C); at −4 MW electrode cooling flux, less energy is being withdrawn from the electrode while the same amount of energy is introduced to the inlet (arc) zones. This process increases the temperatures of the exposed electrode surface as well as the average temperatures of the gaseous product, the slag, and the average packed bed.

Results from Figure 19 show that the increased energy that remains in the furnace decreases the P₂O₅ loss in the slag, which decreases the gaseous product deviation and subsequently lowers the SPC.

Fig. 19

The influence of changing the cooling flux to the top part of the electrode on P₂O_5(slag) and SPC values

After the completion of the project, it was revealed that the actual value of the electrode cooling flux was not −5 MW but rather the significantly lower value of −1 MW. For this reason, the electrode cooling flux was adapted and a revised scenario 3 was solved again as scenario 11, which is the same as scenario 3 but with a negative energy sink = −1 MW

In Figures 20 and 21, the results of the revised scenario 30 were superimposed onto Figures 14 and 15.

Fig. 20

The model vs the actual values of the SPC and P₂O_5(slag) for scenarios 1 through 5 and scenario 11 with a 1-MW electrode cooling rate

Fig. 21

The model vs the actual values of the gaseous product flow rates and the gaseous outlet temperatures for scenarios 1 through 5 and scenario 11 with a 1-MW electrode cooling rate

Scenario 11 provided a P₂O₅ value of 2.1 pct (2.8 pct − 2.1 pct = 0.7 pct) and, therefore, has resulted in the gaseous product flow rate increasing to 3.23 kg/s—very close to the actual value of 3.26 kg/s. This increase, in turn, also led to the decrease in SPC to 14.8 MWh/ton P₄—also a lot closer to the actual value of 14.3 MWh/ton P₄. The reason for the improved results is the extra 4 MW introduced to the furnace. This energy allowed more apatite to react, thus producing the additional gaseous product as well as consuming more energy inside the packed bed. However, the gaseous outlet temperature did increase by 100 K (100 °C), which resulted in an additional 0.32 MW of energy leaving the furnace. Closer predictions are expected through optimization of the enthalpy (Figure 5) and heat capacity values (Figure 6) by using the expanded Factsage[14] databases currently available.

Conclusions

The article presents a model that successfully incorporates accurate, multifield thermodynamic, kinetic, and industrial data with computational fluid dynamic calculations. It is 3D and uses industrial measurements, laboratory experiments, and a combination of empirical and thermodynamic data to obtain accurate boundary conditions, initial values, and material specifications. It accounts for fully developed gas flows of gaseous products generated from within the packed bed; the energy associated with chemical reactions, heating, and melting; as well as thermal conductivity and the particle–particle radiation within the burden. The following conclusions were made:

• 1. The model proves the existence of a narrow, gas–solid reduction zone where the bulk of the phosphorus is produced. This is important an insight

• 2. The upper boundary of this reduction zone mainly is determined by the chemical composition of the pellet, bed conductivity, and bed permeability (i.e., the major variables controlling conduction and convection phenomena within the bed).

• 3. The lower boundary of this reduction zone (and therefore the thickness thereof) mainly is determined by the structural qualities of the sintered pellet, the porosity of the pellet, and to a lesser extent, the reaction rate. The liberated P₂O_5(g) must diffuse through the pellet pores and react with the carbon on the coke surface. This is the rate-limiting step, not the reaction rate. If within the reduction zone, the pellets were to collapse mechanically, then this diffusion step and the important convection phenomena transporting the energy would be impeded, and the lower boundary would be reached.

• 4. Fast reaction rates in a narrow reaction band combined with long residence times diminish the influence changing reaction rates have on the process. Only when the reaction rate is slowed excessively (artificially or naturally, e.g., secondary feed material) does the lower boundary of the reduction zone start touching the slag interface and increase P₂O_5(slag).

• 5. Most heat exchanged between new pellets entering the furnace and the gaseous product gas produced in the reduction zone takes place in the top 0.5 m of the furnace bed.

• 6. The gaseous product and velocity information show low, recirculating gaseous flow velocity areas that cause dust accumulation and thus result in an increased pressure measurement in the outlet duct.

• 7. For the Base Case scenario (closer to the normal production), the model provided values closer to reality in three of the four major output variables. This finding concludes that the model is robust, nonstiff (in a differential equation sense), and reflects reality (within the fundamental limitation).

The true power of this study is the extensive and methodical validation that ensures industrially endorsed results.