**2. Model Establishment**

The geometric model of a blast furnace was established based on blast furnace B in Bayi Steel of Baowu Group in China, and it was a one-sixth (60 degree section) sector of the blast furnace. The total height of the model was 27.43 m, as shown in Figure 1, where, for the sake of simplification, the part below the taphole in the blast furnace hearth is not included.

**Figure 1.** Geometric and mesh models and volume fraction distribution of solid phase of blast furnace B in Bayi Steel of Baowu Group (**a**) geometric model; (**b**) mesh model; (**c**) volume fraction distribution of solid phase.

Both gas and solid phases were considered. The gas phases were the same as those in practice, including CO, CO2, H2, H2O and N2. The density of gas phase was calculated by the ideal gas law. The solid phase was simplified to include Fe2O3, Fe3O4, FeO, Fe, and C. Other components were not considered for the sake of simplicity. The density of solid phase was 4190 kg·m−3, the same as the apparent density of the original solid material. Both gas phase and solid phase were considered as continuous phases using the Eulerian method. The mass, energy, and species transfer can be described by Equation (1) under steady state [20,21].

$$\nabla(a\_{\mathcal{V}} \cdot \rho\_{\mathcal{V}} \cdot \phi \cdot \overrightarrow{V}\_{\mathcal{V}}) = \nabla(a\_{\mathcal{V}} \cdot \Gamma \cdot \nabla(\phi)) + \mathcal{S}\_{\phi} \tag{1}$$

where ρ is the phase. Γ and *S* are the effective diffusivity and the source respectively, varying with respect to the different variables φ as listed in Table 1 [21–25].


**Table 1.** Parameters in Equation (1).

*Rn* refers to different reactions.

$$\overline{\overline{\tau}}\_{p} = \varepsilon\_{p} \cdot \mu\_{p} \cdot \left[\nabla \cdot \overrightarrow{\overline{\upsilon}}\_{p} + \left(\nabla \cdot \overrightarrow{\overline{\upsilon}}\_{p}\right)^{T}\right] - \frac{2}{3} \cdot \varepsilon\_{p} \cdot \mu\_{p} \cdot \left(\nabla \cdot \overrightarrow{\overline{\upsilon}}\_{p}\right) \cdot \overline{\overline{I}}$$

$$\overrightarrow{F}\_{\mathcal{S}^{s}} = -[150 \cdot \frac{\left(1 - \varepsilon\_{\mathcal{S}}\right)^{2} \cdot \mu\_{\mathcal{S}}}{\varepsilon\_{\mathcal{s}}^{3} \cdot d\_{\mathcal{s}}^{2}} + 1.75 \cdot \frac{\rho\_{\mathcal{S}} \cdot \varepsilon\_{\mathcal{s}} \cdot \left|\overrightarrow{\overline{\upsilon}}\_{\mathcal{S}} - \overrightarrow{\overline{\upsilon}}\_{\mathcal{S}}\right|}{d\_{\mathcal{s}}}] \cdot \left(\overrightarrow{\overline{\upsilon}}\_{s} - \overrightarrow{\overline{\upsilon}}\_{\mathcal{S}}\right)$$

$$E\_{\mathcal{S}^{s}} = -\frac{6 \cdot k\_{\mathcal{S}} \cdot \varepsilon\_{\mathcal{S}} \cdot \varepsilon\_{\mathcal{s}}}{d\_{\mathcal{s}}^{2}} \cdot \left(2.0 + 0.6 \cdot \text{Re}\_{\text{s}}^{1/2} \cdot \text{Pr}\_{\text{\mathcal{S}}}^{1/3}\right) \cdot \left(T\_{\mathcal{S}} - T\_{s}\right)$$

The chemical reactions between gas phase and solid phase are as follows, including the indirect reduction of iron oxide (Reactions 1–6), carbon solution loss reaction (Reactions 7), water gas reaction (Reaction 8), water gas shift reaction (Reaction 10), combustion reaction of C and H2 (Reaction 9 and 11), and so on.


The three-interface unreacted core model was used to calculate the chemical reaction rates, and the physical chemistry data were taken from Perry et al.'s book [26]. The reaction rate constants of the indirect reduction of iron ore by CO or H2, carbon solution loss reaction, water gas reaction, combustion of carbon and water gas shift reaction were taken from other work [27–30]. The combustion rate of H2 with O2 was taken from Kuwabara et al.'s work [31]. The effective diffusion coefficients were taken from other work [27–30]. The viscosity and thermal conductivity of gas were obtained from the literature [18,32–34]. The chemical reaction rates are as follows.

$$R\_1 = \frac{A\rho\_\mathcal{g}}{\mathcal{W}} \sum\_{m=1}^{\mathcal{G}} a\_{1,m} (K\_m \frac{w\_{\text{CO}\_2\mathcal{G}}}{M\_{\text{CO}}} - \frac{w\_{\text{CO}\_2\mathcal{G}}}{M\_{\text{CO}\_2}}) \tag{2}$$

$$R\_2 = \frac{A\rho\_\mathcal{g}}{W} \sum\_{m=1}^{\mathcal{G}} a\_{2,m} (K\_m \frac{w\_{\text{CO},\text{\(g\)}}}{M\_{\text{CO}}} - \frac{w\_{\text{CO}\_2,\text{\(g\)}}}{M\_{\text{CO}\_2}}) \tag{3}$$

$$R\_3 = \frac{A\rho\_\mathcal{g}}{W} \sum\_{m=1}^3 a\_{3,m} (K\_m \frac{w\_{\text{CO},\text{g}}}{M\_{\text{CO}}} - \frac{w\_{\text{CO}\_2,\text{g}}}{M\_{\text{CO}\_2}}) \tag{4}$$

$$R\_4 = \frac{A\rho\_\mathcal{J}}{W} \sum\_{m=4}^6 a\_{4,m} (K\_m \frac{w\_{\text{H}\_2\text{-}\mathcal{J}}}{M\_{\text{H}\_2}} - \frac{w\_{\text{H}\_2\text{O}\_\mathcal{J}}}{M\_{\text{H}\_2\text{O}}}) \tag{5}$$

$$R\_{\mathsf{F}} = (\frac{6\alpha\_{\mathsf{s}} \cdot \varepsilon\_{\mathsf{s}}}{\varphi\_{\mathsf{s}} \cdot d\_{\mathsf{s}}}) \frac{\rho\_{\mathsf{g}}}{\mathcal{W}} \sum\_{m=4}^{6} a\_{\mathsf{5},m} (K\_{m} \frac{w\_{\mathsf{H}\_{2}\mathsf{g}}}{M\_{\mathsf{H}\_{2}}} - \frac{w\_{\mathsf{H}\_{2}\mathsf{O}\_{\mathsf{A}}}}{M\_{\mathsf{H}\_{2}\mathsf{O}}}) \tag{6}$$

$$R\_6 = \frac{A\rho\_\mathcal{g}}{W} \sum\_{m=4}^6 a\_{6,m} (K\_m \frac{w\_{\text{H}\_2\text{g}}}{M\_{\text{H}\_2}} - \frac{w\_{\text{H}\_2\text{O}\_{\text{Al}}}}{M\_{\text{H}\_2\text{O}}}) \tag{7}$$

$$R\_{\heartsuit} = \frac{k\_{7,1} P\_{\text{CO}\_2} \omega\_{\text{C}} \rho\_{\text{s}} \varepsilon\_{\text{s}}}{1 + k\_{7,2} P\_{\text{CO}} + k \gamma\_{\text{\gamma}} P\_{\text{CO}\_2}} \tag{8}$$

$$R\_8 = \frac{k\_{8,4} P\_{\text{H}\_2\text{O}} \omega\_{\text{C}} \rho\_{s\text{E}s}}{1 + k\_{8,2} P\_{\text{CO}} + k\_{8,3} P\_{\text{CO}\_2} + k\_{8,5} P\_{\text{H}\_2\text{O}}} \tag{9}$$

$$R\_{\theta} = \frac{1}{1 + 2500 \exp\left(-\frac{12.400}{1.987 T\_{\odot}}\right)} \left[\frac{\varepsilon\_{\circ} \rho\_{\circ} \omega \nu\_{\bullet}}{M\_{\text{O}\_{2}}}\right] \left[\frac{d\_{s} \rho\_{s}}{AD\_{\text{O}\_{2}, \text{N}\_{2}}} \frac{1}{\text{S} \text{t}}\right]^{-1} \tag{10}$$

$$R\_{10} = \varepsilon\_5 (k\_{10,1} + K\_{10} k\_{10,2}) \left( P\_{\rm CO} P\_{\rm H\_2O} - P\_{\rm CO\_2} P\_{\rm H\_2} / K\_{10} \right) P\_{\rm H\_2O} \tag{11}$$

$$R\_{11} = \frac{\rho\_{\mathcal{S}} \alpha\_i}{M\_i} \left( \frac{1}{Ak\_{f11}} + \frac{1}{k\_{11,1}} \right)^{-1} \tag{12}$$

In the numerical solution of the arising differential equations, a structured grid was applied. The average mesh size is 40 mm. The sensitivity study of mesh size was carried out with average size of 100 mm, 80 mm, 60 mm, 40 mm, 30 mm, and 20 mm. The difference of top gas temperature between 60 mm and 40 mm is 4.9% while that between 40 mm and 30 mm is within 0.5%. This suggests that the mesh size of 40 mm is reasonable and confirms the mesh independence. The model contains 261,131 nodes and 248,496 hexahedral cells.

The assumptions for this model were as follows: (1) The powder phase was not considered; (2) Other chemical reactions such as flux decomposition and reduction of non-ferrous compounds were ignored; (3) The melting of solids was ignored.

The main operational parameters of the blast furnace in 2018 are listed in Table 2, where HM is the abbreviation of hot metal.


**Table 2.** Main operational parameters of the blast furnace.

The total Fe content in the mixed iron-bearing materials was 56.20%, and the total C contents in the coke and coal were 87.38% and 62.46%. Based on the material and energy balances, Fe in the mixed iron-bearing materials was converted to Fe2O3 and only C in the coke was considered. The coal injection and the oxygen of the blast at the tuyere were converted to CO. Therefore, the mass fraction of Fe2O3 and C were 77.05% and 22.95% at the top of blast furnace, respectively. The volume fraction of the solid phase was fixed, as shown in Figure 1c. The burden velocity distribution was firstly calculated with a simple model without taking into account heat and energy transfer. The results of burden velocity were then loaded to the present complex model taking into consideration all the transfers and reactions. The top gas pressure in the blast furnace was 194 kPa. The gas compositions at the tuyere were N2-74.92 vol%, O2-15.34 vol%, CO-9.15 vol%, H2O-0.59 vol%. The conservation equations were solved numerically by the finite volume method with commercial software ANSYS FLUENT (release 17.0) [35]. The Eulerian multiphase module was used in this model. The first order upwind scheme was used for discretization of density, momentum, volume fraction, energy, gas and solid species, and so on. The Phase Coupled SIMPLE method was applied [21,30]. The simulation was considered to have converged when the residuals for each variable were less than 10−5. The solution flow chart is shown in Figure 2.

**Figure 2.** Solution flow chart of the simulation.
