*2.7. 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):

$$\dot{\mathbf{x}}(t) = f(\mathbf{x}(t), \mathbf{z}(t), \mathbf{u}(t), \mathbf{p}) \tag{55}$$

$$0 = h(\mathbf{x}(t), \mathbf{z}(t), \boldsymbol{\mu}(t), \boldsymbol{p}) \tag{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:

$$y(a) = \frac{1}{a} \tag{57}$$

where *a* → 0, were rewritten as:

$$y(a) = \frac{1}{a + \epsilon} \tag{58}$$

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:

$$y(\mathfrak{c}) = \begin{cases} y\_1(\mathfrak{c}) & a > b \\ y\_2(\mathfrak{c}) & a \le b \end{cases} \tag{59}$$

were rewritten using hyperbolic tangent functions:

$$\bar{y}(c) = y\_1(0.5\tanh(\gamma(a-b)) + 0.5) + y\_2(0.5\tanh(\gamma(b-a)) + 0.5) \tag{60}$$

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 *tij*, 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 *rij* of the flux added at time *tij*, Equation (41) for lime dissolution rate can be reformulated as:

$$\frac{dr\_{ij}}{dt} = \begin{cases} 0 & t < t\_{ij} \\ k\_L \frac{\rho\_s}{100p\_L} (\% \text{CaO}\_s - \% \text{CaO}\_{\text{sat}}) & t \ge t\_{ij} \end{cases} \tag{61}$$

which was implemented using a hyperbolic tangent function.

#### **3. Results and Discussion**

#### *3.1. System Parameters and Input Data*

The parameters *α<sup>p</sup>* 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 *α<sup>p</sup>* are given in Table 1.



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/N2 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.

#### *3.2. 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].


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].

**Figure 7.** (**a**) Comparison between measured [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:


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 CO2 at approximately 14 min. The percentage of CO2 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.

**Figure 8.** (**a**) Total decarburization rate and decarburization rate at the emulsion zone and (**b**) composition of the off-gas stream exiting the BOF.

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.

**Figure 9.** (**a**) Control profile for Cicutti et al. [13]'s data and (**b**) scaled control profiles for a heat from Plant A.

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 *RB*. Sarkar et al. [8] modified the droplet generation rate (Equation (46)) by a factor of 15. Using a modified correlation for *RB*, 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 SiO2, 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 SiO2 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 SiO2 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].

**Figure 10.** Evolution of slag composition for the Cicutti data [13] and model prediction: (**a**) FeO, (**b**) SiO2, (**c**) CaO, (**d**) MgO.

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.


**Table 3.** Simulation time for Cicutti et al. [13]'s data required in different studies.
