Preliminary Investigation on the Capability of eXtended Discrete Element Method for Treating the Dripping Zone of a Blast Furnace

Abstract

The role of molten iron and slag in the dripping zone of a blast furnace is very critical to reach a stable operational condition. The existence of several fluid phases and solid particles in the dripping zone of a blast furnace, makes the newly developed eXteneded Discrete Element Method (XDEM) as an Eulerian-Lagrangian approach, suitable to resolve the dripping zone of a blast furnace. In the proposed model, the fluid phases are treated by Computational Fluid Dynamics (CFD) while the solid particles are solved by Discrete Element Method (DEM). These two methods are coupled via momentum, heat and mass exchanges. The main focus of current study is to investigate the influence of packed properties such as porosity and particle diameters, calculated by the XDEM, on the fluid phases for isothermal. In order to present the capability of the XDEM for this application. The validity of the proposed model is demonstrated by comparing model prediction results with the available experimental data.

1. Introduction

Blast Furnace iron making is one of the largest reactor in the world and still the main method to produce molten raw iron, for steel making and foundry purposes. The blast furnace is a vertical cylinder whose diameter varies by its height.¹⁾ In the Blast furnace, the raw materials such as coke and iron ore are charged from the top alternatively while the hot blast gas is introduced through tuyeres from the bottom sides. The iron ore starts reducing partially and where the temperature reaches to its melting point, iron ore starts melting in the cohesive zone. Then, the produced liquid phases start trickling toward the dripping zone. Two liquid phases, molten iron and slag, with different properties exist in the dripping zone. They flow downwards by gravity force through piles of coke particles while hot gas ascends upward. The dripping zone can be classified as a counter-current and cross-current trickle bed reactor.²⁾

The flow of molten iron and slag plays a critical role in the dripping zone of the blast furnace while transferring mass and energy with the gas phase.^3,4) These phases not only affect the gas phase distribution and composition but also affect the pressure and temperature inside the blast furnace. Therefore, it is necessary to understand the flow behavior of liquid phases in this particular zone, where the pressure and temperature is so high that experimental studies are either inapplicable or very expensive to carry out. In addition, an acceptable flow of hot metal and slag is essential in the lower part of the blast furnace to reach a stable operation in order to improve the production rate.⁵⁾

Many studies have investigated the flow behavior of the liquid phase with the closed condition to the dripping zone of a blast furnace experimentally and numerically. Some early experimental studies on the flow behavior of gas and liquid through a packed bed of particles gave some preliminary principles applicable to the dripping zone of a blast furnace.^6,7,8,9,10) According to these experiments, the distribution of the liquid flow is highly affected by the gas flow insofar as the gas flow tended to displace the liquid phase at the vicinity of its entrance where the dry zone was constructed. They also concluded that the size of dry zone increases by increasing the gas flow rate.^6,9) Other hydrodynamic parameters such as liquid hold up, gas pressure drop and flooding were also studied.⁷⁾ The effect of different particle diameters and radial distribution of the liquid phase was investigated by an experimental study conducted by Eto et al.⁸⁾ while the interaction between all phases was studied by Wang and his colleagues.¹⁰⁾

From the numerical point of view, mainly two methods were used to model dripping zone: Discrete and Continuous methods. In the first method, the discrete phase either liquid or solid particles were treated as discrete entities while in the second method all phases are consider as continuous media. In 1993, Yagi¹¹⁾ reported a review over a continuous flow behavior of four phases of gas, powders, liquid and packed particles to obtain a numerical solution for the whole blast furnace process. Yagi’s¹¹⁾ model was used and modified by Austin et al.^12,13) toward more actual blast furnace operational conditions. Aside from the above mentioned models, which are based on continuous approach, Gupta et al.⁹⁾ proposed a force balance model for the liquid phase, which considers liquid flow as rivulets and drops under the three forces of gas drag, gravity and bed resistance. In Gupta’s method the liquid phase is treated as a discrete phase since there is non-wetting flow condition between liquid and solid particles. Wang et al.^10,14,15) have added stochastic motion related to the packing geometry to Gupta’s model. Their treatment considered only the downward liquid flow without addressing other possible liquid directions and may be justified by the fact that the structure of the packed bed is reasonably isotropic. The force balance method was also used by Chew et al.^3,16) to focus more on the interaction between phases which consider both counter-current and non-counter-current flows. Recently, Kon et al.¹⁷⁾ developed a new numerical model based on Moving Particle semi-implicit (MPS) method combined with Computational Fluid Dynamics (CFD) to predict and understand in furnace phenomena. They applied their model to the dripping zone of the blast furnace and investigated the liquid physical properties in different conditions.

The mathematical treatment of phases in both discrete and continuous approaches seems to have their own strengths and weaknesses, which are both accepted as a promising way of predicting fluid flow behavior in the blast furnace technology. From the theoretical and practical points of view, it might be more reliable to consider the discrete approach for the liquid phase inside the cohesive zone where the liquid phase is generated. But a continuous approach for the dripping zone where the liquid phase has been generated and has a continuous behavior. Whereas, the discrete particle method, mainly based on Discrete Element Method (DEM), offer a more convenient way to track the movement of the solid coke particles in the dripping zone of a blast furnace.¹⁸⁾ In this regard, the combined CFD-DEM approaches are highly appreciated to fulfill this gap. Despite several studies in CFD-DEM modeling of different blast furnace part,^19,20) the dipping zone including liquid phase has not been yet treated with this combined approach. Therefore, in this study all fluid phases such as gas and liquid phases are treated using Computational Fluid Dynamics as a continuous method while the solid particles are treated using eXtended Discrete Element Method (XDEM) as a discrete method. The coupling of the model is performed by introducing momentum transfer between the discrete and continuous phases in order to consider interaction between all phases.

The Extended Discrete Element Method is a multi-physics, numerical simulation method, which is based on classical discrete element method as enhanced by supplementary properties such as the thermodynamic state or stress/strain for each element. This method can be coupled to continuous methods such has CFD and/or Finite Element Analysis (FEA) by covering the interaction between discrete and continuous phases such as fluids or solid structures at the same time. The XDEM has been developed by Peters²¹⁾ and applied to many Multiphysics applications such as combustion,^22,23) heat up processes,²⁴⁾ drying,^25,26,27) pyrolysis^28,29,30,31) as well as reduction^32,33,34,35) and many other applications.^36,37,38) So far, the upper part of the blast furnace known as the shaft, where the reducing gas is the only fluid phase was modeled by XDEM method. This method has been validated not only for the rheology of the gas phase through layered packed bed of iron ores and coke particles in the shaft of a blast furnace but also for the drying and different reduction of iron ores.^33,35,39)

The main target of this paper is to present an Eulerian-Lagrangian mathematical model for the dripping zone of a blast furnace by introducing liquid phases into the XDEM. This model gives more detailed information for the distribution of the liquid phases through packed bed of coke particles at different gas velocities, which cannot be realized by the continuous methods. Almost all continuous methods either use constant porosity distribution^1,2) or benefit the correlation^13,40) for porosity calculation, while the XDEM itself calculates the porosity distribution by considering the position of the particles. Although this model has the capacity to consider several liquid phases, in the current study the validation was done solely for one liquid phase. The results have been validated using experimental study conducted by Szekely and Kajiwara⁶⁾ and the same cases with the molten iron and slag properties were repeated to study the effect of their density and viscosity on different parameters such as liquid and gas velocities as well as liquid hold up.

2. The Mathematical Model and Numerical Solution Based on XDEM

The dripping zone of an iron making blast furnace is a counter-current reactor that involves several fluid and solid phases. The existence of continuous and particulate phases in this zone makes the coupled Eulerian-Lagrangian methods valuable enough to be investigated. In this section, the mathematical model of the XDEM for the dripping zone is described.

2.1. The Mathematical Model for Continuous Phases

On the one hand, tracking the local properties of each phase is very complex and in most cases not possible, especially at the interfaces in multiphase systems.⁴¹⁾ On the other hand, the macroscopic aspect of the multiphase flow is more important to the design and operation of a multiphase system. Therefore, an appropriate averaging approach could yield the mean value of flow and thermal properties. Therefore, the widely-used Eulerian volumetric average⁴¹⁾ was applied to the governing equations, where the volume fraction of each fluid phase, i, is porosity times saturation (ε_i=ϕα_i). In this regard, the conservation of mass can be written as:   

$$\frac{\partial }{\partial t}\left(\varphi {\alpha }_i{\rho }_i\right)+\nabla \cdot \left(\varphi {\alpha }_i{\rho }_i{\stackrel{\rightharpoonup }{v}}_i\right)=\sum _{j=1,j\ne i}^n{\dot{m}}_{ji}$$(1)

The right-hand side of Eq. (1) represents mass transfer per unit volume from all other phases to the i^th phase. The total mass equation can be derived by summation over mass conservation for all phases:   

$$\frac{\partial }{\partial t}\left(\sum _{i=1}^n\varphi {\alpha }_i{\rho }_i\right)+\nabla \cdot \sum _{i=1}^n\varphi {\alpha }_i{\rho }_i{\stackrel{\rightharpoonup }{v}}_i=0$$(2)

With the assumption that all phases share the same pressure field, the pressure equation can be derived using total mass Eq. (2) by substituting densities with pressure using equation of state for each phase. The ideal gas law was considered for the gas phase as ${\rho }_g= \frac{p}{RT}$ and constant speed of sound for the liquid phase as ${\rho }_l= {\rho }_{lo}\left(1+\frac{1}{{\rho }_{l0}{C}^2}(p-p_0)\right)$.

On this subject, the conservation of momentum for each phase can be expressed as follows:   

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\varphi {\alpha }_i{\rho }_i{\stackrel{\rightharpoonup }{v}}_i\right)+\nabla \cdot \left(\varphi {\alpha }_i{\rho }_i{\stackrel{\rightharpoonup }{v}}_i{\stackrel{\rightharpoonup }{v}}_i\right)=\\ -\nabla \cdot \left(\varphi {\alpha }_{i}p\right)+\varphi {\alpha }_i{\rho }_{i}g+\nabla \cdot \left(\varphi {\alpha }_i\overline{{\tau }_i}\right)+\sum _{j=1}^n{F}_i^j\end{array}$$(3)

The last two terms ${⟨\overline{{\tau }_i}⟩}^i$ and $F_i^j$ represent the phase average stress tensor and the momentum exchange between phases respectively and are defined as Eqs. (4) and (5).   

$$\overline{{\tau }_{\text{i}}}={\mu }_i\left[\nabla {\stackrel{\rightharpoonup }{v}}_i+{\left(\nabla {\stackrel{\rightharpoonup }{v}}_i\right)}^T\right]-\frac{2}{3}{\mu }_i\left(\nabla \cdot {\stackrel{\rightharpoonup }{v}}_i\right)\delta$$(4)

  

$$F_i^j=K_{ij}\left({\stackrel{\rightharpoonup }{v}}_j-{\stackrel{\rightharpoonup }{v}}_i\right)$$(5)

In this modelling, drag force is the only force considered between each two pair of phases. In Eq. (5), K_ij represents the drag force coefficient. In this particular application of the XDEM model, which is a one way coupling between the fluid phases and solid particles (since the particles assumed to be stagnant after settling), the very well-known Ergun correlation was used to consider momentum transfer between solid particles and fluid phases and Schiller-Naumann correlation was used to include drag force between the gas and liquid phases.   

$$K_i^j=\frac{150{\left(1-{\epsilon }_j\right)}^2{\mu }_i}{{\epsilon }_i{d}_j^2}+\frac{1.75\left(1-{\epsilon }_j\right){\rho }_i}{d_j}\left|{\stackrel{\rightharpoonup }{v}}_j-{\stackrel{\rightharpoonup }{v}}_i\right|$$(6)

  

$$K_{ij}=\frac{3}{4}{C}_{D,i}\frac{{\epsilon }_i{\rho }_i}{d_i}\left|{\stackrel{\rightharpoonup }{v}}_j-{\stackrel{\rightharpoonup }{v}}_i\right|$$(7)

  

$$C_D= \{\begin{array}{cc}\frac{24}{Re}\left(1+0.15R{e}^{0.687}\right),& Re<1\mathrm{\hspace{0.17em} }000\\ 0.44,& Re\ge 1\mathrm{\hspace{0.17em} }000\end{array}$$(8)

  

$$Re=\frac{d_i {\rho }_i\left|{\stackrel{\rightharpoonup }{v}}_j-{\stackrel{\rightharpoonup }{v}}_i\right|}{{\mu }_i}$$(9)

The porosity and particle diameter must also be known to solve the governing equations for the fluid phases, which are not uniform in the dripping zone of a blast furnace. In order to calculate the drag force due to the solid particles on the fluid phases, the XDEM calculates the volumetric average value of porosity for each CFD cell by knowing the number of particles in each CFD cell based on a kernel-based interpolation procedure proposed by Xiao et al.⁴²⁾ According to this method, the volume of each particle is weighted (η) to all the CFD cells according to its distance to the cell centers. Therefore, the average value of porosity for each CFD cell can be calculated as:   

$${\varphi }_{cell}=1-\frac{1}{V_{cell}}\sum _{p=1}^N{\eta }_{k,cell}{V}_p$$(10)

In order to avoid numerical instability due to the zero porosity, the smallest CFD cell volume must be bigger than the largest particle volume, which is known as unresolved method. The mean particle diameter, D_p, is calculated from the total mean surface area, S, for N particle in the given CFD cell.   

$$\overline{S}=\frac{1}{N}\sum _{p=1}^N{A}_p$$(11)

  

$${\overline{D}}_p=\sqrt{\frac{\overline{S}}{\pi }}$$(12)

2.2. The Mathematical Model for Discrete Entities

The XDEM platform utilizes the traditional DEM in order to tracks the movement of the particles using Newton’s second law for the translational (v_i) and rotational (ω_i) motion of particle i.   

$$m_i\frac{d{v}_i}{dt }=\sum _j{F}_{c,ij}+F_{g,i}$$(13)

  

$$I_i\frac{d{\omega }_i}{dt}= \sum _j\left(M_{c,ij}+M_{ext}\right)$$(14)

Where m_i and I_i are the mass and inertia tensor of particle i, respectively. F_c,ij is the contact force between particle i and j, F_g,i is the gravitational force, M_c,ij is the torque acting on particle i due to the collision with other particles and M_ext any other external torques. A comprehensive and detailed description of the dynamic module of the XDEM is reported by Samiei et al.,^43,44,45) which is beyond the purpose of this study. In this study, linear spring-dashpot model with spherical glass bead properties were used.

2.3. Initial and Boundary Conditions

The geometry of the case, which is based on the experimental study conducted by Szekely and Kajiwara⁶⁾ at room temperature, is shown in Fig. 1. On the upper part of the case, as indicated in Fig. 1, some cells were specified to produce liquid, which is defined as a mass source in the continuity equation of liquid phase. Fixed uniform inlet velocity was specified for the gas phase while no slip boundary condition for the walls and the total pressure for the outlet boundary condition.

Fig. 1.

(a) Experimental setup and dimensions.⁶⁾ (b) Simulation domain and grid. (Online version in color.)

The physical properties of the fluids used in this study are reported in Table 1. By considering the following assumption,

Table 1. Fluid phases properties.

Phase	ρ (Kg/m3)	μ (Kg/ms)
Air	1.18	1.7e−7
Water	1000	1.0e−3
Liquid iron	6600	5.0e−3
Slag	2600	0.1

• all phases are compressible (liquid phases are assumed to be weekly compressible by considering the speed of sound as 1482 m/s for the liquid phase to treat all the fluid phases in the same way)

• isothermal system

• laminar flow

• no inter-phase mass transfer

• Liquid and gas phase share the same pressure (a mixture pressure field based on the continuity equations is used.⁴⁷⁾)

the governing equations were solved based on the algorithm discussed in the next section.

2.4. Solution Algorithm

The initial position and packed structure of the particles were determined using the DEM module of the XDEM discussed in section 2.2. In the coupling module of the XDEM, which is responsible for transfer and exchange of data between the mesh-based Eulerian fields of the continuous phase and the particle-based-Lagrangian quantities of the particulate phase, the porosity field and mean particle diameter is calculated. This data is accessible to the CFD solver to calculate the momentum transfer term. Finally, with the assumptions mentioned in section 2.3, the momentum equation Eq. (3) and mass conservation equation Eq. (1) are solved by a segregated projection algorithm based on the pressure implicit of operators (PISO) algorithm exist in OpenFOAM.^47,48)

The temporal derivative was discretized using first order implicit Euler scheme while the spatial derivatives were discretized by either Gauss linear or limitedlinear (first order), presented in standard OpenFOAM source code. Some preliminary numerical cases were set up to investigate mesh independence of the solver. The optimal uniform grid size of 6 mm in each direction was used, where no significant changes were observed by refining the mesh. The residuals of less than 10⁻⁴ were accepted by using adjustable time step to satisfy Courant number of less than 0.5.

3. Results and Discussion

In the first step, the code was validated using the experimental study conducted by Szekely and Kajiwara⁶⁾ on the maldistributed counter-current flow of gases and liquids. The main motivation of their work was its relevance to the flow phenomena in the lower part of the iron blast furnace. In the next step, the model was applied to liquid iron and slag to compare the results and to investigate the effect of the liquid properties such as density and viscosity on the calculated parameters.

3.1. Code Validation

Szekely and Kajiwara⁶⁾ examined the maldistribution pattern of the gas and liquid streams by designing the counter-current flow of air and water through a packed bed of glass spheres. The experimental set up and dimensions are shown in Fig. 1(a). The liquid phase is introduced into the domain through a distributer located on the top and the gas phase with a certain flow rate, was blasted through the tuyere. In Fig. 1(b), the simulation domain with CFD grid are shown, 60×150×8 cells in x, y and z directions, where some cells on the top were specified to produce liquid with the same mass flow rate as indicated in the experimental study in terms of cubic meter per minute.

The importance of porosity distribution or in another word bed permeability on the operation of a blast furnace has been previously discussed.^2,18,49) The permeability of the gas phase and liquid phases hold up, which govern the flow of gas and liquid phases, are key factors to smooth and stable blast furnace operation. These parameters are highly sensitive to the porosity distribution. In order to calculate the porosity field, the column was first filled with particles of two different sizes (3 and 6 mm). The position of the particles as shown in Fig. 2(a), was introduced to the CFD cells, then the mean particle diameter and the porosity field were calculated by XDEM as shown in Figs. 2(b), 2(c), respectively. More detailed information on the calculated porosity field and the mean particle diameter can be found in Fig. 3. As shown in this figure, the CFD cells including both size of particles will receive an average value for particle diameter. Higher porosity was calculated for d_p=6 mm, which is more obvious in Fig. 4. In this figure, the porosity distribution over two lines in x direction at z=0.025 m, one in the height of y=0.3 m, where particle diameters are d_p=6 mm and the other one at the height of y=0.7, where the particle diameters are d_p=3 mm are shown to better observe the wall effect on the porosity distribution and in consequence on gas and liquid velocities as well as liquid saturation. Higher porosity is observed at the vicinity of the walls which can effect calculated parameters and predict a non-uniform description of calculated parameters. This XDEM advantage eliminates the requirement of empirical correlations and avoid the errors due to the constant porosity assumption. The average porosity calculated by XDEM in the balk of particles is 0.362 for d_p=3 mm and 0.381 for d_p=6 mm, while in the experimental study the porosity of 0.33, 0.34 were obtained. The deviation of approximately 10% was calculated, which could be due to the small spaces created during simulation between the particles. In the experimental studies, normally the beds are shacked to fill the empty gaps with particles leading to lower porosities.

Fig. 2.

(a) Simulation set up and particle positions. (b) Calculated mean particle diameter by the XDEM. (c) Calculated porosity distribution by the XDEM. (Online version in color.)

Fig. 3.

(a) Particle positions for a part of a domain indicated in Fig. 3, a. (b) Calculated mean particle diameter by the XDEM for the specified zone. (c) Calculated porosity distribution by the XDEM for the specified zone. (Online version in color.)

Fig. 4.

The magnitude of porosity, liquid saturation, gas velocity and liquid velocity over two horizontal line in the central x-y plane, where the dash line is at height of y=0.3 m and solid line is at height of y=0.7 m for $Q_l^0=650 \frac{c{m}^3}{min}$.

The simulation cases were set up for two different liquid flow rates ($Q_l^0=650 \frac{c{m}^3}{min}$ and $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$) as was carried out by the experimental study. In the first step, the variables were computed for both liquid flow rates without introducing the gas phase as shown in Fig. 5. Figure 5(a) shows the iso-surfaces for liquid saturation while Fig. 5(b) displays the vectors of the liquid phase. Through the results shown in Fig. 5, higher liquid hold up and liquid velocities were observed by increasing the liquid mass flow rate in the figure. The reason for this may be that as the liquid flow rate increases, higher velocities are obtained due to the higher gravity force, which leads to higher momentum transfer and in consequence, more hold up. In addition, the liquid hold up in the zone with d_p=3 mm is higher than in the zone with d_p=6 mm while the liquid velocity is smaller, which is due to the lower solid porosity (lower permeability) in the zone with smaller particle diameters. Therefore, the total resistance to the liquid flow increases as the particle size decreases.⁵⁰⁾ This figure also shows that the liquid distribution is not fully homogeneous due to the porosity distribution inside the bed, which is not uniform as well. For example, higher velocities are predicted near to the wall zone, where porosity is also higher as was also discussed in Fig. 4. This capability makes the current method close to the reality compared to other studies with the assumption of constant porosity distribution.^1,2) The liquid hold up reported in this study is the resultant of the particles drag force and gravitational force without taking into account the capillarity effects. The importance of capillary pressure on the liquid hold up in packed is non-negligible. But its effect on liquid hold up become less dominant in the non-wetting condition where the contact angle is close to 180⁰.¹⁷⁾ Therefore, the reported results for liquid hold up may be higher by considering capillarity pressure.

Fig. 5.

(a) Iso-surfaces of the liquid saturation. (b) Velocity vector of the liquid phase. (Online version in color.)

In the next step, different gas flow rates were introduced to investigate the effect of gas phase on the deviation of liquid flow. The measured liquid flow rate based on the experimental data was compared with the XDEM predictions for the liquid flow rate of $Q_l^0=650 \frac{c{m}^3}{min}$ in Fig. 6 and $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$ in Fig. 7. These figures plot the water flow rate percentage as a function of dimensionless distance in x direction at the outlet for different gas flow rates. The water flow rate percentage is defined as the percentage of deviation from the zero gas inlet condition. Therefore, zero percentage refer to zero deviation from the case in the absence of the gas flow, while the negative value correspond to decrease in local liquid flow rate and positive value increase in the local liquid flow rate at the outlet. Very good agreement between the XDEM prediction results and the experimental data was observed for both liquid flow rates. In Figs. 8 and 9, the calculated outlet liquid mass flow rate based on the deviation obtained in the experimental study against the XDEM prediction results are pointed in order to investigate the error of the proposed method. Although, the absolute mean error of the prediction results seems to slightly increases by enlarging the gas flow rate, the total absolute mean error of 3.27 and 4.68 for $Q_l^0=650 \frac{c{m}^3}{min}$ and $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$ respectively, represents the robustness of the XDEM solver for these cases. Despite higher fluctuations for $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$, the trend of the data seems to be reasonable, which represent more negative value at the interance of the gas phase for the higher gas inlet velocities.

Fig. 6.

Comparison between the XDEM calculated results (dash line) for the deviation of water flow rate at different gas inlet flow rate against experimental results (dots) for $Q_l^0=650 \frac{c{m}^3}{min}$.

Fig. 7.

Comparison between the XDEM calculated results (dash line) for the deviation of water flow rate at different gas inlet flow rate against experimental results (dots) for $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$.

Fig. 8.

Comparison between the XDEM calculated results (orange dots) for the outlet water flow rate at different gas inlet flow rate against experimental results (black squares) for $Q_l^0=650 \frac{c{m}^3}{min}$ with the corresponding absolute mean error%. (Online version in color.)

Fig. 9.

Comparison between the XDEM calculated results (orange dots) for the outlet water flow rate at different gas inlet flow rate against experimental results (black squares) for $Q_l^0=650 \frac{c{m}^3}{min}$ with the corresponding absolute mean error%. (Online version in color.)

In Fig. 10, the iso-surfaces of liquid saturation for different inlet gas flow rates are shown while the liquid flow rate is $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$. The results show that the increase in gas flow rate leads to higher liquid hold up, which is due to the higher interaction area between phases which gives rise to higher resistance for the flow. The liquid saturation is nearly zero at the entrance of the gas phase called the dry zone, meaning that liquid is pushed away from the raceway zone due to the high interaction between the gas and liquid phase.

Fig. 10.

Iso-surfaces of the liquid phase saturation when $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$ for different inlet gas flow rates. (Online version in color.)

The influence of the inlet gas flow rate on the size of the dry zone as well as the liquid phase vectors is clearer in Fig. 11. The velocity vectors of the liquid phase are boosting by increasing the gas inlet flow rate due to the higher interaction between these two fluid phases. It can also be concluded from this figure that the area of the dry zone is proportional to the inlet gas flow rate. This fact is also demonstrated in Fig. 12, which illustrates the liquid phase stream lines for inlet liquid flow rate of $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$. In this figure, the deviation of the liquid phase due to the gas flow rate is shown. As was mentioned before, the liquid phase flows downward except at the raceway zone which has radial variation as well. This radial variation is also proportional to the gas flow rate. The vectors and stream lines for the gas phase are shown in Figs. 13 and 14 to compare the magnitude of gas velocities in different conditions. The gas velocity reduces rapidly due to the porosity effects. Higher velocities were observed at the walls where lower liquid volume fraction exist since there is higher permeability for the gas phase to flow.

Fig. 11.

Liquid phase velocity vectors in a part of a domain indicated in Fig. 3, a when $Q_l^0=650 \frac{c{m}^3}{min}$ for different inlet gas flow rates. (Online version in color.)

Fig. 12.

Liquid phase streamlines when $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$ for different inlet gas flow rates. (Online version in color.)

Fig. 13.

Gas phase velocity vectors when $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$ for di_erent inlet gas flow rates. (Online version in color.)

Fig. 14.

Gas phase streamlines when $Q_l^0=1\mathrm{\hspace{0.17em} }500\frac{c{m}^3}{min}$ for different inlet gas flow rates. (Online version in color.)

3.2. Slag and Liquid Iron

Although the physical and thermophysical properties of water are different from the molten iron and slag, this preliminary study provides basic information for the rheological behavior of liquid phases and on the other hand, the ability of the XDEM method for the dripping zone of a blast furnace. Therefore, the same cases with the isothermal properties of liquid iron and slag listed in Table 1, have been examined to compare the real liquid phases in the dripping zone.

In order to examine liquid hold up, the same amount of liquid ($Q_l^0=650 \frac{c{m}^3}{min}$) was introduced for both liquid phases (liquid iron and slag), shown in Fig. 15, in the absence of the gas flow. The slag has a very higher viscosity than the liquid iron, which leads to higher resistance due to the solid particles, therefore, higher liquid hold up. On the other side liquid iron density is around three times higher than the slag density which accelerates the downward flow of the liquid phase due to gravitational force, leading to lower liquid hold up for this phase and higher velocities (see Fig. 16). As was explained, higher velocities are observed for liquid iron according to Fig. 16. Although constant viscosity is assumed for the liquid phases, the viscosity of the slag in not constant and changes with its composition and temperature while moving downward whereas the viscosity of molten iron is not changing significantly. The correct calculation of slag viscosity is important since it considerably controls the blast furnace operational conditions such as liquid hold up and dry zone size. This matter should be consider in non-isothermal studied, which is not in the scopes of the current contribution.

Fig. 15.

Iso-surfaces of the liquid saturation for Slag (left) and liquid iron (right) without inlet gas flow rate. (Online version in color.)

Fig. 16.

Velocity vectors of Slag (left) and liquid iron (right) without inlet gas flow rate. (Online version in color.)

In Figs. 17 and 18 some snapshots for different times are shown to compare the speed of liquid iron and slag movement through the packed bed of particles. The results show that liquid iron flows faster than slag due to its higher density and lower viscosity. It is worth to mention that the coexistence of both liquid phases in the dripping zone, leads to an interaction between them. The momentum transfer between liquid iron and slag can accelerate or decelerate each other’s velocities while moving through packed bed of particles, which is not discussed in this study. Individual investigation on the resistant time of each liquid phase predicts that the same amount of slag requires approximately two times more than liquid iron to exit from the bottom.

Fig. 17.

Iso-surfaces and the velocity vectors of liquid iron at different snapshots. (Online version in color.)

Fig. 18.

Iso-surfaces and the velocity vectors of Slag at different snapshots. (Online version in color.)

In order to examine and compare the effect of the gas phase on the liquid iron and slag, the gas inlet flow rate of ($Q_g^0=0.08\frac{N{m}^3}{min}$) was introduced through the gas inlet patch. The liquid mass flow rate percentage as a function of dimensionless distance show in Fig. 19, represent that liquid iron is more affected by the gas phase. Despite the high liquid saturation of slag, which increases the interaction area, it is pushed less away from the raceway zone due to its higher resistance to the gas phase. The liquid saturation of liquid iron and slag at the raceway zone are shown in Fig. 20. Since the liquid iron is pushed more toward the central zone due to the gas phase, bigger dry zone is also predicted for this phase. In addition, in Fig. 20, a liquid stagnant zone right above the race way zone is created which is larger for the liquid iron for the same reason.

Fig. 19.

Liquid flow rate percentage as a function of dimensionless distance at $Q_g^0=0.08\frac{N{m}^3}{min}$ for Liquid iron (solid line) and slag (dash line).

Fig. 20.

Liquid saturation at the entrance of the gas phase with $Q_g^0=0.08\frac{N{m}^3}{min}$ for slag (left) and iron (right). (Online version in color.)

The hydrodynamic parameters of fluid phases are mainly related to the porosity and particle size in a blast furnace, which vary from the central zone toward the wall for constant particle zones. The XDEM method provides a promising method to consider these effects by solving a one way Eulerian-Lagrangian equations. However, it should be noted that the proposed model should be extended to non-isothermal case to consider the effect of temperature on the rheology of the fluid phases. This aspect will be examined in future work.

4. Conclusions

A preliminary investigation on the capability of the eXtended Discrete Element in conjunction with OpenFOAM was carried out in order to predict flow behavior of fluid phases in the dripping zone of a blast furnace iron making industry. The presented model is a multi-scale Eulerian-Lagrangian method, which solves the momentum and mass conservation for each phase using CFD while the solid particles are treated by XDEM.

The results have been validated using experimental data. The results confirmed the effect of packing structure on liquid saturation and velocities. The influence of non-uniform porosity distribution on these parameters shows that the liquid and gas velocities are not uniform. For instance, the gas and liquid phase velocities are higher close to the wall, where the porosity is higher. The drag force due to the gas phase on the liquid phase was also investigated, which shows that the gas phase tends to replace the liquid phase at the entrance of the gas phase. This behavior leads to production of a dry zone, whose size is proportional to the gas inlet flow rate. In addition, the cross flow of the liquid phase was also observed aside from its downward flow at the vicinity of the tuyere.

The liquid phase with constant liquid iron and slag properties was also studied to discuss the importance of viscosity and density on the flow pattern of fluid phases. This analysis shows that there is a strong relation between the fluid properties and its flow behavior through packed bed of solid particles, which is an effective parameter to control the stability of the blast furnace operational condition.

Through this study, it has been found that the XDEM is a reliable method for resolving the flow behavior of fluid phases in the dripping zone of a blast furnace. This Eulerian-Lagrangian approach could provide an excellent opportunity to better understand and control the in furnace process, where measurements are difficult to handle.

Acknowledgment

The authors are thankful to the Luxembourg National Research Fund (FNR) for the financial support of this study and grateful to the High Performance Computers (HPC) team of university of Luxembourg for their support.

Nomenclature

c: Speed of sound in liquid (m.s⁻¹)

C_D: Drag coefficient

d: Particle diameter (m)

D: Mean particle diameter (m)

A: Mean surface area (m²)

A: Surface area (m²)

I: Momenum of Inertial (Kg.m²)

F: Force (J)

g: Gravitational acceleration (m.s⁻²)

K: Momentum transfer (J.s.m⁻¹)

M: Torque (N.m)

$\dot{m}$: Mass flow rate (Kg.m⁻³.s⁻¹)

n: Number of fluid phases

N: Number of particles

p: Pressure (Pa)

Q⁰: Volumetric flow rate (m³.min⁻¹)

R: Universal gas constant (J.mol⁻¹.K⁻¹)

Re: Reynolds number

t: Time (s)

S: Total surface area (m²)

T: Temperature (K)

V: Volume (m³)

Greek symbols

∂: Differential operator

∇: Nabla operator

Δ: Difference

ρ: Density (Kg.m⁻³)

ψ: Variable

ν: Velocity (m.s⁻¹)

ε: Volume fraction

τ: Stress strain tensor (Kg.m⁻¹.s⁻²)

μ: Dynamic viscosity (Kg.m⁻¹.s⁻¹)

ϕ: Porosity

α: Volume saturation

η: Weight of particle

δ: Unit tensor

ω: Angular velocity (Rad.s⁻¹)

Subscripts and Superscripts

f: fluid

g: gas

i,j: phase

l: liquid

p: particle

cell: computational fluid dynamics cell

T: transpose

0: reference

References

• 1)   H.  Nogami,  M.  Chu and  J.  Yagi: Comput. Chem. Eng., 29 (2005), Issue 11–12, 2438.
• 2)   X. F.  Dong,  A. B.  Yu,  J. M.  Burgess,  D.  Pinson,  S.  Chew and  P.  Zulli: Ind. Eng. Chem. Res., 48 (2009), 214.
• 3)   S. J.  Chew,  P.  Zulli and  A.  Yu: ISIJ Int., 41 (2001), No. 10, 1112.
• 4)   S.  Ghosh and  N. N.  Viswanathan and  N. B.  Ballal: Steel Res. Int., 87 (2017), 9999.
• 5)  H. J. Jin and S. Choi: Ironmaking Steelmaking, 37 (2010), 89.
• 6)   J.  Szekely and  Y.  Kajiwara: Metall. Trans. B, 10 (1979), 447.
• 7)   T.  Fukutake and  V.  Rajakumar: Trans. Iron Steel Inst. Jpn., 22 (1982), 355.
• 8)   Y.  Eto,  K.  Takeda,  S.  Miyagawa,  S.  Taguchi and  I.  Hiroshi: ISIJ Int., 33 (1993), 681.
• 9)   G. S.  Gupta,  J. D.  Litster,  V. R.  Rudolph,  T. E.  White and  A.  Domanti: ISIJ Int., 36 (1996), 32.
• 10)   G. X.  Wang,  D. Y.  Liu,  J. D.  Litster,  A. B.  Yu,  S. J.  Chew and  P.  Zulli: Chem. Eng. Sci., 52 (1997), 4013.
• 11)   J.  Yagi: ISIJ Int., 33 (1993), 619.
• 12)   P. R.  Austin,  H.  Nogami and  J.  Yagi: ISIJ Int., 37 (1997), 458.
• 13)   P. R.  Austin,  H.  Nogami and  J.  Yagi: ISIJ Int., 38 (1998), 246.
• 14)   G. X.  Wang,  S. J.  Chew,  A. B.  Yu and  P.  Zulli: ISIJ Int., 37 (1997), 1153.
• 15)   G. X.  Wang,  S. J.  Chew,  A. B.  Yu and  P.  Zulli: Metall. Mater. Trans. B, 28 (1997), 333.
• 16)   S. J.  Chew,  P.  Zulli and  A. B.  Yu: ISIJ Int., 41 (2001), 1122.
• 17)   T.  Kon,  S.  Natsui,  S.  Ueda,  R.  Inoue and  T.  Ariyama: ISIJ Int., 53 (2013), 590.
• 18)   X.  Dong,  A.  Yu,  J.  Yagi and  P.  Zulli: ISIJ Int., 47 (2007), 1023.
• 19)   H.  Kurosawa,  S.  Matsuhashi,  S.  Natsui,  T.  Kon,  S.  Ueda,  R.  Inoue and  T.  Ariyama: ISIJ Int., 52 (2012), 1010.
• 20)   A.  Adema,  Y.  Yang and  R.  Boom: 7th Int. Conf. on CFD in the Minerals and Process Industries, CSIRO, Melbourne, (2009), 1.
• 21)   B.  Peters: Thermal Conversion of Solid Fuels, WIT Press, Ashurst, United Kingdom, (2003), 1.
• 22)   B.  Peters and  J.  Smuła-Ostaszewska: Energy Fuels, 24 (2009), 945.
• 23)   A. H.  Mahmoudi,  X.  Besseron,  F.  Hoffmann,  M.  Markovic and  B.  Peters: Chem. Eng. Sci., 142 (2016), 32.
• 24)   B.  Peters and  A.  Dziugys: Mechanics, 20 (2014), 24.
• 25)   B.  Peters,  X.  Besseron,  A.  Estupinan,  F.  Hoffmann,  M.  Michael and  A. H.  Mahmoudi: IFRF Combust. J., 14 (2014), 1.
• 26)   A. H.  Mahmoudi,  F.  Hoffmann and  B.  Peters: Int. J. Therm. Sci., 75 (2014), 65.
• 27)   M.  Mohseni and  B.  Peters: Chem. Eng. Sci., 152 (2016), 689.
• 28)   B.  Peters: Fuel, 90 (2011), 2301.
• 29)   B.  Peters: J. Anal. Appl. Pyrolysis, 91 (2011), 352.
• 30)   A. H.  Mahmoudi,  F.  Hoffmann and  B.  Peters: J. Anal. Appl. Pyrolysis, 106 (2014), 9.
• 31)   A. H.  Mahmoudi,  F.  Hoffmann,  B.  Peters and  X.  Besseron: Int. Commun. Heat Mass Transf., 71 (2016), 20.
• 32)   A. A.  Estupinan and  B.  Peters: Int. J. Refract. Met. Hard Mater., 49 (2015), 88.
• 33)   B.  Peters and  F.  Hoffmann: Chem. Eng. J., 304 (2016), 692.
• 34)   B.  Peters,  F.  Hoffmann,  D.  Senk,  A.  Babich,  L.  Hausemer and  J.-P.  Simoes: Metall. Ital., 3 (2016), 49.
• 35)   B.  Peters,  F.  Hoffmann,  D.  Senk and  A.  Babich: Chem. Eng. Sci., 140 (2016), 189.
• 36)   M.  Michael,  F.  Vogel and  B.  Peters: Comput. Methods Appl. Mech. Eng., 289 (2015), 227.
• 37)   K.  Samiei and  B.  Peters: Chem. Eng. Sci., 87 (2013), 234.
• 38)   E.  Copertaro,  P.  Chiariotti,  A. A.  Estupinan,  N.  Paone,  B.  Peters and  G. M.  Revel: Powder Technol., 275 (2015), 131.
• 39)   F.  Hoffmann: Ph.D. thesis, University of Luxembourg, (2014).
• 40)   J. A.  Castro,  H.  Nogami and  J.  Yagi: ISIJ Int., 42 (2002), 44.
• 41)   A.  Faghri and  Y.  Zhang: Transport Phenomena in Multiphase Systems, Academic Press, Boston, (2006), 238.
• 42)   H.  Xiao and  J.  Sun: Commun. Comput. Phys., 9 (2011), 297.
• 43)   K.  Samiei and  B.  Peters: Chem. Eng. Sci., 87 (2013), 234.
• 44)   K.  Samiei,  B.  Peters,  M.  Bolten and  F.  Hoffman: Comput. Chem. Eng., 49 (2013), 183.
• 45)   K.  Samiei: Ph.D. thesis, University of Luxembourg, (2012).
• 46)   H. D.  Baehr and  K.  Stephan: Heat and Mass Transfer, Springer Verlag, Berlin, Heidelberg, Berlin, (2011), 357.
• 47)   S. T.  Miller,  H.  Jasak,  D. A.  Boger,  E. G.  Paterson and  A.  Nedungadi: Comput. Fluids, 87 (2013), 132.
• 48)   H.  Jasak: Ph.D. thesis, Department of Mechanical Engineering Imperial College of Science, Technology and Medicine, (1996).
• 49)   Z.  Zhuo,  H.  Zhu,  A.  Yu,  B.  Wright,  D.  Pinson and  P.  Zulli: ISIJ Int., 45 (2005), 1828.
• 50)   I.-H.  Jeong,  H.-S.  Kim and  Y.  Sasaki: ISIJ Int., 53 (2013), 2090.