Dynamic Modeling and Simulation of Basic Oxygen Furnace (BOF) Operation

: Basic oxygen furnaces (BOFs) are widely used to produce steel from hot metal. The process typically has limited automation which leads to sub-optimal operation. Economically optimal operation can be potentially achieved by using a dynamic optimization framework to provide operators the best combination of input trajectories. In this paper, a ﬁrst-principles based dynamic model for the BOF that can be used within the dynamic optimization routine is described. The model extends a previous work by incorporating a model for slag formation and energy balances. In this new version of the mathematical model, the submodel for the decarburization in the emulsion zone is also modiﬁed to account for recent ﬁndings, and an algebraic equation for the calculation of the calcium oxide saturation in slag is developed. The dynamic model is then used to simulate the operation of two distinct furnaces. It was found that the prediction accuracy of the developed model is signiﬁcantly superior to its predecessor and the number of process variables that it is able to predict is also higher.


Introduction
The Basic Oxygen Furance (BOF) is responsible for approximately 70% of the steel production worldwide [1]. A schematic diagram for the BOF in shown in Figure 1. Scrap metal and hot metal are charged to the BOF, and an oxygen jet at supersonic speed is injected from the top through the lance onto the surface of the metal bath. Some of the species in the metal bath are oxidized and form a less dense slag layer. Flux is added to prevent refractory wearing and to contribute to slag formation. Millions of metal droplets are generated at the impact zone due to the impingement of the supersonic oxygen jet on the liquid metal. Carbon in the metal droplets can react with iron oxide in the slag forming carbon monoxide [2]. Experiments using X-ray fluoroscopy have shown that gas bubbles can be formed within the metal droplet itself [2,3] and not only at the metal droplet-slag interface. Carbon monoxide and carbon dioxide are also formed due to decarburization taking place at the impact zone. The droplet-gas-slag mixture, commonly referred to as emulsion, can occupy most of the furnace volume during the main blow. The process has limited automation and is highly dependent on operators' knowledge and past experience. Without a framework that can consistently aid the operators with the decision-making process, some of the more complex interactions between the process variables may be overlooked, leading ultimately to sub-optimal operation. Therefore, an optimization framework for the BOF operation that can consistently provide operators with the economically optimal operating conditions could greatly help steelshops reduce production costs. However, in order to have an optimization framework it is necessary to have an accurate enough dynamic model of the physical process.
Several dynamic models [4][5][6][7][8][9] for the BOF have been developed in recent years. Jalkanen [4] developed the CONSIM 5 simulator for the BOF. All reactions are assumed to take place in one zone and oxygen is partitioned according to its affinity to a certain element. The model accounts for slag formation, energy balance, scrap melting and decarburization.
Sarkar et al. [8] also developed a dynamic model for the BOF. They assumed that the refining reactions take place only in the emulsion zone, between elements dissolved in the droplets and FeO in slag. Oxygen is distributed between the elements in the upper part of the metal bath according to their concentration. The model does not include an energy balance and temperature is assumed to increase linearly during the blow.
Another dynamic model for the basic oxygen furnace was developed by Lytvynyuk et al. [7]. All the supplied oxygen is used to form iron oxide and all other refining reactions take place via reaction with FeO. The rate of the reactions taking place in the metal bath and slag is primarily dependent on the kinetics of mass transfer. For the energy balance, the slag and metal bath are assumed to have the same temperature. The model seems to give excellent prediction of the slag composition, metal bath composition and metal bath temperature at the end of the blow.
A few dynamic models based primarily on empirical relationships have also been developed. Kattenbelt and Roffel [5] used step responses to model the decarburization rate during the main blow. They studied the process response to step changes in lance height, oxygen flow rate and iron ore addition. The model makes extensive use of empirical relationships but provides limited physical insight.
A very comprehensive, first-priniciple dynamic model for the kinetics of decarburization in the BOF was developed by Dogan et al. [6]. Their model accounts for carbon oxidation in two zones: The emulsion and the impact zone (Figure 1). At the impact zone, carbon in the metal bath is assumed to react directly with oxygen and carbon dioxide, and in the emulsion zone carbon in the metal droplets is oxidized by FeO in the slag. The metal bath and slag temperature are assumed to increase linearly with blowing time and the slag composition is required as an input. Rout et al. [9] continued the work of Dogan et al. [6] by modeling the kinetics of manganese, phosphorus and silicon oxidation enabling the model to predict the slag composition.
The current study extends the work of Dogan et al. [6] by incoporating a model for slag formation as well as energy balances. The model for decarburization in the emulsion zone is modified to account for recent fidings [10] and the model for scrap melting is updated based on recent studies [11]. The mathematical model for the BOF is implemented as a system of Differential Algebraic Equations (DAEs) using CasADi [12] with a Python front-end. The resulting framework is used to simulate the Cicutti data [13] as well as 71 heats for the BOF of an industrial collaborator (Plant A).

Mass and Energy Balances
The phenomena taking place in the BOF are quite complex, therefore several assumptions are made for the derivation of the mass and energy balances. Flux, iron ore, metal droplets and gas bubbles are considered to be uniformly dispersed in the slag phase. The metal droplets, gas bubbles and slag form an emulsion for which the continuous phase is assumed to be the slag. The oxidation reactions take place mainly at the impact zone (IZ), the interface between the oxygen jet and the metal bath, and there is no resistance to the diffusion of oxides from the metal bath to the slag. A schematic representation of the flow of material in the BOF is shown in Figure 2. The reactions modelled by the current work are as follows [14], where heats of reaction, ∆H rxn , were obtained from FactSage at 1900 K: Decarburization: Desiliconization: Iron Oxidation: Post-combustion: Argon or nitrogen can be injected at the bottom of the BOF to improve mixing. Even though the injected nitrogen can undergo reaction at BOF operating conditions, the amount of nitrogen compounds formed is small and can be neglected. Therefore, it is assumed that the stirring gas leaves the furnace unreacted. At the impact zone, iron, carbon and silicon react with the injected oxygen forming iron oxide, carbon monoxide and silica as shown in Equations (1)- (5). Some of the CO may further react with O 2 forming CO 2 (Equation (6)).
The mass of dissolved flux and oxides is promptly incorporated by the slag. Similarly, the mass of melted scrap is assimilated by the metal bath. All Fe in the iron ore is assumed to be in the form of magnetite (Fe 3 O 4 ) that is reduced to FeO by carbon in the metal bath (Equation (7)).
∆H rxn = 184.660 kJ/mol (7) The impact of the oxygen jet on the liquid metal causes millions of metal droplets to be ejected in the emulsion zone. Carbon in the metal droplets can reduce FeO in the slag according to the reaction shown in Equation (3). The refined droplets then fall back to the metal bath.
The gases CO, CO 2 , N 2 /Ar and unreacted oxygen O 2 form the off-gas stream. Owing to the high temperatures, the residence time for the gases is assumed to be negligible in this work. A mass balance for the metal bath gives: where W b represents the mass of the metal bath.Ẇ sc is the scrap melting rate,Ẇ C,IO is the rate at which carbon in the metal bath is consumed to reduce magnetite in the iron ore to FeO,Ẇ D,b−e andẆ D,e−b are the mass flow rate of droplets from the metal bath to the emulsion and from the emulsion to the metal bath, respectively.Ẇ C,b−e ,Ẇ Fe,b−e ,Ẇ Si,b−e are the mass flow rate of carbon, iron and silicon out of the metal bath and equals the oxidation rate of the respective element at the impact zone. A similar mass balance on the slag-metal-gas emulsion gives: where W sme is the mass of the slag and droplets in the emulsion.Ẇ DissFlux is the rate of lime and dolomite dissolution,Ẇ FeO,IO denotes the mass flow rate of iron oxide coming from iron ore,Ẇ FeO,b−e , W SiO 2 ,b−e are the mass flow rate of iron oxide and silicon dioxide into the slag and are equal to the rate of formation of the respective component at the impact zone.Ẇ CO,EZ is the rate carbon monoxide is formed at the emulsion zone due to droplet decarburization. Heat is generated in the BOF by the oxidation and post-combustion reactions (Equations (1) and (4)-(6)) and is consumed by scrap melting, to heat the fluxes and iron ore, by decarburization in the emulsion (Equation (3)) and carbon dioxide reduction (Equation (2)). The total heat Q rxn gen generated by the oxidation reactions is given by where Q rxn FeO is the heat from iron oxidation (Equation (5)), Q rxn SiO 2 is the heat from silicon oxidation (Equation (4)), Q rxn CO is the net heat generation from decarburization at the impact zone (Equations (1) and (2)), Q rxn CO 2 is the heat of CO post-combustion. Q rxn i is obtained by multiplying the rate of formation of compound i at the impact zone by the respective heat of reaction.
Individual energy balances were derived for the emulsion and metal bath. To simplify the problem, the following assumptions are made: • The gas hold up in the slag and metal bath is negligible • The gas exits the furnace at the slag temperature • The temperature of the slag-metal-emulsion and metal bath is uniform • A fraction α p l of all the heat generated is assumed to be lost to the environment and the remainder is assumed to be absorbed by the slag.
An energy balance for the slag-metal-gas emulsion yields: with In the above, W s is the mass of slag, W D is the mass of the droplets in emulsion, C P,s , C P,b , C P,i are the heat capacity of the slag, molten metal and component i, respectively. Q Fe,EZ is the heat consumed by the decarburization reaction in the emulsion, Q D,b−e is the heat required to raise the temperature of the incoming droplets to the temperature of the slag-metal-gas emulsion; Q FeO,b−e , Q SiO 2 ,b−e , Q CO,b−e , Q CO,b−e and Q CO 2 ,b−e are the heat required to raise the temperature of the iron, silicon, carbon oxide and carbon dioxide formed at the impact zone to the emulsion temperature, respectively; Q e−b is the heat lost from the emulsion to the metal bath, Q f lux and Q IO are the heat consumed by the flux and iron ore additions, Q N 2 /Ar,b−e is the heat required to raise the temperature of the stirring gas from the bath to the emulsion temperature, Q O 2 ,b−e is the heat required to raise the temperature of the non-reacted oxygen from the bath to the emulsion temperature. The rate of heat transfer between the slag and metal bath is given by: where R is the furnace diameter and h s−b is the heat transfer coefficient, which can be estimated using dynamic data. Applying an energy balance for the metal bath gives: where Q D,e−b is the heat transferred by the droplets coming from the emulsion to the metal bath, and Q sc is heat used to melt the scrap, Q O 2 ,in and Q N 2 /Ar,in are the heat consumed to raise the temperature of oxygen and stirring gas to the bath temperature: where F i , T i,0 are the inlet mass flow rate and temperature of gas i.

The Impact Zone
The impact zone in the BOF is the region where the oxygen jet impinges on the molten metal bath ( Figure 1) causing a deformation on the liquid surface that is assumed to have the shape of a paraboloid of height h c and maximum diameter d c . The number of impact zones is equivalent to the number of nozzles n n in the lance as long as there is no coalescence of the oxygen jet. The area of one impact zone is equal to the surface area of the paraboloid: where d c and h c are calculated using the correlations proposed in [15]. For the current work, the decarburization rates previously used by Dogan et al. [6] are modified by a parameter α p in order to account for the difference between the conditions at which the equation rates were derived and the BOF operating conditions. These parameters can be tuned using process data. Therefore, the decarburization rate via O 2 (Equation (1)) is -2r O 2 ,iz , where: with: The partial pressure in Equation (17) is in atm, and the parameter α p Si,C in Equation (18) accounts for the inhibiting effect that silicon has on the rate of carbon oxidation [14]. Similarly, the decarburization rate via CO 2 (Equation (2)) is given by: In the above expressions, k a and k O 2 are the rate coefficients calculated using the correlations found in Dogan et al. [6] and P is the partial pressure calculated as: where P a is the ambient pressure. In Equation (20), F O 2 and F N 2 /Ar are the inlet molar flow rate of oxygen and bottom stirring gas, and F CO and F CO 2 are the molar flux of carbon monoxide and carbon dioxide formed from the decarburization reaction in Equations (1) and (2) and the post-combustion reaction in Equation (6). When the carbon content of the metal bath [%C] becomes lower than the critical carbon content C C , the decarburization rate is given by [6]: where V b is the volume of liquid metal and k m is the rate constant calculated using the correlation developed by Kitamura et al. [16]. The desiliconization rate is calculated using the expression found in Rout et al. [17] and modified by a parameter α p Si to give: where ρ b is the density and [%Si] is the silicon content of the liquid metal. Rout et al. [9] found that the equilibrium silicon content of the metal bath ([%Si eq ]) is approximately zero, therefore it is neglected in the calculations. For the current study, it is assumed that the rate of iron oxidation is proportional to the partial pressure of oxygen and is given by: All the oxygen injected in the system via the lance that is not used for decarburization, desiliconization or iron oxidation, is assumed to be used in the post-combustion of CO.

Scrap Melting
In the BOF, the heat generated by the oxidation reactions (Equations (1) and (4)-(6)) is much higher than that required to reach the targeted end-point temperature. Therefore, scrap metal is usually added in order to absorb part of the surplus heat via melting. In this section the model used to compute the scrap melting rateẆ sc as well as the heat absorbed Q sc is presented.
Consider a scrap metal plate of half-thickness L, initial temperature T sc0 and melting temperature T m , as shown in Figure 3a. This plate is then submerged in a metal bath at temperature T b > T m and carbon content [%C]. Assuming that heat transfer occurs only in the axial direction and constant physical properties for the plate, a heat balance at the interface between the solid plate and the metal bath gives [11,[18][19][20]: where k sc is the thermal conductivity of the scrap plate, h sc is the heat transfer coefficient, ρ sc is the scrap density, ∆H sc is the latent heat of melting of the scrap, C P,sc is the specific heat of the scrap, A sc is the interfacial area, x is the distance of a point within the scrap from the interface, t is time and T sc (x, t) is the scrap temperature at a position x. A schematic representation can be found in Figure 3b. Equation (24) can be solved using the Quasi-Static approach to yield Equations (25)-(27) [21]: where α = k sc /ρ sc C P sc is the thermal diffusivity. The scrap melting rateẆ sc is then given by: where n sc is the total number of scrap plates. At the beginning of the BOF operation, the temperature of the molten metal may not be high enough to melt the scrap, in which case the interface temperature T m is equal to the bath temperature T b . The interface temperature is otherwise assumed to correspond to the melting temperature. This is shown in Equation (29) : The melting temperature T m is dependent on the carbon content C m at the melting interface and can be calculated using Equation (30) [22]: Dogan [22] assumed C m to be equal to the carbon content in the bulk metal [%C]. This assumption reflects, to some extent, what happens in the BOF: Carbon in the metal bath migrates to the scrap surface lowering its melting temperature [21,23]. One consequence of such an assumption is that scrap types of similar physical properties will all melt at the same time, independent of their carbon content. On the other hand, if C m is assumed to be equal the carbon content of the scrap C sc , low-carbon content scrap types will only melt towards the very end of the blow, which is not observed in practical BOF operations. For the current work, C m is defined as: The heat transfer coefficient h sc in Equation (25) is calculated according to the correlation proposed by Gaye et al. [24] and given by: where˙ is the mixing power due to bottom stirring and top blowing in Wm −3 , and h sc is in units of Wm −2 K −1 . The rate Q sc at which heat is absorbed by the scrap is given by the heat conduction term in Equation (24) before melting starts, and by the convective term once melting starts as shown in Equation (33).

Iron Ore Dissolution
Some steel shops also add iron ore before or during the blow to cool down the metal bath and meet the desired end point temperature. Since the mass of iron ore added is usually small compared to the amount of scrap, a rather simple model is used to account for the cooling effect of iron ore.
Iron ore composition can vary greatly according to the region but the predominant element is usually Fe followed by oxygen. For the current model all Fe in the iron ore is assumed to be in the form of magnetite (Fe 3 O 4 ) that is reduced to iron oxide by carbon in the metal bath according to Equation (7). Furthermore, the iron ore melting point T m is considered to be equal to the FeO melting point (1644 K). Reduction of magnetite and melting of the formed iron oxide are assumed to happen concomitantly.
The iron ore particles are assumed to be spherical and to have uniform temperature, therefore the total heat transfered to the surface of an iron ore particle can be calculated using Equation (34): where T ore is the iron ore temperature, h ore is the heat transfer coefficient, r ore is the radius and T s is the slag temperature. It is assumed that a fraction T ore /T m of the total heat Q conv,ore contributes to melting of the iron ore and the remainder is used for sensible heating. A similar approach was used by MacRosty and Swartz [25] and Bekker et al. [26] to model scrap melting in electric arc furnaces. At the beginning of the process T ore << T m and most of the heat is used to heat up the iron ore. As T ore approaches T m , the fraction of the total heat used for melting increases. An energy balance for a single iron ore particle yields: where W ore is the mass of a single iron ore particle and can be obtained from Equation (36): In the above, ρ ore is the iron ore density and C P,ore is the iron ore heat capacity. A heat balance at the interface of the iron ore particle gives Equation (37) for the rate of dissolution of an iron ore particle: where ∆H ore is the latent heat of melting of FeO, C P,FeO is the heat capacity of iron oxide and A ore is the area of the interface, here equal to the surface area of a sphere of radius r ore . The heat consumed to reduce magnetite to iron oxide is given by: where y Fe,ore is the mass fraction of Fe in the iron ore, ∆H rxn is the heat for the reaction defined in Equation (7), M Fe is the molar mass of iron andẆ ore,melt is the melting rate of iron ore given by: The total heat consumed by iron ore melting and reduction is given by: where n ore is the number of iron ore lumps.

Flux Dissolution
Similarly to the decarburization rates at the impact zone, flux dissolution is modelled after Dogan et al. [27]. The flux particles are assumed to have a spherical shape, and the rate of flux dissolution is proportional to the rate of change of the particles' radius.
The rate of lime dissolution is given by Equation (41) [27,28], where r L is the particle radius, (%CaO) is the concentration of CaO in slag, (%CaO sat ) is the saturation concentration of CaO in the slag, k L is the mass transfer coefficient, ρ s and ρ L are the slag and lime density, respectively.
The parameter α p L is introduced here to account for deviations between the experimental conditions at which k L was derived, and the BOF operating conditions.
For dolomite, the rate of dissolution is given by [27,29]: where M MgO and M CaO are the molar mass of MgO and CaO, respectively, (%MgO) is the concentration of MgO in slag, (%MgO sat ) is the saturation concentration of CaO in the slag and k D is the mass transfer coefficient. The rate constants k L and k D are calculated using the correlations proposed by Dogan et al. [27] with the parameter used to modify the Reynolds number β set to the nominal value of 1. Data for the calcium oxide saturation in slag (%CaO sat ) for different slag compositions were obtained using the Cell Model [30], a Matlab program that gives the saturation concentration of the individual species in the slag as a function of composition and temperature. Kadrolkar et al. [30] validated the Cell Model results against FactSage TM and ThermoCalc TM . The function 'fitnlm' in Matlab was then used to fit a curve to the generated data, yielding Equation (43).
An R 2 of 0.9 was obtained for the nonlinear regression. The goodness of Equation (43) to predict the calcium oxide saturation in slag was evaluated using the Cicutti data [13]. The calculated values for CaO sat using the Cell Model and Equation (43) are shown in Figure 4, where it can be seen that there is excellent agreement between them. The rate at which heat is absorbed by the flux particles is given by: where T s is the temperature of the slag, T i is the flux temperature, h i is the heat transfer coefficient determined using data and n i is the number of flux particles. Assuming a uniform temperature for the flux particle, an energy balance gives Equation (45) for the temperature of an iron ore particle: where C P,i is the heat capacity and W i = 4/3ρ i πr 3 is the mass of one single flux particle.

Decarburization in the Emulsion Zone
In Figure 5, the initial carbon composition and diameter of a metal droplet are represented by C i,0 and D 0 , respectively. The metal droplet is formed from the interaction of the oxygen jet and the metal bath at the impact zone at time t 0 , and ejected to the slag-metal-gas emulsion. The droplet stays in the emulsion zone for R t seconds, then it returns to the metal bath at time t 0 + R t having a composition C i, f and diameter D f . While in the emulsion zone, carbon in the metal droplet can reduce FeO in slag according to reaction 3. The droplet generation rate R B , given by Equation (46), is calculated using the empirical expression proposed by Subagyo et al. [31] modified by a parameter p DG to account for possible differences between the conditions at which the correlation was derived and a given BOF operation: whereV O 2 is the volumetric flow rate of oxygen and N B is the blowing number [31].
Since millions of droplets are generated at every point, it can become computationally expensive to track the composition of the individual metal droplets. Therefore, an algebraic equation to calculate the final carbon content of the metal droplets was developed and the following assumptions were made in order to make the problem tractable: • All droplets decarburize immediately after they are ejected from the impact zone • The carbon content of the metal droplets in the emulsion zone is approximated as the average carbon content of the individual metal droplets • Except for carbon, the mass of other elements in the metal droplets stays constant whilst the droplet travels through the emulsion zone The first-principles model developed by Kadrolkar and Dogan [10] for the droplet decarburization was used to generate data for the final carbon content (C C, f ) of individual droplets with respect to its initial carbon content C C,0 , the slag temperature T s and composition. Data were generated for the range of initial carbon content between 0.3% and 5% and slag temperature of 1623-2153 K, and for the slag composition provided by Cicutti et al. [13]. Equations (47) and (48) give a good description of the generated data with an R 2 of 0.96 and 0.92, respectively.
In Figure 6, the average final carbon content of the droplets reported by Cicutti et al. [13], the values predicted by Kadrolkar and Dogan [10]'s model and what was obtained using Equation (47) are shown. Similar to the first-principles model proposed by Kadrolkar and Dogan [10], the final carbon content of the droplets predicted by Equation (47) is largely within the range of values reported by Cicutti [13] for a real BOF operation. These results indicate that Equation (47) approximates the kinetics of droplet decarburization reasonably well and can therefore be used within the BOF model to describe the kinetics of decarburization in the emulsion zone.
with α p µ,D given by: α whereẆ C,s is the rate of carbon removal, α p µ is a parameter and µ s is the slag viscosity. The initial carbon content of the droplet C C,0 is equal to the carbon content of the metal bath C b at the time of ejection. The term α p µ,D is necessary because the effect of slag viscosity is not taken into account in Equations (47) and (48). However, it is likely that the high slag viscosity at the beginning of the blow significantly decreases the decarburization rate. At high viscosities the slag becomes less fluid, negatively impacting the rate of FeO mass transfer to the droplet surface and potentially leading to FeO depletion in the neighborhood of the droplet, decreasing the decarburization rate. The value of α p µ can be determined using dynamic data.
The mass of carbon W C,s in the emulsion is then given by: where the first term on the right hand side is the amount of carbon entering the emulsion zone minus the carbon removed via decarburization, and the second term is the flow rate of carbon returning to the metal bath. C C,s is the average carbon content of the droplets in the emulsion calculated by: and W D,s is the total mass of droplets in the emulsion determined as: The droplet residence time R t is calculated using the following empirical expression: As long as the residence time R t of the droplets in the emulsion is small, C C,s is approximately equal to C C, f and the simplifying assumption of uniform carbon content for the droplets in the emulsion does not significantly impact the final result, while avoiding the need to track the behavior of individual droplets.

Implementation
The mathematical model presented in the previous sections was implemented as a DAE system using CasADi [12] with Python 3.7. An index-1 DAE in semi-explicit form can be written in the form given by Equations (55) and (56):ẋ where f and h are the differential and algebraic functions, respectively, x and z are the differential and algebraic states, u and p represent the control variables and time-independent parameters, and t is time. Within CasADi, the DaeBuilder class was selected since it is capable of performing model reduction by eliminating algebraic variables that can be explicitly calculated. The resulting DAE system was solved using the variable step size DAE solver IDAS [32]. For implementation of the mathematical model, the following modifications were made where necessary: • Equations of the form: where a → 0, were rewritten as: where is a positive small number. This strategy was used on the submodels for flux dissolution, iron ore and scrap melting to prevent division by zero since either the radius or thickness of the particles continuously decreases with time and can eventually equal zero. • Piecewise functions of the form: were rewritten using hyperbolic tangent functions: where γ is an adjustable parameter that controls the steepness of the continuous switching function approximation. This was used for flux dissolution (Equation (42)), scrap melting (Equation (29)), decarburization in the emulsion zone (Equations (47) and (48)), among others for a smooth transition and to ensure differentiability. • Flux additions: Flux and iron ore can be added at anytime during a blow, and each individual addition is modeled as shown in Section 2.5. To model the indiviudal flux additons a new variable t ij , where i is the flux type (lime, dolomite, iron ore) and j is the addition number (first, second, third), is defined for the flux addition time. Given the radius r ij of the flux added at time t ij , Equation (41) for lime dissolution rate can be reformulated as: which was implemented using a hyperbolic tangent function.

System Parameters and Input Data
The parameters α p introduced during the model development to account for distinct BOF operations and differences between the conditions at which the respective equations were derived and the operating conditions were manually adjusted for a data set available in the literature for a 200 ton furnace [13], as well as for the data provided by Plant A for 70 heats for a 250 ton furnace. The final values of the parameters α p are given in Table 1. Table 1. The values of model parameters for simulation and optimization. The heat transfer coefficients h are given in Wm −2 K −1 , and α p Fe is given in kgm −2 Pa −1 s −1 . All the other parameters are dimensionless.

Parameters
Cicutti In Cicutti et al. [13]'s study, lime is added before the blow starts and at every minute up to 7 min, whereas dolomite is added before the blow starts and again at 7 min. Ar/N 2 gas is continuously injected from the bottom of the furnace to aid with stirring. In Plant A operations, lime and dolomite are added only before the blow starts, scrap selection varies for each heat, and the furnace design does not allow for bottom stirring.

Simulation Results for Cicutti's Operations
Cicutti et al. [13]'s data have been used to validate several dynamic models developed for the BOF [6,8,9,33]. Information regarding the hot metal and scrap compositions is shown in Table 2. For the mass of flux, gravel and iron ore refer to Cicutti et al. [13]. Table 2. Mass and composition of hot metal and scrap types added to the BOF, and average scrap thickness [22].

Input Hot Metal Heavy Scrap Light Scrap Pig Iron External Scrap
Mass ( Figure 7 shows the metal bath temperature and carbon content predicted by the model described in this article, as well as the data published by Cicutti et al. [13].  [13] and predicted values for the temperature of liquid metal and slag. (b) Comparison between measured [13] and predicted values for the carbon content of liquid metal and the returning metal droplets.
In the first few minutes, the bath temperature decreases due to the heat absorbed by the scrap. Thereafter, the temperature increases approximately linearly as heat is released by the oxidation reactions. The oscillations in the slag temperature during the first half of the blow are due to the flux additions, done at every minute. Figure 7b shows that the carbon content predicted by the model for the bulk metal agrees well with the measured data. The predicted final carbon content of the returning droplets is also shown in Figure 7b. Due to the high interfacial area, the rate of carbon refining at the droplet level is significantly higher than that for the bulk metal, thus the lower carbon content.
The decarburization rate in the emulsion zone is shown in Figure 8a. The extent of droplet decarburization is primarily dependent on the initial carbon content of the liquid metal droplet when it is ejected from the impact zone [10]. At low initial carbon contents the liquid metal droplets do not bloat and their residence time in the emulsion zone decreases significantly [10,34]. Due to that the contribution of the emulsion zone to the total decarburization rate decreases significantly towards the end of the blow. It follows that at high carbon contents it is possible to increase the decarburization rate in the emulsion zone by increasing the droplet generation rate.
It is possible to identify the three decarturization periods characteristic of the BOF operation on the total decarburization rate graph shown in Figure 8a: • Period I: As the silicon content of the metal bath decreases, the decarburization rate increases • Period II: Decarburization rate stays approximately constant • Period III: Decarburization rate is controlled by mass transfer of carbon in the metal bath The change from Period II to III is shown by the dashed line in Figure 8a, but it can also be seen in the compostion of the gases exiting the furnace in Figure 8b. As the decarburization at the impact zone decreases, more oxygen becomes available for the post-combustion reaction explaining the increase in the percentage of CO 2 at approximately 14 min. The percentage of CO 2 in Figure 8b is slightly higher than would normally be observed in practice. This can be due to the ideal assumption that all the oxygen not used in the oxidation reactions is consumed in the post-combustion of CO. The profile for the lance height and oxygen flow rate is shown in Figure 9a. The effect of lance height changes on the decarburization rate is clear at 4 min and 7 min in Figure 8a: Lowering the lance height increases the droplet generation rate, as well as the rate constants for the decarburization reactions taking place at the impact zone, which leads to a higher decarburization rate. The contribution of the decarburization in the emulsion to the total decarburization is significantly lower for the present work than previously suggested [8,17]. Rout et al. [17] suggested that 76% of the total decarburization happens in the emulsion, while in the modeling approach adopted by Sarkar et al. [8] decarburization only takes place in the emulsion zone. For the current paper, the emulsion zone was responsible for 15% of the total carbon removed. The reason why the contribution of the decarburization in the emulsion is lower for the present study is because of the significantly lower droplet generation rate R B . Sarkar et al. [8] modified the droplet generation rate (Equation (46)) by a factor of 15. Using a modified correlation for R B , Rout et al. [35] obtained a droplet generation rate similar in magnitude to Sarkar et al. [8]. It is not currently viable to measure how much decarburization occurs at the impact and emulsion zones individually, and it may be the case that a different set of parameter values yields approximately the same total decarburization rate. However, taking into account the gradual decrease in the decarburization rate in the emulsion zone in Figure 8a, it can be inferred that for a very high contribution of the emulsion to the total decarburizaton rate, Period II would no longer be characterized by an approximately constant decarburization rate.
The slag composition throughout the blow is presented in Figure 10 for SiO 2 , CaO, FeO and MgO. There is a good agreement between the values predicted by the model and the data. The large content of silicon dioxide in the slag during the first minute is due to 800 kg of gravel addition. A large SiO 2 content increases the slag viscosity, reducing the decarburization rate at the emulsion and allowing FeO to build up. Flux dissolution slows down significantly at high slag viscosities, but as the blow proceeds SiO 2 gets diluted by FeO and the flux dissolution rate increases. The error between the model prediction and measured data is, most likely, due to the treatment of slag as a homogeneous phase for the density and viscosity calculations [36]. A comparison between the carbon content prediction by the present and previous [6,8,9] studies is shown in Figure 11. Only for the current model, the carbon content prediction starts from time zero, and enregy balances are included; moreover, the quality of the prediction itself is quite good compared with previous works. This is also the only study for which the mathematical model was transposed as a DAE system and integrated using a variable step size solver, whilst in the aforementioned works integration was carried out using a fixed step size. The first main benefit stemming from the current implementation is the reduced computational time as shown in Table 3. Moreover, the convergence is taken care of by the integrator, and convergence studies based on step size are not required. Secondly, the dynamic model can be easily built within an optimization framework to determine the optimal input trajectories. Figure 11. Comparison of the carbon content prediction by different models [6,8,9] and the measured values [13] for a 200-ton furnace.

Simulation Results for Plant A Operations
The developed framework model was also used to simulate 71 heats using the data provided by Plant A for its BOF operation. The model parameters, shown in Table 1, were manually adjusted based on data collected at the end of the blow for the end-point carbon content, slag composition and temperature. The extreme conditions in the BOF often prevent samples from being collected during the blow, therefore there is no data regarding the state of the system during the blow.
The goodness of the prediction was evaluated in terms of the average of the model predictions and process data for the 71 heats. Figure 12 shows the quality of the prediction for the end-point slag composition together with the standard deviation. The results indicate that the model is able to predict the slag composition at the end of the blow reasonably well. The carbon content of the metal bath can be measured using a carbon probe before tapping, or a sample is collected and the measurement is done in a laboratory setting, with the second giving the more accurate results. However, laboratory results are not available for all the heats for the current study. The average end-point carbon content and temperature measured using the probe and the model's prediction for the 71 batches are shown in Figure 13, based on which it is possible to conclude that the model performs fairly well. The model was able to predict the end-point carbon content of 80% of the heats with a precision of ±0.03%, and the end-point temperature of 61% of the heats with a precision of ±30 K. The prediction error can be attributed to several factors. If the furnace is used back-to-back, there is always some slag left over from one heat to another, however its temperature and mass are not measured and the initial mass of slag present in the system is not accurately known. Moreover, the initial temperature of the hot metal charged to the furnace is also not precisely known since it is measured, on average, 40 minutes before the blow starts. Furthermore, due to the lack of data, it is assumed that all scrap types have similar properties, which can lead to inaccurate prediction of the heat needed for complete melting. Another potential reason could be the assumption that the off-gas and slag have the same temperature; however, the data available is not sufficient to estimate the actual off-gas temperature during the process. This, with the absence of sub-models to account for phosphorus oxidation, manganese oxidation and the formation of higher iron oxides, can explain the poorer performance of the model in predicting the end-point temperature.
For each heat, the mass of iron ore, flux, scrap and hot metal charged to the furnace is different, explaining the variation in the predicted and measured values. Since the inputs are always changing, it makes it challenging to identify the primary source for larger deviations between the simulated and actual data. Plant A has three different input profiles for the oxygen flow rate and lance height, however the overall quality of the predictions is not significantly affected by the different control profiles indicating that the model can potentially be used to simulate a wide range of operating conditions.
The trajectories for the temperature, carbon content and slag composition are shown in Figure 14 for one of the heats. We observe that the trajectories are similar to the Cicutti case. The content of SiO 2 is lower for Plant A operation at the beginning of the blow because no gravel is added. Moreover, the rate of silicon oxidation is slower due to the lack of bottom stirring. As mentioned before, at low SiO 2 contents, decarburization in the emulsion is enhanced due to the reduced slag viscosity. This is evident in Figure 14a by the larger difference between the carbon content of the droplets and the carbon content of the metal bath at high FeO contents at the beginning of the blow. As expected, due to the lack of bottom stirring, the FeO content of the slag increases much faster than that for Cicutti's case toward the end of the blow. The rapid increase can also be explained by the higher lance height towards the end of the blow compared to the mid-blow. The scaled control profiles for the oxygen flow rate and lance height for the given heat are shown in Figure 9b.  Figure 14 indicates that the model is able to predict the end-point carbon content, silicon content and temperature of liquid metal reasonably well. It also provides insights about the trajectory followed by the variables and how the change in one or more of the inputs affects the system. The developed framework could potentially aid in finding more profitable operation modes as well as in understanding the more complex relationships between process variables.

Conclusions
The dynamic mathematical model presented in this study extends the work of Dogan et al. [6] by incorporating slag formation and energy balances. The models for the decarburization in the emulsion zone and scrap melting were updated with more recent findings, and an empirical relationship to calculate the calcium oxide saturation was obtained. Therefore, the main phenomena taking place in the BOF were accounted for, making the model suitable for further studies.
The mathematical model was translated as a DAE system to an open source environment (Python with CasADi [12]). The current implementation allowed for a significantly shorter simulation time compared with previous studies [6,8,9]. Moreover, as more complex and detailed models for the phenomena taking place in the BOF are developed, they can be incorporated relatively easily using the current framework.
It was shown that the updated model gives a better prediction of the carbon content trajectory for Cicutti et al. [13]'s data compared with the previous version, and more recently developed dynamic models for the BOF [8,9]. The dynamic model was also used to simulate 71 heats for a real industrial BOF operation. The model predictions for the end-point carbon content, slag composition and temperature agree reasonably well with the data for a wide range of operating conditions.