Law of Mass Action Based Kinetic Approach for the Modelling of Parallel Mass Transfer Limited Reactions: Application to Metallurgical Systems

Abstract

The objective of this paper was to present a new law of mass action based rate expression for mass transfer limited reversible reactions. A simple reaction model was derived for parallel oxidation of silicon, chromium and carbon under conditions relevant to the argon-oxygen decarburization (AOD) process. Our hypothesis is that when the forward rate coefficients approach infinite values, the composition at the reaction surface approaches a constrained equilibrium. In numerical analysis, however, only finite numbers are allowed and therefore only finite values are accepted for rate coefficients. In order to circumvent this problem, additional residual affinity constraints were introduced. This assures that the affinities of all the reactions at the reaction front reach a pre-defined non-zero residual affinity and the rate coefficients remain finite. The calculated equilibrium composition is essentially the same as that obtained with the equilibrium coefficient method. In the case of effective gas side mass transfer, the component having the highest mole fraction or the highest mass transfer rate on the liquid side consumes most of the oxygen. When the gas side mass transfer rate is decreased, the mass transfer rate of oxygen begins to limit the overall rate and the partial pressure of oxygen at the reaction interface decreases. Then, the role of interfacial equilibrium becomes important as the species start to compete for the oxygen. The proposed method provides a transparent and direct solution of the mass transfer limited reaction rates and is thus suitable for process simulators and CFD software.

1. Introduction

In the primary and secondary metallurgy of steelmaking, temperature of the metal bath is typically in the range of 1500–1800°C, but local temperatures can be much higher. In basic oxygen furnace, for example, the temperature of the gas jet impact area may exceed 2000°C.^1,2,3) Comparable circumstances are found also in other industrial applications: the temperature in the modern combustion furnaces varies from 1088–1173 K (815–900°C) in fluidized boilers, 1923 K (1650°C) in pulverized coal boilers and up to 2200–2603 K (1927–2330°C) in premixed gaseous burners.⁴⁾ In many of the aforementioned applications, the local kinetic reaction rate and the rate expression at the reaction interface is unknown and virtually impossible to measure. The resistance of the interfacial reaction decreases as the temperature increases. The temperature dependency of the reaction rate coefficient k can be described with the Arrhenius equation, Eq. (1), or the Eyring equation, Eq. (2).⁵⁾   

$$k=Aexp\left(-\frac{E_{\text{a}}}{RT}\right)$$(1)

  

$$k=\frac{k_{\text{B}}T}{h}exp\left(\frac{\mathrm{\Delta }{S}^{*}}{R}\right)exp\left(-\frac{\mathrm{\Delta }{H}^{*}}{RT}\right)$$(2)

As seen from Eqs. (1) and (2), the reaction rate coefficient becomes large at high temperatures, and thus the mass transfer or the mixing becomes the rate limiting step. It is thus reasonable to assume that the reaction front can reach the chemical equilibrium composition that is subject to mass transfer constraints. However, there are examples of reactions having slow kinetic rates at high temperatures. In particular, the presence of surface active elements can hinder the interfacial reaction to the extent that the interfacial reaction rate becomes the rate limiting step. A well-known example is the de-nitrification reaction of liquid steel.⁹⁾

The boundary condition for the rate expression of a single mass transfer limited reaction can be determined directly with the equilibrium constant method. One of the first qualitative descriptions of the boundary layer theory was proposed by Nernst¹⁰⁾ in 1904. In the case of parallel reactions, however, the description of mass transfer in the boundary layer needs to be coupled with a description of the competitive thermodynamic equilibrium at the interface.^7,11) The equilibrium constant method is not applicable as such, particularly if the rate limiting step is not known a priori or it varies during the process.

The mass transfer resistances of the considered phases depend on the observed geometry. Typical geometries in the primary and secondary steelmaking processes include gas jet in contact with surface of a steel bath (Fig. 1(a)), a liquid metal droplet entrained in gas flow (Fig. 1(b)) as well as slag droplet in liquid metal or liquid metal droplet in slag (Fig. 1(c)). Small dots on the reaction surface present existing small slag particles. The mass transfer rates of the reagents and reaction products are interrelated through the reaction stoichiometry.

Fig. 1(a).

Gas jet in contact with surface of a steel bath.

Fig. 1(b).

Steel droplet entrained in gas flow.

Fig. 1(c).

Gas bubble immersed in liquid steel.

In this work, a novel approach based on the Law of Mass Action was studied and developed further in order to provide a thermodynamically consistent treatment of parallel mass transfer controlled reactions in the context of steelmaking. The method can be used in metallurgical and chemical processes involving mass transfer controlled reactions. The method discussed in this paper has already been employed by the authors in different cases for defining the mass transfer limited rates of reactions. The studied applications in the AOD process include mathematical modelling of reactions in a single gas bubble in liquid steel, during side-blowing decarburization and during the reduction of slag. The method has also been applied for mathematical modelling of chemical heating in the CAS-OB process. The objective of this paper is to re-derive, discuss and further develop the law of mass action based method for multiple parallel reversible reactions. A generalized approach for implementation into practice is also provided and discussed.

2. Law of Mass Action Based Approach

The Law of Mass Action (LMA), proposed by Waage and Guldberg^13,14) in 1864, is a mathematical expression of solutions in dynamic equilibrium. More specifically, it relates the thermodynamic equilibrium to the reaction kinetics and implies that forward and backward reaction rates must be equal at equilibrium. The reaction rate of an arbitrary reaction can be expressed by:   

$$R^{″}=\underset{\text{Forward}\mathrm{\hspace{0.17em} }\text{reaction}}{\underbrace{k_{\text{f}}{\displaystyle \prod }_r{x}_r^{{\nu }_r}}}-\underset{\text{Backward}\mathrm{\hspace{0.17em} }\text{reaction}}{\underbrace{k_{\text{b}}{\displaystyle \prod }_p{x}_p^{{\nu }_p}}}=k_{\text{f}}\left({\displaystyle \prod }_r{x}_r^{{\nu }_r}-\frac{{{\displaystyle \prod }}_p{x}_p^{{\nu }_p}}{K}\right)$$(3)

where k_f is the forward reaction rate coefficient, k_b is the backward reaction rate coefficient x is the molar fraction, r denotes reactants and p denotes reaction products. In order to be generally applicable, it is necessary to replace the concentrations with activities, which describe effective concentration:   

$$R^{″}=k_{\text{f}}\left({\displaystyle \prod }_r{a}_r^{{\nu }_r}-\frac{{{\displaystyle \prod }}_p{a}_p^{{\nu }_p}}{K}\right)$$(4)

Here, the expression shown in Eq. (4) is referred to as the modified Law of Mass Action. In the following, a case of n parallel reactions between oxygen and dissolved elements in the liquid steel is studied. The term key component refers to the species that has been used to define the rate expression and has molar stoichiometric coefficient of –1. The reaction equations and equilibrium constants can be written either in terms of the oxidized elements i, Eq. (5), or in terms of the O₂ as the key component, Eq. (6).   

$$i+{\nu }_i{\text{O}}_2\rightleftharpoons i{\mathrm{O}}_{2{\nu }_i}\hspace{1em}\text{for}\hspace{1em}i=1,\mathrm{\hspace{0.17em} }2\dots n$$(5)

  

$$\frac{1}{{\nu }_i}i+{\text{O}}_2\rightleftharpoons \frac{1}{{\nu }_i}i{\mathrm{O}}_{2{\nu }_i}\hspace{1em}\text{for}\hspace{1em}i=1,\mathrm{\hspace{0.17em} }2\dots n$$(6)

where ν_i is the molar stoichiometric coefficient of O₂ in the oxidation reaction of species i and iO_2νi refers to the oxide of species i. The equilibrium constants corresponding to Eqs. (5) and (6) are obtained from Eqs. (7) and (8), respectively. In this work, an atmospheric pressure of 1 atm was assumed.   

$$K_i=exp\left(-\frac{\mathrm{\Delta }{G}_i^{\circ }}{RT}\right)=\frac{a_{i{\mathrm{O}}_{2{\nu }_i}}}{a_i{p}_{{\mathrm{O}}_2}^{{\nu }_i}}$$(7)

  

$$K_i^{\text{G}}=K_i^{1/{\nu }_i}=exp\left(-\frac{\mathrm{\Delta }{G}_i^{\circ }}{RT{\nu }_i}\right)=\frac{a_{i{\mathrm{O}}_{2{\nu }_i}}^{\frac{1}{{\nu }_i}}}{a_i^{\frac{1}{{\nu }_i}}{p}_{{\text{O}}_2}}$$(8)

The source terms for all other species that participate in the reaction are calculated simply by multiplying the rate of the key component by their stoichiometric coefficients. In this work, the oxidation of Fe, i.e. the solvent, is neglegted. FeO is often an intermediate product and acts effectively as oxygen storage. It is reasonable to expect that in an actual process, the supplied oxygen first oxidizes iron to form iron oxide, which is then reduced by species with higher oxygen affinity. In a transient case, the composition at the reaction surface approaches asymptotically the constrained chemical equilibrium. The form of the chemical rate expression is irrelevant as long as the equilibrium composition and the mass transfer control can be achieved. Therefore, it is convenient to employ the modified law of mass action, which satisfies the equilibrium condition as a limit, to formulate the reaction rate expressions.¹⁵⁾ The law of mass action based approach has been shown to hold for elementary reactions of radical gas species containing only a single reaction step and to yield the correct equilibrium composition.^15,16) The surface mole balance for the elements i = 1...n diffusing from the bulk liquid steel and reacting at the surface can be written as follows:   

$${\beta }_{\text{L}}\frac{c_{\text{L}}}{M_{\text{L}}}\left(x_i^{\infty }-x_i\right)={\beta }_{\text{L}}\frac{c_{\text{L}}}{M_{\text{L}}}\left(x_i^{\infty }-\frac{a_i}{{\gamma }_i}\right)=k_{\text{f}}\left(a_i{p}_{O_2}^{{\nu }_i}-\frac{a_{i{\mathrm{O}}_{2{\nu }_i}}}{K_i }\right)$$(9)

where x denotes the mole fraction and γ_i is the Raoultian activity coefficient of species i. The rate expression is written here using the law of mass action for Eq. (9). For the sake of simplicity it is assumed that the transport rate or the amount of oxide species is not rate-controlling; the activities of the oxides were treated as constants. Here, a value of 0.5 was employed. The mass balance for the oxygen diffusing from the gas phase to reaction surface and distributing between all the species can be defined as follows:   

$${\beta }_{\text{G}}\frac{p_{\text{G}}}{RT}\left(x_{{\text{O}}_2}^{\infty }-\frac{p_{{\text{O}}_2}}{{\gamma }_{{\text{O}}_2}}\right)={\displaystyle \sum }_{i=1}^r{k}_{\text{f}}\left(a_i^{\frac{1}{{\nu }_i}}{p}_{{\text{O}}_2}-\frac{a_{i\mathrm{O}}{}_{2{\nu }_i}^{\frac{1}{{\nu }_i}}}{K_i^G }\right)$$(10)

By using the mass balances shown in Eqs. (9) and (10), the activities of dissolved elements and oxygen can be solved from Eqs. (11) and (12), respectively.   

$$a_i=\frac{\frac{ {\beta }_{\text{L}}{x}_i^{\infty }}{k_{\text{f}}}\frac{c_{\text{L}}}{M_{\text{L}}}+\frac{a_{i{\mathrm{O}}_{2{\nu }_i}}}{K_i }}{\frac{ {\beta }_{\text{L}}}{k_{\text{f}}{\gamma }_i}\frac{c_{\text{L}}}{M_{\text{L}}}+p_{{\text{O}}_2}^{{\nu }_i}}\hspace{1em}\text{for}\hspace{1em}i=1,\mathrm{\hspace{0.17em} }2\dots n$$(11)

  

$$p_{O_2}=\frac{\frac{{\beta }_{\text{G}}{x}_{{\mathrm{O}}_2}^{\infty }}{k_{\text{f}}}\frac{p_{\text{G}}}{RT}+{{\displaystyle \sum }}_{i=1}^r\left(\frac{a_{i\mathrm{O}}{}_{2{\nu }_i}^{\frac{1}{{\nu }_i}}}{K_i^{\text{G}} }\right)}{\frac{{\beta }_{\text{G}}}{k_{\text{f}}{\gamma }_{{\text{O}}_2}}\frac{p_{\text{G}}}{RT}+{{\displaystyle \sum }}_{i=1}^r\left(a_i^{\frac{1}{{\nu }_i}}\right)}$$(12)

As seen from Eqs. (11) and (12), the system consists of n + 1 variables and n + 1 equations, and hence numerical solution is possible. The main hypothesis of the law of mass action based approach is that as the rate of the surface reaction is infinitely fast and the reaction front is able to reach equilibrium. Mathematically, this means that if k_f → ∞, the following limits are obtained for the activities of species i and O₂:   

$$\underset{k_{\text{f}}\to \infty }{lim}{a}_i=\frac{a_{i{\mathrm{O}}_{2{\nu }_i}}}{K_i{p}_{{\mathrm{O}}_2}^{{\nu }_i}}=a_{i,\text{eq}}\hspace{1em}\text{for}\hspace{1em}i=1,\mathrm{\hspace{0.17em} }2\dots n$$(13)

  

$$\underset{k_{\text{f}}\to \infty }{lim}{p}_{{\mathrm{O}}_2}=\frac{{{\displaystyle \sum }}_{i=1}^r\left(\frac{a_{i\mathrm{O}}{}_{2{\nu }_i}^{\frac{1}{{\nu }_i}}}{K_i^G}\right)}{{{\displaystyle \sum }}_{i=1}^r\left(a_i^{\frac{1}{{\nu }_i}}\right)}=p_{{\text{O}}_2,\mathrm{\hspace{0.17em} }\text{eq}}$$(14)

It should be noted that Eq. (13) yields the equilibrium activities for the parallel reactions, while Eq. (14) defines the equilibrium activity oxygen, which has to be same for all reactions. The final surface composition can only be obtained after solving the system of non-linear equations defined by Eqs. (9) and (10), which include the mass transfer rates. It is important to understand that the composition solved at the reaction interface is not the full chemical equilibrium that can be determined with traditional methods, such as the equilibrium constant method or Gibbs energy minimization, but rather a constrained chemical equilibrium: the equilibrium permitted by the mass transfer rates onto and from the reaction interface. In essence, the LMA method studied in this paper is a rate expression for mass transfer limited dynamic equilibrium.

Based on the analytical expressions for the activities of oxygen and oxidizing species at the reaction front it can be deduced that their values approach the equilibrium values when the forward rate coefficient k_f → ∞. This means that the results obtained from the law of mass action method approach the solution obtained with the traditional equilibrium constant method. In comparison to the equilibrium constant method, the main advantage the LMA method is that its starting point is the kinetic rate expression and the control volume approach. As discussed broadly by Rao,¹²⁾ the equilibrium constant method is numerically very stiff and the use of the standard Newton’s method might be numerically unstable.

2.1. Numerical Implementation

The implementation of the law of mass action based approach requires a finite value for the forward rate coefficient k_f. Therefore, it is necessary to develop a robust way to define the values of k_f such that the requirements for accuracy and numerical stability are fulfilled. An analytical expression is not available due to non-linearity, and numerical methods need to be used. It should be noted that the employed approach functions correctly also in the case that the forward rate coefficient is known a priori. Therefore, the approach is applicable also for reactions, which are not limited by mass transfer.

As an illustrative numerical example, a simple case of liquid stainless steel surface exposed to O₂–CO gas is studied; a schematic illustration is shown in Fig. 1(a). In the simplified setting, the injected oxygen may oxidize carbon, chromium and silicon according to Eqs. (15), (16) and (17), respectively. As pointed out by Wei and Zhu,¹⁷⁾ FeO is essentially an intermediate product and hence the oxidation of iron was not accounted for. The values of the changes in standard Gibbs free energy were taken from Rao.¹²⁾ More specifically, the reference states of Si, Cr and C were taken as pure liquid, pure liquid and graphite, respectively. The reference states of SiO₂ and Cr₂O₃ were taken as β-cristobalite and pure solid, respectively. Ideal gas law was employed for the gaseous species.   

$${\text{Si}}_{\left(\text{l}\right)}+{\text{O}}_{2\left(\text{g}\right)}\rightleftharpoons {\text{SiO}}_{2\left(\text{l}\right)}\hspace{1em}\mathrm{\Delta }{G}_0^{\circ }=-938\mathrm{\hspace{0.17em} }913+193.719\mathrm{\hspace{0.17em} }T$$(15)

  

$$\begin{array}{l}{\text{Cr}}_{\left(\text{l}\right)}+\frac{3}{4}{\text{O}}_{2\left(\text{g}\right)}\rightleftharpoons \frac{1}{2}{\text{Cr}}_2{\text{O}}_{3\left(\beta -\text{cristobalite}\right)}\\ \mathrm{\Delta }{G}_1^{\circ }=-566\mathrm{\hspace{0.17em} }934+128.323\mathrm{\hspace{0.17em} }T\end{array}$$(16)

  

$${\text{C}}_{\left(\text{graphite}\right)}+\frac{1}{2}{\text{O}}_{2\left(\text{g}\right)}\rightleftharpoons {\text{CO}}_{\left(\text{g}\right)}\hspace{1em}\mathrm{\Delta }{G}_2^{\circ }=-119\mathrm{\hspace{0.17em} }025+83.482\mathrm{\hspace{0.17em} }T$$(17)

The corresponding equilibrium constants were defined according to   

$$K=\text{exp}\left(-\frac{\mathrm{\Delta }{G}^{\circ }}{RT}\right)=\frac{{{\displaystyle \prod }}_p{a}_p^{{\nu }_p}}{{{\displaystyle \prod }}_r{a}_r^{{\nu }_r}}$$(18)

For the sake of simplicity, it was assumed that all species behave ideally i.e. γ_i = 1 for all species i. Morever, the Raoultian activities of SiO₂ and Cr₂O₃ were assumed to be 0.5. In the case of transient applications also the bulk phase compositions need to be accounted for. However, the calculations presented in this paper are related to a single moment in time only, and hence the composition of the species in the bulk phases including the slag was assumed to be constant. Employing the stationary medium approach, the conservation equations can be written as follows:   

$${\beta }_{\text{L}}\frac{{\rho }_{\text{L}}}{M_{\text{L}}}\left(x_{\text{Si}}^{\infty }-a_{\text{Si}}\right)-k_{\text{f},0}\left(a_{\text{Si}}{p}_{{\text{O}}_2}-\frac{a_{{\text{SiO}}_2}}{K_0}\right)=0=f_0$$(19)

  

$${\beta }_{\text{L}}\frac{{\rho }_{\text{L}}}{M_{\text{L}}}\left(x_{\text{Cr}}^{\infty }-a_{\text{Cr}}\right)-k_{\text{f},1}\left(a_{\text{Cr}}{p}_{{\text{O}}_2}^{0.75}-\frac{a_{{\text{Cr}}_2{\text{O}}_3}^{0.5}}{K_1}\right)=0=f_1$$(20)

  

$${\beta }_{\text{L}}\frac{{\rho }_{\text{L}}}{M_{\text{L}}}\left(x_{\text{C}}^{\infty }-a_{\text{C}}\right)-k_{\text{f},2}\left(a_{\text{C}}{p}_{{\text{O}}_2}^{0.5}-\frac{1-p_{O_2}}{K_2}\right)=0=f_2$$(21)

  

$$\begin{array}{l}{\beta }_{\text{G}}\frac{p_{\text{G}}}{RT}\left(x_{{\text{O}}_2}^{\infty }-p_{{\text{O}}_2}\right)-{\beta }_{\text{L}}\frac{{\rho }_{\text{L}}}{M_{\text{L}}}\left(x_{\text{Si}}^{\infty }-a_{\text{Si}}\right)\\ -0.75{\beta }_{\text{L}}\frac{{\rho }_{\text{L}}}{M_{\text{L}}}\left(x_{\text{Cr}}^{\infty }-a_{\text{Cr}}\right)-0.5{\beta }_{\text{L}}\frac{{\rho }_{\text{L}}}{M_{\text{L}}}\left(x_{\text{C}}^{\infty }-a_{\text{C}}\right)=0=f_3\end{array}$$(22)

In addition, the system was set subject to the following constraints:   

$$a_{{\text{SiO}}_2}=x_{{\text{SiO}}_2}=0.5$$(23)

  

$$a_{{\text{Cr}}_2{\text{O}}_3}=1-x_{{\text{SiO}}_2}$$(24)

In addition to these it is necessary to find sufficiently large values for rate coefficients k_f,k that satisfy the equilibrium requirement. At the equilibrium, the affinity i.e. the driving force of each reaction should be equal to zero, resulting in definitions of equilibrium constants. As discussed earlier, an accurate equilibrium can be solved only if the forward rate coefficient approaches infinity. This is not computationally possible. Instead of having zero affinities at the reaction front, a concept of residual affinity is introduced and employed to move the reaction front thermodynamically to a small distance away from the equilibrium, making the solution numerically possible. In addition, the same value of residual affinity was employed for all reactions which means that the thermodynamic driving forces of all reactions are balanced at the reaction front. The absolute value of the residual affinity is very small in comparison to the standard Gibbs free energy of reaction; for the studied reactions the absolute values of ΔG° are in the order of 10⁵–10⁶ J/mol at 1873 K. The residual affinity of a reaction can be expressed as follows:   

$$-A=\mathrm{\Delta }{G}^{\circ }+RTln\left(\frac{{{\displaystyle \prod }}_p{a}_p^{{\nu }_p}}{{{\displaystyle \prod }}_r{a}_r^{{\nu }_r}}\right)$$(25)

The residual affinities of reactions 0, 1 and 2 need to be re-written for implementation as residuals:   

$$a_{\text{Si}}{p}_{{\text{O}}_2}-\frac{a_{{\text{SiO}}_2}}{exp\left(-\frac{A+\mathrm{\Delta }{G}_0^{\circ }}{RT}\right)}=0=f_4$$(26)

  

$$a_{\text{Cr}}{p}_{{\text{O}}_2}^{0.75}-\frac{a_{{\text{Cr}}_2{\text{O}}_3}^{0.5}}{exp\left(-\frac{A+\mathrm{\Delta }{G}_1^{\circ }}{RT}\right)}=0=f_5$$(27)

  

$$a_{\text{C}}{p}_{{\text{O}}_2}^{0.5}-\frac{1-p_{{\text{O}}_2}}{exp\left(-\frac{A+\mathrm{\Delta }{G}_2^{\circ }}{RT}\right)}=0=f_6$$(28)

With the exception of the residual affinities added to the change standard Gibbs free energy of reaction, Eqs. (26), (27) and (28) correspond to the traditional equilibrium constant method. The resulting relative error in the equilibrium constant can be calculated by comparing the restricted equilibrium coefficient $K_k^{*}$ to the true equilibrium coefficient, K_k:   

$$\text{Error}=1-\frac{K_k^{*}}{K_k}=1-\frac{\text{exp}\left(-\frac{A+\mathrm{\Delta }{G}_k^{\circ }}{RT}\right)}{\text{exp}\left(-\frac{\mathrm{\Delta }{G}_k^{\circ }}{RT}\right)}=1-\text{exp}\left(-\frac{A}{RT}\right)$$(29)

Assuming a temperature of T = 1873 K (1600°C), and assigning A with values of 100, 1, 0.01 and 0.0001 J/(mol·K) yields an relative error of 6.4×10,^–3 6.4×10,^–5 6.4×10^–7 and 6.4×10,^–9 respectively. Thus it can be concluded that the error in the equilibrium constant is negligible. The above 7×7 non-linear system can be solved effectively with numerical methods, for example the Newton’s method. Mathematically, this can be expressed in matrix form as follows:   

$$\mathbf{J}\times \mathrm{\Delta }\mathbf{y}=-\mathbf{f}$$(30)

where J is the Jacobian matrix, Δy is the column vector of activities and rate coefficients and f is the column vector of residuals. The Jacobian matrix consists of the first-order differentiates of the residuals with respect to free variables (see Table 1). The vectors Δy and f are defined as follows:   

$$\mathrm{\Delta }\mathbf{y}={\left[\mathrm{\Delta }{a}_{\text{Si}},\mathrm{\Delta }{a}_{\text{Cr}},\mathrm{\Delta }{a}_{\text{C}},\mathrm{\Delta }{p}_{{\text{O}}_2},\mathrm{\Delta }{k}_{\text{f},0},\mathrm{\Delta }{k}_{\text{f},1},\mathrm{\Delta }{k}_{\text{f},2}\right]}^{\text{T}}$$(31)

  

$$\mathbf{f}={\left[f_0,f_1,f_2,f_3,f_4,f_5,f_6\right]}^{\text{T}}$$(32)

Table 1. Elements of the Jacobian matrix, $J_{i,j}=\frac{\partial f_i}{\partial J_j}$ .

Corrections to variables y are obtained by solving Eq. (30) and updating them according to Eq. (33).   

$${\mathbf{y}}^{\mathit{n}+1}={\mathbf{y}}^{\mathit{n}}+\mathrm{\Delta }\mathbf{y}$$(33)

A small numerical constant, ε, was required for the numerical stability in the case of negative exponents in the Jacobian matrix. Here, a value ε = 10⁻²⁰⁰ was adopted. Initial guess value of activities was 0.5 for all species and the iteration was continued until the root-mean square (RMS) of the correction vector was less than 10⁻¹⁶. This required typically only 5–20 iterations. The value of A was set equal to 0.001 J/mol and thus the driving residual force was equal for all the studied reactions.

An additional and important aspect of the new method is the possibility to study reversible reactions. As discussed above, the interfacial composition is shifted a small distance away from the equilibrium; this distance is measured using the residual affinity A. The value of A is always set to be negative which means that in principle, the direction of the reaction is always from left to right for a positive value of k_f. Then, in the case of a reversed reaction, it would be necessary to change the sign of the residual affinity. But in practice, for this approach we would need to monitor the condition by a conditional if-statement, which would cause an undesirable discontinuity. Instead, we used the following simple approach to treat reversed reactions.

As an example, the oxidation reaction of Cr is studied in the following. As discussed earlier, the forward rate coefficient in the law of mass action rate expression should approach to ∞ in order to achieve equilibrium condition and fulfill the definition of equilibrium coefficient. To avoid the need to use infinite k values we defined a slightly falsified equilibrium, namely Eq. (27), having the residual affinity (Gibbs free energy) of finite value on the negative side, if the reaction proceeds from left to right. We now solve the “falsified” Cr₂O₃ activity from Eq. (27) and substitute it to the rate of law of mass for Cr oxidation, resulting in Eq. (34) that then defines the rate of oxidation.   

$${R_f^{″}}_{}=k_{\text{f},1}{a}_{\text{Cr}}{p}_{{\text{O}}_2}^{0.75}\left(1-exp\left(-\frac{A}{RT}\right)\right)$$(34)

Following the definition of the residual affinity, Eq. (25), the value of this rate is always positive with positive values of A and k_f,1. It should also be noted that the value in the parenthesis is equal to the relative error in the equilibrium constant, see Eq. (29). Physically, the numerical values of k must to be always positive. But as a trick, being numerically a free variable, we could let k have negative values too in the case of reversed reaction. Assuming the same values for activities, but an opposite sign of A, results in a reversed reaction. More specifically, the driving force is moved to the other side of the equilibrium.   

$${R_b^{″}}_{}=k_{\text{f},1}{a}_{\text{Cr}}{p}_{{\text{O}}_2}^{0.75}\left(1-exp\left(\frac{A}{RT}\right)\right)$$(35)

This rate is always negative with positive values of A and k_f,1. The condition on the surface can be such that the mass transfer wants to turn the reaction opposite but with the Eq. (34), this would not be possible and Eq. (35) should be used. Our question is now: if we set the sign of A to be positive and let the sign of the k_f change in Eq. (34), do we obtain the same result as by using Eq. (35)? The ratio of negative rate (Eq. (34)) and the accurate reversed reaction (Eq. (35)) can be written as follows:   

$$\begin{array}{l}\frac{-R_{\mathrm{f}}^{″}}{R_{\text{b}}^{″}}=\frac{exp\left(-\frac{A}{RT}\right)-1}{1-exp\left(\frac{A}{RT}\right)}=\frac{exp\left(-X\right)-1}{1-exp\left(X\right)}\\ \hspace{1em}\hspace{1em}\hspace{0.17em}=\frac{1-\text{X}+\frac{{\text{X}}^2}{2}-\frac{X^3}{6}\dots -1}{1-1-\text{X}-\frac{{\text{X}}^2}{2}-\frac{X^3}{6}\dots }=\frac{-\text{X}+\frac{{\text{X}}^2}{2}-\frac{X^3}{6}\dots }{-\text{X}-\frac{{\text{X}}^2}{2}-\frac{X^3}{6}\dots }\cong 1\end{array}$$(36)

When X = $\frac{A}{RT}$ approaches zero, the value of $\frac{-R_{\text{f}}^{″}}{R_{\text{b}}^{″}}$ approaches unity. Numerically, at 1600°C, substituting A = 1000 J/mol, 1 J/mol or 0.001 J/mol yields 0.937801, 0.999358 and 0.999994, respectively. Therefore, the use of negative values for k_f is justified when the value of A is sufficiently small.

3. Results and Discussion

By using the new LMA approach, we calculated the momentary mass transfer controlled rate of reactions between O₂, Si, C and Cr and composition at the steel surface. The studied cases consider a single moment of time and hence bulk phase composition was assumed to be constant.

3.1. Effect of Mass Transfer Coefficients

The effect of the ratio of gas and liquid side mass transfer on the rate limiting mechanism was studied at carbon mole fractions of 0.04 and 0.01. These correspond to the process conditions during the early stage of carbon removal in the AOD process. The model parameters were taken to be as β_L = 5×10⁻⁴ m/s, T = 1873 K (1600°C), x_Cr = 0.17, x_Si = 0.002 and A = 0.001 J/mol. The employed liquid phase mass transfer coefficient corresponds roughly to the highest experimental value presented by Chatterjee et al.,⁴⁸⁾ who measured apparent liquid phase mass transfer coefficients during top-blowing of oxygen on liquid silver. Figures 2 and 3 present the rate control of different mechanisms and the selectivity of oxygen, i.e. the share of oxygen consumed by each of the studied reactions. The values calculated by the method of Wei and Zhu¹⁷⁾ are also presented in Figs. 2 and 3, as defined later (and calculated) by Eqs. (41), (42), (43). These are not dependent on the mass transfer coefficients. This example illustrates the flexibility of the proposed method as well as its sensitivity to boundary conditions.

Fig. 2.

O₂ selectivity and rate limiting mechanisms, x_C=0.04.

Fig. 3.

O₂ selectivity and rate controlling mechanisms, x_C=0.01.

The simplified model was derived in such a way that it solves automatically the mass transfer controlled rates, the surface composition, and also the values for rate coefficients required to reach a pre-defined residual affinity A. The hypothesis was that if the values of rate coefficients are increased, the rates become mass transfer limited and the surface will reach equilibrium. Consequently, the next task is to determine the values of rate coefficients k_f so that they satisfy the requirement set for the residual affinity when the gas side mass transfer coefficient is changed. Figures 4 and 5 show these rate coefficients obtained as a solution of the system of equations.

Fig. 4.

Rate coefficients required for 0.001 J/mol residual affinities, x_C=0.04.

Fig. 5.

Rate coefficients required for 0.001 J/mol residual affinities, x_C=0.01.

Three different regimes can be observed: a) the rate is controlled by the mass transfer of oxygen, b) the rate is controlled by the mass transfer of the liquid phase and c) Cr₂O₃ is reduced and the rate is control by the mass transfer of the liquid phase. As seen in Figs. 2 and 3 when the gas side mass transfer is very effective typically β_G > 1.3 m/s, there is an excess of O₂ available at the surface. All species on the liquid side are then driven by the maximum diffusion rates towards the surface and selectivity of O₂ is defined by the molar composition of the bulk liquid and the stoichiometry. When the gas side mass transfer coefficient decreases below this critical value, the amount of O₂ at the surface decreases and the mass transfer of O₂ becomes the rate-determining step. In that case, the selectivity of O₂ is defined by the combination of the equilibrium and the liquid side mass transfer rates and a different reaction mechanism is obtained. As the gas mass transfer coefficient decreases the selectivity of the carbon and the silicon increases. If the mass transfer rate of O₂ is close to or smaller than the mass transfer rates of C and Si, Cr is no longer oxidized. As the oxygen partial pressure decreases to low values, Cr₂O₃ starts to reduce and oxidizes Si and C. Because the transport of oxides was not assumed to be limiting, the thermodynamic equilibrium defines the selectivity; this can be seen as the constant selectivity levels at the left side of the Figs. 2 and 3. Employing the same boundary conditions, the approach presented by Wei and Zhu¹⁷⁾ yields entirely different results, as seen in Figs. 2 and 3, selectivity values given in the text boxes. The calculated selectivity in the Wei and Zhu¹⁷⁾ method is not dependent on the mass transfer rate and only weakly on the composition too. It should be noted that rate coefficients are not the same for all reactions for reaching the same affinity level, as shown in Figs. 4 and 5. The sharp drop in the curve of the Cr oxidation reaction results from the onset of the reduction and changing the direction of the reaction reverse when delivery rate of O₂ is limited. In the analytical expressions (13) and (14), only the value k_f → ∞ resulted in the equilibrium for all reactions. It is important to understand that large, but finite values of k_f are required.

3.2. Effect of Residual Affinity

The effect of residual affinity A on the reaction rates, was also studied. The following parameters were employed in the calculations:β_L = 5×10⁻⁴ m/s, 1873 K (1600°C), x_C = 0.02, x_Cr = 0.17 and x_Si = 0.002. Figures 6 and 7 show the calculated rate coefficients required to reach pre-defined affinities for parallel oxidation reactions of Si, Cr and C. The studied range of residual affinities is 10⁻⁵ ≤ A ≤ 10⁴ J/mol. It should be noted that the upper range of values is considerably in excess of the general requirement suggested for Gibbs free energy minimization methods, which is typically 10^–3 J/mol.¹²⁾ Figures 8 and 9 present the effect of the residual affinity on the actual reaction rates. More specifically, it is illustrated which value of the residual affinity required to obtain a sufficiently high accuracy in terms of the reaction rates.

Fig. 6.

Rate coefficients required for achieving different residual affinities, β_G=2 m/s.

Fig. 7.

Rate coefficients required for achieving different residual affinities, β_G=0.02 m/s.

Fig. 8.

The effect of residual affinity on reaction rates, β_G=2 m/s.

Fig. 9.

The effect of residual affinity on reaction rates, β_G=0.02 m/s.

As can be seen, the required residual affinity i.e. the accuracy requirement has a considerable effect on the required rate coefficient values. The closer the solution is to the equilibrium, the larger are the rate coefficient values. The reverse, of course, is also true. Numerically, the higher the value of k_f is the more difficult it is to obtain convergence and numerical solution. Therefore, in order to have to have highest possible computational efficiency, the maximum values of the affinity that give correct results should be used.

The results shown in Figs. 6, 7, 8, 9 leave room for some interesting observations. In the case of the highest gas side mass transfer rate, Fig. 8, the value of residual affinity has no effect on the reaction rates within the studied residual affinity range of 10^–5…10⁴ J/mol. The overall rate is limited entirely by the liquid side mass transfer and consequently, the concentrations of oxidizing species are very low and there is plenty of O₂ available on the surface. As noted earlier, the rate mechanism is not controlled by the equilibrium and the liquid side mass transfer defines the reaction rates. The component having the highest mole fraction or the highest rate of mass transfer gets most of the oxygen. When the gas side mass transfer begins to limit the rate, see Fig. 9, this means that the maximum O₂ delivery rate approaches to the mass transfer rate of minor species, here C and Si. These simple examples given above show that in order to have a proper solution for the case of parallel gas-liquid oxidation reactions, it is necessary to consider both mass transfer and restricted chemical equilibrium. Because the mass transfer conditions may vary during different process stages in steelmaking, it is more preferable to use a more general model that incorporates all the relevant mass transfer and reaction mechanisms instead of predefining the rate-determining step, as suggested also by Oeters.¹⁸⁾

3.3. Applications and Comparison to Existing Approaches

The main advantage of the employed approach is that the local reaction rates and the equilibrium composition are solved simultaneously. Moreover, it provides a convenient approach to define the source term of chemical reactions in relation to volumetric unit. This is advantageous for mathematical and computational fluid dynamics (CFD) based models, which employ the control volume method to solve the conservation equations for momentum, mass and energy. In comparison to traditional Gibbs energy solvers, an additional benefit is that the source terms of species can be differentiated with respect to the mass transfer rates and the composition of the fluid flow, increasing the efficiency of the Newton’s method.

The method discussed also in this paper has already been employed in different cases for defining the mass transfer limited rates of reactions based on the law of mass action. The studied applications in the AOD process include mathematical modelling of reactions in a single gas bubble in liquid steel,¹⁹⁾ during side-blowing decarburization^20,21) and during the reduction of slag.^22,23) The method has also been applied for mathematical modelling of chemical heating in the CAS-OB process.²⁴⁾

The traditional equilibrium coefficient method would be to insert these activities as the boundary condition at the reaction surface and then solve the non-linear system of equations. In practice, this becomes very challenging as there is no rate expression in the system that could be differentiated and used as a method to stabilize the convergence of the solution. For CFD applications, the equilibrium coefficient method is not acceptable as the system on the surface should be somehow interlinked to flow field composition solution and the rates should be properly differentiated. We will discuss below, how to implement the LMA method to a practical problem and how to avoid infinite values for rate coefficients that would result in numerical overflow.

The most common approach to consider mass transfer controlled chemical reactions in metallurgical applications, and especially in the steel converters, is to simply derive the mass transfer limited local rate for some selected geometry, see Figs. 1(a), 1(b), 1(c), and then use the equilibrium composition, solved from equilibrium constant, at the reaction front as a boundary condition. A typical example could be the oxidation of carbon particles, by the reaction 2C + O₂ ⇌ 2CO.⁴⁾ The local mass transfer controlled oxidation rate of carbon particles $R_{\mathrm{C}}^{″}$ (in mol/(m²·s)) can be obtained simply as follows.   

$$\begin{array}{l}{R}_{\mathrm{C}}^{″}=2\frac{{\beta }_{\text{G}}{p}_{\text{G}}}{RT}\left(x_{{\text{O}}_2}^{\infty }-x_{{\text{O}}_2}^s\right)=2\frac{{\beta }_{\text{G}}{p}_{\text{G}}}{RT}\left(x_{{\mathrm{O}}_2}^{\infty }-x_{{\text{O}}_2}^{eq}\right)\\ \hspace{1em}\hspace{0.17em}=2\frac{{\beta }_{\text{G}}{p}_{\text{G}}}{RT}\left(x_{{\text{O}}_2}^{\infty }-\frac{p_{\text{CO}}^2}{K}\right)\cong 2\frac{{\beta }_G{p}_{\text{G}}}{RT}{x}_{{\text{O}}_2}^{\infty }\end{array}$$(37)

where x is the mole fraction. Overall rate can be obtained just by multiplying Eq. (37) by the specific surface area S_a [m²/m³] of carbon particles within the control volume. This is very useful approach for single reactions or simple reaction systems and can be easily implemented into process simulators or even to CFD.^25,26) However, if the system involves multiple and parallel reactions and the rate limiting step is not known beforehand or it varies during the process, this approach is not viable as the overall model would become a vast combination conditional of expressions. Solving the equilibrium composition directly from an equilibrium constant based non-linear system of equations is also complicated, because the system is stiff and the convergence is difficult to obtain.¹²⁾

To overcome the challenges in heavy computational CFD problems, Lu and Pope²⁷⁾ and Ren and Goldin²⁸⁾ developed the in situ adaptive tabulation (ISAT) method to speed up extremely time consuming computational calculations of reactive flows. As an example, they studied the possibility to speed up the solution of the chemical equilibrium in connection to the CFD simulation. They used a relaxation model wherein the rate is assumed to be proportional to the difference between the actual local and equilibrium concentrations, i.e. $R_i^{″}$ ~(c_i – c_i,eq). In general, it is necessary to solve the equilibrium composition using restricted equilibrium Gibbs free energy solvers, which is time consuming. However, the equilibrium composition at the surface is not the true equilibrium composition, but instead a restricted equilibrium for which, the traditional methods are not suitable. They were able to use the ISAT method successively and showed that the calculation time was reduced by a factor of hundreds.

A short review of other methods comparable to the LMA method is provided in the following.

3.3.1. Constrained Gibbs Free Energy Minimization

Koukkari et al.^29,30) proposed a method for calculating the constrained chemical equilibrium by means of introducing rate constraints to traditional Gibbs energy minimization routine. One possibility is to use constraints for rate-controlled (slow) reactions, which then leads to an effective combined thermodynamic-kinetic algorithm for the calculation of the chemical state of the system. This method is useful in process simulations, where the input and the output flow rates can be described without diffusion or a more complex mass transfer effects. In the applications studied in this study, the starting point is the species conservation equations such as Eq. (37) having a complex interlinked equations for chemical rate, mass transfer, diffusion, convection, turbulence effects and all these often in three dimensions. For this approach, explicit and differentiable rate expressions are required.

3.3.2. Affinity Based Selectivity Method

Wei and Zhu¹⁷⁾ calculated parallel competing reactions in AOD converter during combined top- and side-blowing and assumed that the O₂ is distributed to different dissolved elements, such as Si, Cr and C in linearly to chemical affinities, A_k = $-\mathrm{\Delta }{G}_k^{\circ }$ of oxidation reactions. Positive sign refers to reactions going from left to right. The main assumption in this method is that all the species on the liquid side have equal availability and the selectivity is linearly dependent on the affinity of the reactions, i.e. their thermodynamic driving force. According to this method, O₂ is distributed to three parallel oxidation reactions of Si, Cr and C with fractions Γ_k:   

$$\text{Si}+{\text{O}}_2\rightleftharpoons {\text{SiO}}_2\hspace{1em}{\mathrm{\Gamma }}_0=\frac{A_0}{A_0+1/0.75A_1+1/0.5A_2}$$(38)

  

$$\text{Cr}+0.75{\text{O}}_2\rightleftharpoons 0.5{\text{Cr}}_2{\text{O}}_3\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }{\mathrm{\Gamma }}_1=\frac{1/0.75A_1}{A_0+1/0.75A_1+1/0.5A_2}$$(39)

  

$$\text{C}+0.5{\text{O}}_2\rightleftharpoons \text{CO}\hspace{1em}{\mathrm{\Gamma }}_2=\frac{1/0.5A_2}{A_0+1/0.75A_1+1/0.5A_2}$$(40)

where the affinities A₀...A_k are calculated by the Eqs. (41), (42), (43) according to the bulk composition.   

$$-A_0=\mathrm{\Delta }{G}_0^{\circ }+RTln\left(\frac{a_{{\text{SiO}}_2}}{a_{\text{Si}}{p}_{{\text{O}}_2}^{\infty }}\right)$$(41)

  

$$-A_1=\mathrm{\Delta }{G}_1^{\circ }+RTln\left(\frac{a_{{\text{Cr}}_2{\text{O}}_3}^{0.5}}{a_{\text{Cr}}{p}_{{\text{O}}_2}^{0.75}}\right)$$(42)

  

$$-A_2=\mathrm{\Delta }{G}_2^{\circ }+RTln\left(\frac{a_{\text{CO}}}{a_{\text{C}}{p}_{{\text{O}}_2}^{0.5}}\right)$$(43)

Generally, the term in the parenthesis is referred to as the reaction quotient, Q_k. When the reaction reaches equilibrium, Q_k = K_k, and the affinity is equal to zero.¹⁶⁾ The drawback of this method is that it does not define the reaction rate directly. For this reason, it is not suitable to be used with the control volume method as such. The advantage of this approach is its simplicity and that it does not require a separate solver for the equilibrium composition. Unfortunately, no experimental validation was given for most of the components, and the detailed derivation of the model was not presented.

3.3.3. Constrained Solubility Method

Ding et al.³¹⁾ proposed a mathematical model for the vacuum oxygen decarburization (VOD) process. Reactions were assumed to take place in two reaction zones (metal-gas zone and metal-slag zone), while the volumetric flow rate of species from and to the reaction zones was estimated based on the plant data. Two different oxygen distribution models were investigated for the metal-gas zone. In both models, every time step Δt is divided into n smaller time steps Δt_s and the calculation proceeds so that only the reaction with the largest affinity takes place and consumes all the oxygen available during the time step Δt_s. The procedure is repeated in an iterative manner until the Gibbs free energy change of all reactions is zero or positive.

3.3.4. Effective Equilibrium Constant Method

The effective equilibrium constant method, also known as coupled reaction model, has been applied extensively for mathematical modelling of hot metal dephosphorization,^32,33,34) desiliconization³⁵⁾ and decarburization.^36,37,38) In this approach, the equilibrium constants are modified to effective equilibrium constants, which in combination with electro-neutrality condition, are employed to solve a system of parallel mass transfer limited reactions.

3.3.5. Coupled Gibbs Free Energy Minimization and Control Volume Method

Ersson et al.³⁹⁾ and Andersson et al.⁴⁰⁾ coupled CFD models of BOF and AOD converters, respectively, with computational thermodynamics. The local control volume equilibrium composition is achieved from the present chemical species after each time step. The obtained rate of phase mass change from this is then used as a constant rate for solving the rate at the new time step. Therefore, in addition to solution of the flow field, equilibrium calculations have to be carried at every time step. In order to save computational effort, equilibrium calculations were carried out only in the cells that have more than one phase present. One major advantage of this approach is that it does not require information on the interfacial surface area as the whole cell is set to equilibrium if an interphase exists. On the other hand, a drawback of the method is that the size of the computational cells at the reaction front should be infinitesimal to describe the mass transfer correctly.

3.3.6. Incremental Step Method

Jalkanen⁴¹⁾ and Virrankoski et al.⁴²⁾ presented a mathematical model for the basic oxygen furnace (BOF) process. The principal approach in their model is that the supplied oxygen is used momentarily only by a single reaction that has the highest affinity, as defined by Eqs. (41), (42), (43). Gaseous oxygen is supplied to the reaction surface by small step amounts and the composition of the reaction surface is constantly updated. When one species is rapidly depleted from the surface, its affinity decreases and in the new situation, another species that has then the highest affinity will continue to consume O₂ and so forth. The model has been reported to be very sensitive to the employed parameters, especially to the amount and size of the O₂ steps supplied to the system.⁴²⁾ In general, a numerical model should always be well posed and should converge to same final result independently on the numerical model parameters.^25,26) Table 2 summarizes the main features of all models.

Table 2. Comparison of methods to calculate parallel oxidation reactions.

Method	References	Applications	Advantages	Drawbacks
Effective equilibrium constant method
The concentration changes of elements are determined based on effective equilibrium constant and electrical neutrality equation.	32,33, 34,35, 36,37, 38)	Parallel oxidation and reduction reactions in the Hot Metal Dephosphorization and BOF processes	Transparent, easy to implement, suitable for on-line applications.
Constrained Gibbs free energy minimization method
Kinetic control considered through additional constraints in the Gibbs free energy minimization routine.	29,30)	Numerous applications, including combustion, chemical engineering and metallurgy.	Commercial software available, combined with commercial library on thermodynamic properties, derived from physical principles, well-suited for process simulators.	Difficult implementation to CFD or complex mass transfer limited systems with diffusion, convection and turbulence effects. Unsuitable for on-line applications.
Affinity based selectivity method
Selectivity of oxygen corresponds to the oxygen affinity of species, additional constraints for mass transfer limitation.	17,43, 44,45, 46,47)	Parallel oxidation and reduction reactions in the AOD process.	Transparent, easy to implement, suitable for on-line applications.	Cannot be derived from physical principles.
Constrained solubility method
Dissolved oxygen content is determined based on its solubility limit. Thereafter, the remaining oxygen is distributed so that the reaction with the largest oxygen affinity consumes all gaseous oxygen during a sub-time step.	31)	Parallel oxidation and reduction reactions in the VOD process.	Transparent, easy to implement, suitable for on-line applications.	Cannot be derived from physical principles. Solution is sensitive to model parameters, such as volumetric flow rates into and from the reaction zones.
Coupled Gibbs free energy minimization and volume element method
Volume cells with multiple phases are set into equilibrium with traditional Gibbs free energy minimization routine.	39,40)	Parallel oxidation and reaction reactions in the AOD and BOF processes.	State of the art combination of CFD and equilibrium calculations, uses commercial software in CFD and equilibrium solution.	Treatment of the interface in unclear, grid size dependency of the solution, large computational expense prevents on-line applications.
Incremental step method
Time steps are divided into sub-steps, during which the reaction with the largest affinity consumes all available oxygen	41,42)	Parallel oxidation and reduction reactions in the BOF process.	Transparent, easy to implement, possibly suitable for online applications.	Solution is sensitive for model parameters and cannot be derived from physical principles. The implementation includes non-physical tuning parameters.
Law of Mass Action based method
Rate expressions are calculated from the law of mass action, predefined residual affinity or equilibrium number determines unknown forward rate coefficients in an iterative procedure.	19,20, 21,22, 23,24)	Parallel oxidation and reduction reactions in the AOD and CAS-OB processes.	Derived from physical principles, non-sensitive to model parameters, easy to implement, simultaneous solution of the mass transfer limited rate and equilibrium composition, suitable for online applications.	No experience on suitability for CFD applications.

4. Conclusions

The objective of this paper was to present a novel law of mass action law based rate method for modeling parallel reversible mass transfer limited reactions in metallurgical systems. At first, we discussed the available methods from open literature for comparison. Then, we re-derived the new method for n parallel reactions in general form and discuss its practical implementation.

As a numerical example, a simple reaction model was derived for the case of a liquid steel surface exposed to O₂–CO gas mixture. The model was employed for studying the parallel oxidation of Si, Cr and C under conditions corresponding to those of the AOD process. To be transparent, the simple model developed here included only the conservation equations of all species wherein the rate expressions were based on a modified law of mass action. As an additional requirement, the affinities of the reactions at the reaction front were formulated such that they reach a pre-defined value, referred to as the residual affinity. This parameter assures that all the reactions are at the same distance from the constrained equilibrium. Our simulations showed that when the residual affinity is smaller than 10 J/mol, a sufficient degree of accuracy for the equilibrium is obtained. As a practical matter, a residual affinity of A = 0.001 J/mol was employed for all reactions in the final calculations.

Despite the relative simplicity of the employed model, the obtained results provide an in-depth analysis on the oxidation kinetics of molten stainless steel. In the case of the highest gas side mass transfer the value of residual affinity has virtually no effect of the reaction rates, within the range of 10^–5…10⁴ J/mol. This case is always limited by the liquid side mass transfer and there is plenty of residual O₂ available on the surface. The rate mechanism is not controlled by the equilibrium and the liquid side mass transfer defines the reaction mechanisms. The component having the highest mole fraction or the highest mass transfer rate consumes most of the oxygen.

When the gas side mass transfer begins to control the rates, partial pressure of O₂ at the reaction surface starts to decrease. Then, the role of the surface equilibrium becomes important and the species start to compete for the oxygen. Under such conditions the oxides of species with lower oxygen affinity will be reduced by species with higher oxygen affinity.

The results obtained in this work show that in order to have a proper solution for the case of parallel gas liquid oxidation reactions, it is an absolute requirement to consider both the mass transfer and the restricted chemical equilibrium. As the conditions vary during different process stages in steelmaking, it is preferable to employ a model that has all the relevant reaction steps included. The proposed new method is a transparent and direct solution of the reaction rates and therefore, it is in principle possible to couple the method with process simulators and CFD software.

Acknowledgements

This research is part of the FIMECC SIMP, a research program coordinated by the Finnish Metals and Engineering Competence Cluster (FIMECC). SSAB Europe Oy, Outokumpu Stainless Oy and the Finnish Funding Agency for Technology (TEKES) are acknowledged for funding this work. The financial support of the Academy of Finland (projects 258319 and 26495), Graduate School in Chemical Engineering, Finnish Foundation for Technology Promotion and Tauno Tönning Foundation is gratefully acknowledged. Ari Kruskopf is acknowledged for valuable discussions regarding this work.

List of symbols

A_k: Thermodynamic affinity of reaction k [J/mol]

A: Residual affinity [J/mol], frequency factor

a: Activity

c: Molar concentration [mol/m³]

E_a: Activation energy, [J/mol]

f: Residual

ΔG°: Change in standard Gibbs free energy [J/mol]

ΔH^*: Enthalpy of activation [J/mol]

h: Planck constant [J·s]

J: Jacobian matrix

k_b: Bolzmann constant [J/K]

k_f: Forward reaction rate coefficient [mol/(m²·s)]

K: Equilibrium constant

M: Molar mass [kg/mol]

p: Partial pressure [atm]

p_G: Total pressure [Pa]

R″: Reaction rate [mol/(m²·s)]

R: Gas constant [J/(mol·K)]

S_a: Specific surface area [m²/m³]

ΔG°: Change in standard Gibbs free energy [J/mol]

ΔS^*: entropy of activation [J/molK]

T: Absolute temperature [K]

t: Time [s]

u: Velocity [m/s]

x: Mole fraction

y: Vector of unknown variables

z: Axial coordinate [m]

Greek symbols

β: Mass transfer coefficient [m/s]

γ: Raoultian activity coefficient

ρ: Density [kg/m³]

ν_i,k: Stoichiometric coefficient of species i in reaction k

Subscripts

G: Gas phase

L: Liquid metal phase

s: Surface

S: Slag phase

∞: Bulk condition

Indices

i: Species

k: Reaction

n: Number of species i

r: Number of reactions k

References

• 1)   K.  Koch,  F.  Wolfdietrich and  P.  Valentin: Arch. Eisenhüttenwes., 47 (1976), 659.
• 2)   K.  Koch,  F.  Wolfdietrich and  P.  Valentin: Arch. Eisenhüttenwes., 49 (1978), 163.
• 3)   K.  Koch,  F.  Wolfdietrich and  P.  Valentin: Arch. Eisenhüttenwes., 49 (1978), 109.
• 4)   G. L.  Borman and  K. W.  Ragland: Combustion Engineering, McGraw-Hill International Editions, Singapore, (1998), 199.
• 5)   O.  Levenspiel: Chemical Reaction Engineering, 2nd ed., John Wiley & Sons, Inc., New York, (1972), 578.
• 6)   F.  Oeters: Kinetik metallurgischer Vorgänge bei der Stahlherstellung, ed. by W. Dahl, K. W. Lange and D. Papamantellos, Verlag Stahleisen m.b.H., Düsseldorf, (1972), 121.
• 7)   F.  Oeters: Metallurgy of Steelmaking, Verlag Stahleisen mbH, Düsseldorf, (1994), 512.
• 8)   R. I. L.  Guthrie: Engineering in Process Metallurgy, Clarendon Press, Oxford, UK, (1989).
• 9)   J.  Riipi,  T.  Fabritius, E.-.P. Heikkinen, P. Kupari and A. Kärnä: ISIJ Int., 49 (2009), 1468.
• 10)   W.  Nernst: Z. physik. Chem., 47 (1904), 52.
• 11)   R.  Taylor and  R.  Krishna: Multicomponent Mass Transfer, John Wiley & Sons, Inc., New York, (1993), 579.
• 12)   Y. K.  Rao: Stoichiometry and Thermodynamics of Metallurgical Processes, Cambridge University Press, Cambridge, UK, (1985), 957.
• 13)   P.  Waage and  C. M.  Guldberg: Forhandlinger i Videnskabs-Selskabet i Christiania, (1864), 35.
• 14)   P.  Waage and  C. M.  Guldberg: J. Chem. Educ., 63 (1986), 1044.
• 15)   K. K.  Kuo: Principles of Combustion, John Wiley & Sons, Singapore, (1986), 810.
• 16)   R. H.  Petrucci,  F. G.  Herring,  J. D.  Madura and  C.  Bissonnette: General Chemistry: Principles and Applications, Pearson, Toronto, (2011), 670.
• 17)  J.-.H. Wei and D.-.P. Zhu: Metall. Mater. Trans. B, 33 (2002), 111.
• 18)   F.  Oeters: Metallurgie der Stahlherstellung, Verlag Stahleisen mbH, Düsseldorf, (1989), 503.
• 19)   M.  Järvinen,  A.  Kärnä and  T.  Fabritius: Steel Res. Int., 80 (2009), 429.
• 20)   M. P.  Järvinen,  S.  Pisilä,  A.  Kärnä,  T.  Ikäheimonen,  P.  Kupari and  T.  Fabritius: Steel Res. Int., 82 (2011), 638.
• 21)   S. E.  Pisilä,  M. P.  Järvinen,  A.  Kärnä,  T.  Ikäheimonen,  T.  Fabritius and  P.  Kupari: Steel Res. Int., 82 (2011), 650.
• 22)  V.-.V. Visuri, M. Järvinen, P. Sulasalmi, E.-.P. Heikkinen, J. Savolainen and T. Fabritius: ISIJ Int., 53 (2013), 603.
• 23)  V.-.V. Visuri, M. Järvinen, J. Savolainen, P. Sulasalmi, E.-.P. Heikkinen and T. Fabritius: ISIJ Int., 53 (2013), 613.
• 24)   M.  Järvinen,  A.  Kärnä, V.-.V. Visuri, P. Sulasalmi, E.-.P. Heikkinen, K. Pääskylä, C. De Blasio, S. Ollila and T. Fabritius: ISIJ Int., 54 (2014), 2263.
• 25)   S. V.  Patankar: Numerical heat transfer and fluid flow, Hemisphere publishing, New York, NY, (1980), 150.
• 26)   J. H.  Ferziger and  M.  Peric: Computational Methods for Fluid Dynamics, Springer, Berlin, (1996), 113.
• 27)   L.  Lu and  S. B.  Pope: J. Comput. Phys., 228 (2009), 361.
• 28)   Z.  Ren and  G. M.  Goldin: Combust. Flame, 158 (2011), 1977.
• 29)   P.  Koukkari,  R.  Pajarre and  K.  Hack: Z. Metallkd., 92 (2001), 1151.
• 30)   R.  Pajarre,  P.  Blomberg and  P.  Koukkari: Comput. Aided Chem. Eng., 25 (2008), 883.
• 31)   R.  Ding,  B.  Blanpain,  P. T.  Jones and  P.  Wollants: Metall. Mater. Trans. B, 31 (2000), 197.
• 32)   S.  Ohguchi,  D. G. C.  Robertson,  B.  Deo,  P.  Grieveson and  J. H. E.  Jeffes: Ironmaking Steelmaking, 11 (1984), 202.
• 33)  S.-.Y. Kitamura, T. Kitamura, E. Aida, R. Sakomura, T. Kaneko and T. Nuibe: ISIJ Int., 31 (1991), 1329.
• 34)  S.-.Y. Kitamura, T. Kitamura, K. Shibata, Y. Mizukami, S. Kukawa and J. Nakagawa: ISIJ Int., 31 (1991), 1322.
• 35)   M.  Ishikawa: ISIJ Int., 44 (2004), 316.
• 36)   F.  Pahlevani,  S.  Kitamura,  H.  Shibata and  N.  Maruoka: Steel Res. Int., 81 (2010), 617.
• 37)   Y.  Lytvynyuk,  J.  Schenk,  M.  Hiebler and  A.  Sormann: Steel Res. Int., 85 (2014), 537.
• 38)   Y.  Lytvynyuk,  J.  Schenk,  M.  Hiebler and  A.  Sormann: Steel Res. Int., 85 (2014), 544.
• 39)   M.  Ersson,  L.  Höglund,  A.  Tilliander,  L.  Jonsson and  P.  Jönsson: ISIJ Int., 48 (2008), 147.
• 40)   N. Å. I.  Andersson,  A.  Tilliander,  L. T. I.  Jonsson and  P. G.  Jönsson: Ironmaking Steelmaking, 40 (2013), 390.
• 41)   H.  Jalkanen: Sohn Int. Symp., Advanced Processing of Metals and Materials, Vol. 2 – Thermo and Physicochemical Principles: Iron and Steelmaking, The Minerals, Metals & Materials Society, Warrendale, PA, (2006), 541.
• 42)   R.  Virrankoski,  T.  Oksanen and  H.  Jalkanen: CONSIM 4.0 – Simulation Software for Oxygen Converter. TKK-MK-121, Helsinki University of Technology, Espoo, Finland, (2001).
• 43)   P.  Sjöberg: Doctor thesis, Royal Institute of Technology, Stockholm, (1994), 157.
• 44)  J.-.H. Wei and D.-.P. Zhu: Metall. Mater. Trans. B, 33 (2002), 121.
• 45)  H.-.L. Zhu, J.-.H. Wei, G.-.M. Shi, J.-.H. Shu, Q.-.Y. Jiang and H.- B. Chi: Steel Res. Int., 78 (2007), 305.
• 46)  G.-.M. Shi, J.-.H. Wei, H.-.L. Zhu, J.-.H. Shu, Q.-.Y. Jiang and H.-.B. Chi: Steel Res. Int., 78 (2007), 311.
• 47)  J.-.H. Wei, Y. Cao, H.-.L. Zhu and H.-.B. Chi: ISIJ Int., 51 (2011), 365.
• 48)   A.  Chatterjee,  D. H.  Wakelin and  A. V.  Bradshaw: Metall. Trans., 3 (1972), 3167.