Mathematical Model of Solid Flow Behavior in a Real Dimension Blast Furnace

Abstract

A mathematical model based on the continuum mechanic concept has been developed to describe the profile of solid particles in an industrial scale blast furnace with respect to the in-furnace conditions and its characteristics such as the shape and size of the deadman. The Navier-Stokes differential equation for multi-phase multi-dimensional space has been used to describe the behavior of existing phases. The surface stress tensor has been defined as an extra term and added to the Navier-Stokes equation to describe the particle-particle interactions. This extra term in the Navier-Stokes equation behave as a breaking force when the particles are sliding down. It is shown that the particles change their profile from a V-shape to a W-shape due to the characteristics of the deadman. Moreover, the velocity magnitude is higher at the outer surface of the deadman for higher grid-slabs in this region than the near-wall cells. However, the situation changes as solid particles moving to even lower level of grid-slabs at the outer surface of the deadman in comparison to near-wall cells. It has also been shown that an increase in the magnitude of the effective pressure reduces the velocity magnitude of descending particles.

1. Introduction

The blast furnace vessel is a very complex operating systems in the field of the process metallurgy. Even though the blast furnace process was assumed to be at its end during 80s, the rate of production of the liquid iron by blast furnace increased dramatically during the last two decades.¹⁾ This increment eventually spotlighted the blast furnaces as the main alternative in the ironmaking root in comparison to direct reduced iron DRI and other methods. Therefore, it is important to have a good understanding of the process and phenomena inside this metallurgical vessel. Meanwhile, the operational conditions and mostly the high temperature make it if not impossible, extremely difficult to retrieve any data from the in-furnace conditions for full scale operational blast furnaces. Hence, alternative methods for predicting the condition of in-furnace processes became the focus of numerous studies. Improvement in numerical methods and developments of computer hardware simultaneously, brought attention to modeling of metallurgical processes. Consequently, mathematical models of blast furnaces have been developed to study and investigate the in-furnace phenomena.²⁾ Since early 90s, many mathematical models have been presented to illustrate the behavior of this vessel.^{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30)}

In general, there are two major approaches in modeling of a blast furnace: continuum and discrete models. The continuum approach assumes that all the phases in the furnace behave like a continuum medium and is based on a local average principle.^{3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21)} On the other hand, the discrete method is based on the Newton’s second law of motion and is focused on the behavior of each and all individual particles in the system.^{22,23,24,25,26,27,28,29,30)} The discrete element method (DEM) is originally developed by Cundall and Strack³¹⁾ and is very effective to model granular materials. However, it is still very time-expensive to use this method for simulations of a blast furnace, even considering the development of computers hardware. One of these limitations is the number of particles that can be handled by the computers, which hardly would be comparable to the number of particles in an actual blast furnace in today’s simulations.

The continuum method, as mentioned earlier, assumes that each phase in the system is a continuum medium. Therefore, the conservation equations can be used to model and predict the behavior of different phases which coexist in the furnace. Since there are more than one medium in this process, a multi-phase model must be used. Thus, the conservation equations for each and all phases must be solved simultaneously, whilst the effect of each phase on other phases in the process must be taken into account. One of the benefits of this method is the simplicity of expanding a two-dimensional model to a three-dimensional model.

In addition to the methods described above, there are other attempts to model a blast furnace. One of these methods is based on the hypo-plasticity theory. This theory was developed at the University of Karlsruhe.^32,33) However, there are not many reports on modeling a blast furnace based on the hypo-plasticity theory yet.^34,35,36)

This study is based on a continuum model. A three-dimensional conservation equation for a multi-phase system has been solved numerically with the finite volume method. The existing phases are gas and solid. This article presents and introduces the mathematical model used to describe solid particles behavior in an industrial scale blast furnace process. The geometry has been adopted from the blue print of SSAB’s blast furnace M4 in Oxelösund. Note that, due to the nature of the solid particles, an extra term has been added to the Navier-Stokes equation to model the surface stress tensor to describe the particle-particle interactions due to collision and packing. A comparison study with results reported⁸⁾ in the literature verifying the functionality is reported in a related paper.³⁷⁾ In the future there might also be a possibility to validate the model using actual data from the experimental furnace hosted by Swerea MEFOS AB.

2. Theory and Model Description

2.1. Governing Equations: General Conservation Equations

As it was mentioned earlier, one of the approaches in modeling a blast furnace process is to assume that all the phases in a blast furnace behave as a continues phase. Therefore, based on the above assumption, the blast furnace process can be modeled based on the continuum mechanics theory in an Eulerian frame of reference. In this study, two separate sets of Navier-Stokes equations are solved and linked together by interfacial forces. The general conservation equation in an Eulerian farame of reference can be written as follows:   

$$\frac{\partial \left({\mathit{\epsilon }}_{\text{i}}{\rho }_{\text{i}}{\mathrm{\Phi }}_{\text{i}}\right)}{\partial \text{t}}+\nabla .\left({\mathit{\epsilon }}_{\text{i}}{\rho }_{\text{i}}{\overrightarrow{\text{u}}}_{\text{i}}{\mathrm{\Phi }}_{\text{i}}\right)=\nabla .\left({\mathit{\epsilon }}_{\text{i}}{\mathrm{\Gamma }}_{{\mathrm{\Phi }}_{\text{i}}}\nabla {\mathrm{\Phi }}_{\text{i}}\right)+{\text{S}}_{{\mathrm{\Phi }}_{\text{i}}}$$(1)

where ε_i is the volume fraction, ρ_i is the density and ${\overrightarrow{\mathrm{u}}}_{\mathrm{i}}$ is the velocity vector for phase i. Moreover, the effective diffusive transfer coefficient, ${\mathrm{\Gamma }}_{{\mathrm{\Phi }}_{\mathrm{i}}}$ , is represented in different ways due to the dependent variable being considered. For example, in the velocity field, it represents the dynamic viscosity of the phase considered, μ_i.¹¹⁾ However, in order to adjust the reality of a solid phase discontinuity, Eq. (1) must be extended, so that it can describe the behavior of this phase more accurately.

Zhang et al.,⁸⁾ presented the following modified equations for the flow of solid particles in a blast furnace:   

$$\frac{\partial {\mathit{\epsilon }}_{\text{s}}}{\partial \text{t}}+\nabla .\left({\mathit{\epsilon }}_s{\overrightarrow{\text{u}}}_{\text{s}}\right)=0$$(2)

  

$$\frac{\partial \left({\mathit{\epsilon }}_s{\rho }_s{\overrightarrow{\text{u}}}_{\text{s}}\right)}{\partial \text{t}}+\nabla .\left({\mathit{\epsilon }}_s{\rho }_s{\overrightarrow{\text{u}}}_{\text{s}}{\overrightarrow{\text{u}}}_{\text{s}}\right)=\nabla .{\tau }_s+{\mathit{\epsilon }}_s{\rho }_s\overrightarrow{g}$$(3)

where $\overrightarrow{g}$ is the gravitational acceleration vector. In fact, particles can either behave as a fluid media or act as a stagnant fluid, providing a resistance to particle movement. Therefore, the surface stress τ_s in Eq. (3) may be considered to consist of two parts:^{38,39,40,41,42)}   

$${\tau }_s= {\tau }_{rd}+ {\tau }_{ri}$$(4)

where, as it is introduced by Zhang et al.,⁸⁾ τ_rd is the rate-dependent part of the surface stress and τ_ri is the rate-independent part. The variable τ_rd originates from the particle-particle kinetic interaction due to collision. Moreover, τ_ri originates from particle-particle contact interaction due to packing, including normal contact and shear friction.

In this study, the three-dimensional framework of Eqs. (2), (3), (4) for a system of two-phase gas-solid has been used which can be defined as follows:

• For the solid phase:   

$$\frac{\partial {\mathit{\epsilon }}_{\text{s}}}{\partial t}+\nabla .\left({\mathit{\epsilon }}_s{\overrightarrow{\text{u}}}_{\text{s}}\right)=0$$(5)

  

$$\begin{array}{l}\frac{\partial \left({\mathit{\epsilon }}_s{\rho }_s{\overrightarrow{\text{u}}}_{\text{s}}\right)}{\partial t}+\nabla .\left({\mathit{\epsilon }}_s{\rho }_s{\overrightarrow{\text{u}}}_{\text{s}}{\overrightarrow{\text{u}}}_{\text{s}}\right)=-{\mathit{\epsilon }}_s\nabla p_s-\beta \left({\overrightarrow{\text{u}}}_{\text{s}}-{\overrightarrow{\text{u}}}_g\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\hspace{0.17em}+\nabla .{\tau }_{rd}+\nabla .{\tau }_{ri}+{\mathit{\epsilon }}_s{\rho }_s\overrightarrow{g}\end{array}$$(6)

• For the gas phase:   

$$\frac{\partial {\mathit{\epsilon }}_{\text{g}}}{\partial t}+\nabla .\left({\mathit{\epsilon }}_g{\overrightarrow{\text{u}}}_{\text{g}}\right)=0$$(7)

  

$$\begin{array}{l}\frac{\partial \left({\mathit{\epsilon }}_g{\rho }_g{\overrightarrow{\text{u}}}_{\text{g}}\right)}{\partial t}+\nabla .\left({\mathit{\epsilon }}_g{\rho }_g{\overrightarrow{\text{u}}}_g{\overrightarrow{\text{u}}}_g\right)=-{\mathit{\epsilon }}_g\nabla p_g-\beta \left({\overrightarrow{\text{u}}}_{\text{g}}-{\overrightarrow{\text{u}}}_s\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} } \mathrm{\hspace{0.17em} }+\nabla .{\tau }_g+{\mathit{\epsilon }}_g{\rho }_g\overrightarrow{g}\end{array}$$(8)

where,   

$${\mathit{\epsilon }}_{\text{g}}+{\mathit{\epsilon }}_{\text{s}}=1$$(9)

Equation (9) represents the relationship between the volume fractions of different phases that exist in the system. Since the voidage of a bed in a real blast furnace can vary, the volume fractions should be considered as unknown parameters. In Eqs. (6) and (8), the parameter β represents the momentum transfer coefficient and is dependent on the gas volume fraction. For ε_g ≤ 0.8 then β should be evaluated based on the classical Ergun equation.⁴³⁾ On the other hand, for ε_g > 0.8 Richardson and Zaki’s formulation should be used.⁴⁴⁾ Hence, β could be represented as:   

$$\beta =\{\begin{array}{cc}150\frac{{\mathit{\epsilon }}_s^2}{{\mathit{\epsilon }}_{\text{g}}}\frac{{\mu }_{\text{g}}}{{\left({\phi }_s{d}_{\text{p}}\right)}^2}+1.75{\mathit{\epsilon }}_{\text{s}}\frac{{\rho }_{\text{g}}}{{\phi }_s{d}_{\text{p}}}\left|{\overrightarrow{\text{u}}}_{\text{g}}-{\overrightarrow{\text{u}}}_s\right|\hfill & {\mathit{\epsilon }}_{\text{g}}\le 0.8\hfill \\ \frac{3}{4}{C}_d\frac{{\mathit{\epsilon }}_g^{-1.65}{\mathit{\epsilon }}_{\text{s}}}{{\phi }_s{d}_{\text{p}}}{\rho }_{\text{g}}\left|{\overrightarrow{\text{u}}}_{\text{g}}-{\overrightarrow{\text{u}}}_s\right| \hfill & {\mathit{\epsilon }}_{\text{g}}>0.8\hfill \end{array}$$(10)

where the parameter φ_s represents the sphericity of the particle, d_p is the diameter and C_d is the drag coefficient and is related to the particle Reynolds number Re_p as follows:   

$$C_d=\{\begin{array}{cc} \frac{24}{{\text{Re}}_{\text{p}}}\left[1+0.15R{e}_p^{0.687}\right]\hfill & {\text{Re}}_{\text{p}}\le 1\mathrm{\hspace{0.17em} }000\hfill \\ 0.44 \hfill & {\text{Re}}_{\text{p}}>1\mathrm{\hspace{0.17em} }000\hfill \end{array}$$(11)

The particle Reynolds number in Eq. (11) is defined in the following manner:   

$$R{e}_p=\frac{{\mathit{\epsilon }}_{\text{g}}{\rho }_{\text{g}}\left|{\overrightarrow{\text{u}}}_{\text{g}}-{\overrightarrow{\text{u}}}_s\right|d_{\text{p}}}{{\mu }_{\text{g}}}$$(12)

2.2. Governing Constitutive Equations: Surface Stress

To be able to fully describe the above equations, it is necessary to determine the constitutive equations. The constitutive equations need to be presented for both phases, i.e. gas and solid, in the system. In addition, other related parameters like the solid viscosity, elastic modulus and friction coefficient should be determined. The constitutive equation for the gas phase is referred to as the viscous stress and it can be expressed as:   

$${\tau }_g=2{\mathit{\epsilon }}_{\text{g}}{\mu }_{\text{g}}{D}_{\text{g}}-\frac{2}{3}{\mathit{\epsilon }}_{\text{g}}{\mu }_{\text{g}}tr{D}_{\text{g}}I$$(13)

where μ_g is the gas viscosity, D_g is gas stretching tensor, trD_g is the trace of the stretching tensor and I is the identity matrix. This equation is based on the Newtonian expression.

As already stated, the solid surface stress consists of a rate-dependent and a rate-independent part. The rate-dependent part, similarly to the gas viscous stress equation, is the Newtonian part of the constitutive equation for the solid phase. Therefore, the viscous part of the solid surface stress equation is expressed as:   

$${\tau }_{rd}=2{\mathit{\epsilon }}_s{\mu }_s{D}_s-\frac{2}{3}{\mathit{\epsilon }}_s{\mu }_{s}tr{D}_{s}I$$(14)

where μ_s is the solid viscosity and D_s is the solid stretching tensor associated with the motion of particles. The definitions of the stretching tensor, identity matrix and trace of the stretching tensor for solid phase are given below:   

$$D_s^{i,j}=\frac{1}{2}\left(\frac{\partial u_s^i}{\partial x_j}+\frac{\partial u_s^j}{\partial x_i}\right)$$(15)

  

$$I=\left[\begin{array}{ccc}1& & \\ & 1& \\ & & 1\end{array}\right]$$(16)

and   

$$tr{D}_{\text{s}}={{\displaystyle \sum }}^{}{D}_s^{i,i}$$(17)

which are the same for the gas phase.

The rate-independent part of the solid surface equation is based on the Coulomb frictional relation and can be expressed as follows:   

$${\tau }_{ri}=-P_{e}I+2{\eta }_s\frac{D_s}{\sqrt{tr{D}_s^2}}$$(18)

where P_e is the effective pressure and η_s is the coefficient of plastic modulus. The variable η_s is based on Norem et al.:⁴⁵⁾   

$${\eta }_{\text{s}}=\frac{1}{\sqrt{2}}\left({\alpha }_1+{\alpha }_2{P}_e^m\right)$$(19)

where α₁ is equal to zero for cohesionless materials and m is equal to 1 corresponding to the Coulomb frictional model. In addition, the friction coefficient α₂ can be expressed as follow:   

$${\alpha }_2=\frac{\left|\frac{\partial w_s}{\partial x} \frac{1}{w_s}\right|}{{\left|\frac{\partial w_s}{\partial x} \frac{1}{w_s}\right|}_{Max}}tan\varphi$$(20)

where w_s is the components of the velocity vector for the solid phase in vertical direction, ϕ is the internal friction angle and Max is the maximum value in the flow field domain.

The effective pressure P_e, in Eqs. (18)–(19) is reported by Zhang et al. as follows and is employed according to the normal component of the surface stress due to Orr:⁴⁶⁾   

$$P_e=\frac{G_0}{c}{e}^{c\left({\mathit{\epsilon }}_{\text{s}}-{\mathit{\epsilon }}_s^{*}\right)}$$(21)

where G₀ is the normalizing unit factor, c is the compaction modulus and ${\mathit{\epsilon }}_s^{*}$ is the packing density, i.e. the solid volume fraction measured under the standard loose or poured packing condition. The parameters G₀ and c are constants which are reported as 1.0 Pa and 150 respectively by Zhang et al.⁸⁾ The packing density ${\mathit{\epsilon }}_s^{*}$ and the solid phase density are set to 0.73 and 3700 kg/m³, respectively.

The effective pressure P_e defined by Eq. (21) and used in Eqs. (18)–(19) simulates the deformation of the particles due to the loading of the furnace whilst p_s in Eq. (6) is the dynamic pressure associated with the Eulerian based calculation of the flow field.

2.3. Source Term

The source term ${\text{S}}_{{\mathrm{\Phi }}_{\text{i}}}$ in Eq. (1) is due to external forces such as body forces, which arise through the gravitational forces. Furthermore, ${\text{S}}_{{\mathrm{\Phi }}_{\text{i}}}$ is also due to inter-phase interactions which arise through chemical reactions^5,6,7,9) or surface interactions. Chen et al.⁴⁾ included the friction force between particles and wall or deadman to the source term in addition to the body force. These surface interactions can be due to the effect of gas on solid particles, which is formulated based on Ergun’s equation⁴³⁾ and Richardson and Zaki.⁴⁴⁾ Moreover, ${\text{S}}_{{\mathrm{\Phi }}_{\text{i}}}$ also includes solid-solid friction-collision which is described by constitutive equations by Zhang et al.⁸⁾ In this study, the source term is based on the body force, the effect of gas on solid particles, and particle-particle interaction.

2.4. Initial and Boundary Conditions

To solve any differential equation it is necessary and critical to specify appropriate boundary conditions. In fact, when a differential equation is to be solved, the solution would be highly dependent on the boundary conditions defined in the domain. Therefore, to be able to find a solution for Eqs. (5), (6), (7), (8), appropriate boundary conditions for the blast furnace process have to be identified and introduced to the solver.

In a blast furnace process simulation, the essential boundary conditions are actual furnace operational conditions. For example, the velocity of gas at the tuyeres, the top gas pressure, the solid mass flow at the top of the furnace which is called throat and the removal rate of solid particles at the bottom of the domain. Moreover, a blast furnace flow is a counter-current flow which means that the existing phases in the system travel in opposite directions. Therefore, the top layer acts as an inlet for the solid phase while it serves as an outlet for the gas phase. The solid outlet is located at the bottom of the domain so that the solid particles could leave the domain from the bottom. In summary, the following boundary conditions have been used in this study:

• A pressure boundary is used at the top of the domain to remove the gas.

• Solid particles are assumed to be charged homogeneously at the top of the domain with a velocity giving a mass flow rate of 90 kg/s.

• A solid velocity is set at the solid outlet based on the mass flow rate of the solid phase at the top.

In addition, the deadman in a blast furnace plays an important role. Therefore, the information on the characteristics of the deadman in any blast furnaces modeling is essential. In most of the blast furnace modeling, the characteristics of the deadman are set in the model. However, Zhang et al.⁸⁾ presented a method to define the boundary of the deadman based on the velocity of the solid phase in the system. Nevertheless, the method described in their study only considers the velocity of the particles. Also, the outlet area should be set in accordance with the inlet condition. In the authors study, the characteristics of the deadman e.g. shape, size and height are set based on a field study during a relining which was done on the SSAB blast furnace in Oxelösund, Sweden.⁴⁷⁾ Furthermore, no-slip wall functions has been used at the deadman surface.

The no-slip boundary condition cannot be applied at the impenetrable solid wall.^48,49) In fact, partial slip of solid particles at the wall is possible. Therefore, the solid tangential velocity at the wall is assumed to be proportional to its gradient at the wall. This is due to the fact that the particle diameter is usually larger than the length scale of the surface roughness of the rigid wall. Therefore, the solid tangential velocity at the wall can be expressed by:⁸⁾   

$${{\text{w}}_{\text{s}}|}_{\text{wall}}=-{\lambda }_{\text{p}}\frac{\partial {\text{w}}_{\text{s}}}{\partial \text{n}}|{ }_{\text{wall}}$$(22)

where w_s is the solid velocity in vertical direction, n is the normal direction to the wall and the slip parameter λ_p can be described by:^48,49)   

$${\lambda }_{\text{p}}=\frac{1}{{\mathit{\epsilon }}_{\text{s}}^{1/3}}{\text{d}}_{\text{p}}$$(23)

The λ_p denotes the mean distance between particles. This equation shows that for small particles, the boundary condition is close to the no-slip condition.

As for the initial conditions, the important factor is the volume fraction of phases considered in the model. Since there are two physical phases present in this model, i.e. solid and gas, the ratio between these phases must be initially set according to the reality of the model by using Eq. (9).

The current model is a two-phase model of the gas-solid movements in the upper part of the blast furnace. However, to further simplify the model in this study, the velocity of the gas at the tuyers level is assumed to be zero. However, a back flow for gas is permitted at the solid outlet so that the material balance would be fulfilled. Note that even though the gas velocity is set to zero; the equations have been solved for the volume fraction to calculate the voidage of the bed in the furnace.

2.5. Geometry and Mesh Structure

In this study, the simulations have been done for an industrial blast furnace. The geometry has been made based on the blue print of M4 blast furnace of SSAB, located in Oxelösund, Sweden. Figure 1 shows the geometry and size of the furnace. The calculation has been done on a quarter of the furnace. The diameter of the belly is equal to 10.6 m and its height is equal to 20.0 m. As seen in Fig. 1 the constructed domain is only from the tuyeres level up to the throat, since the focus of this study is on the solid flow profile in the upper part of the furnace. The coordinate origin is located at the bottom-left corner of the domain where the symmetric faces intersect as illustrated in Fig. 1. The cut-cell method was used to construct the mesh. Figure 2(a) shows the mesh structure for the furnace geometry. The simulation mesh can be seen in Fig. 2(b). The grid dimension was [53 × 53 × 100], respectively.

Fig. 1.

The size and shape of the SSAB M4 blast furnace.

Fig. 2.

a) Constructed mesh for the furnace geometry. b) Constructed mesh for the domain of the calculation.

In the current model, commercial CFD code PHOENICS was used to solve a two-phase, three-dimensional flow. A staggered grid was employed with the differencing scheme SMART⁵⁰⁾ to solve the partial differential equations by converting it to algebraic finite volume equations. The equations are solved iteratively using the IPSA algorithms.⁵¹⁾

3. Results

The results are presented in two sections. In the first section, the solid phase velocity profile for different heights has been plotted to show the behavior of solid particles while descending through the furnace geometry. In the second section, a comparison has been made between three cases: one case without the extra term in the Navier-Stokes equations represented by Eqs. (14), (15), (16), (17), (18), (19), (20), (21), the other two with the extra term in the Navier-Stokes equation included and the compaction modulus set to 150.0 and 230.0, respectively. Figure 3 shows the position and distance of each grid-slab from the furnace throat.

Fig. 3.

The chosen slabs and their distance from the top.

3.1. Velocity Profile at Different Heights

Figures 4(1)–4(4) show velocity profile for slabs 96, 80, 65 and 45 located in the upper part of the furnace. For these slabs, it can be seen that the maximum of velocity appears along the center axis. Moreover, the magnitude of velocity reduces downwards from slab 96 to 45. The magnitude of the velocity at the center axis for slab 96 is around 0.87 mm/s, while this magnitude is reduced to 0.41 mm/s in slab 45 (11.6 m below the top). However, Fig. 4(5) shows the velocity profile at the belly of the furnace. It can be seen that the velocity profile around the center axis changes its shape and flattens out. It can also be seen that the velocity profile close to the furnace wall changes its shape in comparison to the above slabs. Figure 4(6) shows the velocity profile for solid phase at slab 22. This slab is located 1.6 m above the deadman. This figure shows that the solid phase profile changes its shape and now the maximum velocity magnitude appears at the sides, above where the outlet is located. The velocity profile for slab 19, located 16.7 m below the top in Fig. 4(7), shows the same behavior. Moreover, it is shown that the magnitude of the velocity at the center axis is even lower than the wall-neighboring cells. It should also be mentioned that the velocity magnitude increases at slab 22 in comparison to the higher located slab 30 at the belly of the furnace. This increase in magnitude of the velocity can be seen at all the slabs located below the belly. Slab number 10 is located below the top of the deadman. Therefore, the velocity profiles are divided into two separate sections. At this slab, the magnitude of the velocity for the cells close to the deadman boundary is higher than the wall-neighboring cells. However, this behavior is reversed at the slab number 3, which is located 0.4 m above the outlet. At this slab, the magnitude of the velocity for wall-neighboring cells are higher than for those cells located close to the deadman boundary.

Fig. 4.

(1) Velocity profile for slab 96: (8.6, 8.72)×10^–4 [m/s]. (2) Velocity profile for slab 80: (5.75, 6.05)×10^–4 [m/s]. (3) Velocity profile for slab 65: (4.48, 4.68)×10^–4 [m/s]. (4) Velocity profile for slab 45: (4.03, 4.13)×10^–4 [m/s]. (5) Velocity profile for slab 30: (4.34, 4.42)×10^–4 [m/s]. (6) Velocity profile for slab 22: (4.65, 4.80)×10^–4 [m/s]. (7) Velocity profile for slab 19: (4.8, 5.0)×10^–4 [m/s]. (8) Velocity profile for slab 10: (5.65, 5.90)×10^–4 [m/s]. (9) Velocity profile for slab 3: (6.74, 6.86)×10^–4 [m/s].

3.2. Comparison of Velocity Profiles

As mentioned above, the results from three different cases have been studied: one case without the extra term in Navier-Stocks equations represented by Eqs. (14), (15), (16), (17), (18), (19), (20), (21) and two cases with these terms included and with compaction modulus c equal to 150.0 and 230.0 in Eq. (21), respectively. The results are shown in Figs. 5(1)–5(9). Figure 5(1) shows the velocity profiles for these three cases at slab number 96. The dashed line represents the case without the extra term and the solid and stared lines shows the cases for c equal to 150.0 and 230.0, respectively. Figure 5(1) shows that the case without forces has the highest velocity magnitude and the case with the highest value of c (equal to 230.0) has the lowest magnitude of the velocity. Also, noted is that the magnitude of the difference between the first and the second case is less than the differences in the magnitude of the velocity from the second to the third. The same behavior can be seen for all the other slabs in the lower level. The same comparison of the magnitude of the velocity for Figs. 5(5)–5(6) shows that the case with c equal to 230.0 has a deeper curvature than the case with c equal to 150.0; while for slab 30 in Fig. 5(5) the velocity profiles for the case without forces almost represent a flat line. Moreover, the velocity profile for c equal to 230.0 in Fig. 5(6) shows a stronger effect regarding to the boundary of the deadman. The same behavior can be seen for slabs 19, 10 and 3 too.

Fig. 5.

(1) Solid phase velocity for slab 96 wrt P_e: (8.2, 9.0)×10^–4. (2) Solid phase velocity for slab 80 wrt P_e: (4.0, 7.5)×10^–4. (3) Solid phase velocity for slab 65 wrt P_e: (2.0, 6.5)×10^–4. (4) Solid phase velocity for slab 45 wrt P_e: (1.5, 6.0)×10^–4. (5) Solid phase velocity for slab 30 wrt P_e: (2.5, 6.0)×10^–4. (6) Solid phase velocity for slab 22 wrt P_e: (3.0, 6.5)×10^–4. (7) Solid phase velocity for slab 19 wrt P_e: (3.0, 6.5)×10^–4. (8) Solid phase velocity for slab 10 wrt P_e: (4.6, 6.6)×10^–4. (9) Solid phase velocity for slab 3 wrt P_e: (6.4, 7.0)×10^–4.

4. Discussion

As explained in the previous section, the maximum velocity for slab 96, 80, 65 and 45, which is shown in Figs. 4(1)–4(4), appears at the center axis of the domain. This shows that the solid phase up to the belly of the furnace behave more like a pipe flow. Also, noted is that the velocity magnitude is reduced in this region with lower heights (from slab 96 towards 45). Moreover, based on the simulated mass flow rate the inlet velocity at the top of the furnace is set as a constant value normal to the throat face. Therefore, an increase in the magnitude of the area of any horizontal grid-slab, reduces the magnitude of the velocity since the velocity has an inverse relation to the area of the cross section in calculation of the mass flow. However, Fig. 4(5) shows that for the slab number 30 the velocity at the center axis is flattened out and does not have a sharp edge at this axis. This shows the effect of the deadman boundary on the velocity profile. In other words, the solid particles are reducing their velocity at the center axis due to the fact that there is a blockage area in the middle of the domain. In Fig. 4(6), this can also easily be seen on the profile of the velocity for slab 22. It is seen that the maximum of the velocity is shifted from the center line towards the side above the location of the outlet. Figure 4(7) shows a similar behavior with a higher magnitude of effect for slab 19. In fact, the center axis velocity at this slab is even lower than the magnitude of the velocity of the solid phase close to the furnace wall due to the location of the deadman at the middle-bottom of the furnace. Figures 4(8)–4(9) show the velocity profile of the solid phase at slabs 10 and 3. As seen, the velocity profiles for these two slabs are divided into two different disconnected parts. This is of course due to the size and height of the deadman. Figure 4(8) shows that at slab number 10, the velocity of the solid phase close to the boundary of the deadman is higher in the magnitude than it is in the wall neighboring cells at the same slab. This can be due to the shape of the domain at this slab, which shows that the deadman boundary has a lower angle than the furnace wall. On the other hand, Fig. 4(9) shows the opposite effect. The solid phase has a higher velocity at the wall neighboring cells than the deadman boundary cells at slab number 3. This can be due to the location of the slab. As it can be seen in Fig. 3 slab number 3 is located at the tuyeres level of the furnace geometry. At the outer face it has a cylindrical shape while at the deadman boundary face it follows the slope of the deadman. Therefore, the changes in the magnitude of the velocity for wall neighboring and the deadman boundary cells can be due to the geometry of the domain at this level as well.

As mentioned in the previous section, Figs. 5(1)–5(9) show the velocity profile for solid phase in three different cases: without forces, with forces for compaction modulus equal to 150.0 and the last one with forces for compaction modulus 230.0. Figure 5(1) shows the velocity profiles for these three cases. The magnitude of the velocity is higher for the first case (without forces) than for the case with the compaction modulus equal to 150.0. This means that the extra terms in the Navier-Stocks equation defined by Eqs. (14), (15), (16), (17), (18), (19), (20), (21) behave as a breaking force. Of course, due to the shape and characteristic of beds in a blast furnace, particles are in contact with their neighboring particles which prevent particles from sliding through. In other words, it can be said that this contact force behaves like a friction force and counteracts the vertical velocity. This behavior can be seen in Fig. 5(1) which shows that Eqs. (14), (15), (16), (17), (18), (19), (20), (21) simulate the particle-particle contact forces. In addition, Fig. 5(1) also shows the effect of increase in the magnitude of compaction modulus on the velocity profile of the solid phase. In fact, an increased magnitude of the compaction modulus reduces the velocity profile magnitude. This means that the interaction forces between particles increases when the value of the compaction modulus increases. It should also be mentioned that the case with a larger compaction modulus has a deeper curvature than the other two cases. The case with compaction modulus equal to 230.0 has a deeper curvature than the case with a compaction modulus equal to 150.0. Moreover, the same trend can be seen between the case with the compaction modulus of 150.0 and the case with no forces. This means that the difference between the minimum and maximum of the velocity profile for the case with compaction modulus 230.0 is higher than the other cases. Figures 5(2)–5(9) also show this behavior. In Fig. 5(5) which shows the velocity profile at slab number 30, a stronger effect of this curvature behavior can be seen. The velocity profile for the case with no forces looks almost like a flat line, while for the case with a compaction modulus of 230.0 some curvature can be identified. This is even more clear for slabs 22 and 19. The same trend with larger effect can also be seen for slabs 10 and 3. It is shown that an increase in the magnitude of the compaction modulus c leads to the solid descending velocity decreases. For all three cases, the solid mass flow rate was kept constant. This means that the solid fraction ε_s increases with an increase in the magnitude of the compaction modulus c. Moreover, the solid fraction is given by solving Eq. (2), the solid phase conservation equation. Thus, seen as the core of this newly proposed method, the effect of the magnitude of the compaction modulus c propagates to the solid effective pressure P_e, which effects the rate-independent stress tensor τ_ri, in turn effecting the solid velocity profile (through Eq. (3)), which is coupled to the solid fraction. An increase in the volume fraction of the solid phase ε_s will correspond to a decrease in the magnitude of the solid velocity and thus an increase in the resident time of the particles inside a blast furnace.

5. Conclusion

In this study, the blast furnace process was modeled based on the concept of continuum mechanics. Therefore, it was assumed that the existing phases in the model behave as a continuum media. Hence, the Navier-Stokes differential equation was used to describe the system. The Navier-Stokes equation was solved for two-phase solid-gas system in a three-dimensional domain. The geometry was based on one of the SSAB’s operational furnace located in Oxelösund, Sweden. Since, the solid particles are not exactly a continuum media, extra terms should be added to the Navier-Stokes equation. Therefore, the concept of the solid surface stress was added to the original differential equation to describe the interactions between particles due to the packing, including normal and shear forces. The results showed that the additional extra terms in the Navier-Stokes equation behave as a breaking force in the continuum model to simulate the particle-particle interaction due to the collision and packing. The profile of the deadman has been assumed to have a specified shape and size. It was shown that the particles have a ‘V’ shape profile at the top of the furnace and they follow the same profile from the throat up to the belly. However, they have been affected by the characteristic of the deadman after the belly and their velocity profile change the shape from ‘V’ to ‘W’. It was also shown that an increase in the compaction modulus increases the effective pressure. An increase in the effective pressure increases the magnitude of the surface stress tensor, which in turn reduces the velocities magnitude and affects the velocity profiles accordingly.

Nomenclature

c: Compaction modulus

C_d: Drag coefficient [–]

d_p: Particle diameter [m]

D: Stretching tensor [ ]

$\overrightarrow{\text{g}}$ : Gravity [m/s²]

G₀: Normalizing units factor [–]

I: Identity matrix [–]

p: Pressure [Pa]

P_e: Effective pressure [Pa]

Re_p: Particle Reynolds number [–]

S: Source term

$\overrightarrow{\text{u}}$ : Velocity vector [m/s]

Greek letter

α₁: Cohesion [–]

α₂: Friction coefficient [–]

β: Gas-particle momentum transfer coefficient [kg/m³.s]

Γ: Effective diffusive transfer coefficient

ε: Volume fraction [–]

${\mathit{\epsilon }}_s^{*}$ : Volume fraction of solid particles (packing density) [–]

η_s: Coefficient of plastic modulus [–]

λ_p: Slip parameter [–]

μ: Viscosity [kg/m.s]

ρ: Density [kg/m³]

τ: Stress tensor

Φ: General dependent variable (momentum: ${\overrightarrow{\mathrm{u}}}_{\mathrm{i}}$ )

ϕ: Internal friction angle [o]

φ_s: Sphericity of the particle [–]

Subscripts

i: Phase index

g: Gas

s: Solid

ri: Rate-independent

rd: Rate-dependent

References

• 1)  International Iron and Steel Institute: Steel Statistical Yearbook, IISI, Brussels, (1989).
• 2)   Y.  Omori: Blast Furnace Phenomena and Modelling, Elsevier, London, (1987).
• 3)   J.-i.  Yagi: ISIJ Int., 33 (1993), 619.
• 4)   J.  Chen,  T.  Akiyama,  H.  Nogami,  J.-i.  Yagi and  H.  Takahashi: ISIJ Int., 33 (1993), 664.
• 5)   P. R.  Austin,  H.  Nogami and  J.-i.  Yagi: ISIJ Int., 37 (1997), 748.
• 6)   P. R.  Austin,  H.  Nogami and  J.-i.  Yagi: ISIJ Int., 37 (1997), 458.
• 7)   P. R.  Austin,  H.  Nogami and  J.-i.  Yagi: ISIJ Int., 38 (1998), 246.
• 8)   S. J.  Zhang,  a. B.  Yu,  P.  Zulli,  B.  Wright and  U.  Tüzün: ISIJ Int., 38 (1998), 1311.
• 9)   K.  Takatani,  T.  Inada and  Y.  Ujisawa: ISIJ Int., 39 (1999), 15.
• 10)   J. A. D.  Castro,  H.  Nogami and  J.-i.  Yagi: ISIJ Int., 40 (2000), 637.
• 11)   J. A. D.  Castro,  H.  Nogami and  J.-i.  Yagi: ISIJ Int., 41 (2001), 18.
• 12)   J. A. D.  Castro,  H.  Nogami and  J.-i.  Yagi: ISIJ Int., 42 (2002), 44.
• 13)   J. A. D.  Castro,  H.  Nogami and  J.-i.  Yagi: ISIJ Int., 42 (2002), 1203.
• 14)   S. J.  Zhang,  A. B.  Yu,  P.  Zulli,  B.  Wright and  P.  Austin: Appl. Math. Model., 26 (2002), 141.
• 15)   H.  Nogami and  J.-i.  Yagi: ISIJ Int., 44 (2004), 1826.
• 16)   S. S.  Mondal,  S. K.  Som and  S. K.  Dash: J. Phys. D: Appl. Phys., 38 (2005), 1301.
• 17)   X.  Dong,  A.  Yu,  J.-i.  Yagi and  P.  Zulli: ISIJ Int., 47 (2007), 1553.
• 18)   X. F.  Dong,  a. B.  Yu,  J. M.  Burgess,  D.  Pinson,  S.  Chew and  P.  Zulli: Ind. Eng. Chem. Res., 48 (2009), 214.
• 19)   G. L.  Qing,  L.  Ma,  X. S.  Zhang,  J. L.  Zhou and  M.  Kuwabara: Ironmaking Steelmaking, 37 (2010), 546.
• 20)   S.  Ueda,  S.  Natsui,  H.  Nogami,  J.-i.  Yagi and  T.  Ariyama: ISIJ Int., 50 (2010), 914.
• 21)   J. A.  de Castro,  A. J.  da Silva,  Y.  Sasaki and  J.-i.  Yagi: ISIJ Int., 51 (2011), 748.
• 22)   H.  Kawai and  H.  Takahashi: ISIJ Int., 44 (2004), 1140.
• 23)   T.  Nouchi,  T.  Sato,  M.  Sato,  K.  Takeda and  T.  Ariyama: ISIJ Int., 45 (2005), 1426.
• 24)   Z.  Zhou,  H.  Zhu,  A.  Yu,  B.  Wright,  D.  Pinson and  P.  Zulli: ISIJ Int., 45 (2005), 1828.
• 25)   H. P.  Zhu,  Z. Y.  Zhou,  a. B.  Yu and  P.  Zulli: Granular Matter., 11 (2009), 269.
• 26)   S.  Natsui,  S.  Ueda,  M.  Oikawa,  Z.  Fan,  J.  Kano,  R.  Inoue and  T.  Ariyama: ISIJ Int., 49 (2009), 1308.
• 27)   S.  Natsui,  S.  Ueda,  Z.  Fan,  N.  Andersson,  J.  Kano,  R.  Inoue and  T.  Ariyama: ISIJ Int., 50 (2010), 207.
• 28)   Z. Y.  Zhou,  H. P.  Zhu,  B.  Wright,  a. B.  Yu and  P.  Zulli: Powder Technol., 208 (2011), 72.
• 29)   Z.  Fan,  S.  Igarashi,  S.  Natsui,  S.  Ueda,  T.  Yang,  R.  Inoue and  T.  Ariyama: ISIJ Int., 50 (2010), 1406.
• 30)   Z.  Fan,  S.  Natsui,  S.  Ueda,  T.  Yang,  J.  Kano,  R.  Inoue and  T.  Ariyama: ISIJ Int., 50 (2010), 946.
• 31)   P. A.  Cundall and  O. D. L.  Strack: Géotechnique, 29 (1979), 47.
• 32)   D.  Kolymbas: Arch. Appl. Mech., 61 (1991), 143.
• 33)   G.  Gudehus: Soils Found., 36 (1996), 1.
• 34)   S. A.  Zaïmi,  T.  Akiyama,  J.-B.  Guillot and  J.-I.  Yagi: ISIJ Int., 40 (2000), 322.
• 35)   S. A.  Zaïmi,  T.  Akiyama,  J.-B.  Guillot and  J.-I.  Yagi: ISIJ Int., 40 (2000), 332.
• 36)   S. a.  Zaïmi,  J.-B.  Guillot and  H.  Biausser: Ironmaking Steelmaking, 30 (2003), 475.
• 37)   R.  Safavi Nick,  A.  Tilliander,  L.  Jonsson and  P.  Jönsson: Steel Res. Int., 84 (2013), 1. DOI: 10.1002/srin.201200283
• 38)   K.  Runesson and  L.  Nilsson: Bulk Solids Handling, 6 (1986), 877.
• 39)   P. C.  Johnson and  R.  Jackson: J. Fluid Mech., 176 (1987), 67.
• 40)   P. C.  Johnson,  P.  Nott and  R.  Jackson: J. Fluid Mech., 210 (1990), 501.
• 41)   K.  Hutter and  K. R.  Rajagopal: Continuum Mech. Therm., 6 (1994), 81.
• 42)   S. B.  Savage: Proc. U.S./Japan Seminar on New Models and Constitutive Relations in the Mechanics of Granular Materials, ed. by J. T. Jenkins and M. Satake, Elsevier, New York, (1983).
• 43)   S.  Ergun: Chem. Eng. Prog., 48 (1952), 89.
• 44)   J.  Richardson and  W.  Zaki: Trans. Inst. Chem. Eng., 32 (1954), 35.
• 45)   H.  Norem: Proc. the Davos Symp., Avalanche Formation, Movements and Effects, IAHS, Oxfordshire, (1987).
• 46)   C.  Orr: Particulate Technology, Macmillan, New York, (1966).
• 47)   M.  Swartling: Master Thesis, in Applied Process Metallurgy, KTH Royal Institute of Technology, Stockholm, (2006).
• 48)   S. M.  Eldighidy,  R. Y.  Chen and  R. A.  Comparin: J. Fluids Eng., 99 (1977), 365.
• 49)   M.  Gad-el-Hak: J. Fluids Eng., 121 (1999), 5.
• 50)   P. H.  Gaskell and  A. K. C.  Lau: Int. J. Numer. Method. Fl., 8 (1988), 617.
• 51)   D. B.  Spalding: Numerical Computation of Multi-phase Fluid Flow and Heat Transfer, in Recent Advances in Numerical Methods in Fluids, ed. by C. Taylor and M. Morgan, Pineridge Press, Swansea, (1980), 139.