The Extented Discrete Element Method (XDEM): An Advanced Approach to Model Blast Furnace

The blast furnace iron making is the oldest but still the main method to produce liquid iron through sequential reduction processes of iron ore materials. Despite the existence of several discrete and continuous numerical models, there is no global method to provide detailed information about the processes inside the furnaces. The extended discrete element method known as XDEM is an advance numerical tool based on Eulerian–Lagrangian framework which is able to cover more information about the blast furnace process. Within this platform, the continuous phases such as gas and liquid phases are coupled to the discrete entities such as coke and iron ore particles through mass, momentum and energy exchange. This method has been applied to the shaft, cohesive zone, dripping zone and hearth of the blast furnace. In this chapter, the mathematical and numerical methods implemented in the XDEMmethod are described, and the results are discussed.


Introduction
A blast furnace is charged with ore and coke from the top, so that a packed bed is formed. Hot gas also referred to as blast is introduced to the furnace through the tuyeres in the lower part. The combustion of coal and coke takes place that generates reducing gas. It ascends mainly in a countercurrent flow reducing iron-bearing materials to form liquid iron and slag in the cohesive zone. The liquids percolate through the coke bed towards the hearth. Unburned coal powder may leave the raceway entrained in the gas flow due to high injection rates. The multiphase flow of solids, liquids, gas and fines in conjunction with reaction kinetics, heat and mass transfer reveals the complexity of a blast furnace and the difficulties to comprehend the process. Since measurements in a blast furnace are very difficult to carry out due to size and limited accessibility, mathematical modelling based on computational fluid dynamics, reaction kinetics and transport phenomena is an applicable tool to investigate into the complex characteristics of a blast furnace.
Multiphase flow phenomena consisting of a fluid, e.g. gas or liquid and a solid phase, could be classified into main categories: On a macroscopic level, all phases are considered as a continuum that has resulted in the well-known two-fluid model [1]. It is a preferred approach due to computational efficiency and convenience for a variety of engineering applications such as process engineering. However, inherent to the averaging concept in the continuum approach, major features of the particulate phase, e.g. material properties, size distribution or shape of individual particles, are lost. Therefore, this gap of information on the particle level usually is compensated for by additional closure or constitutive relations.
Alternatively, the particles are treated as discrete entities, while the fluid phase in the void space between the particles is still considered as a continuum and, therefore, is referred to as the Combined Continuum and Discrete Model (CCDM) [2,3]. It bears the valuable advantage that constitutive relations are avoided and the particulate phase is resolved on a fundamental level. An analysis of results related to particles offers a finer resolution than a continuum method and, therefore, leads to a deeper understanding of particle processes which is confirmed by Zhu et al. [4,5] in a review for the applicability of the CCDM approach. It has developed significantly during the last two decades and represents the dynamics of the particles by the discrete element method (DEM), while the remaining continuum phases are described by differential conversation equations. Hence, CCDM is established as a powerful tool to describe the complex interaction between particles and a fluid as stated by Yu and Xu [6], Feng and Yu [7] and Deen et al. [8]. Although CCDM is the method of choice to reveal underlying physics for challenging multiphase flow applications as reviewed by Zhu et al. [4,5], current methodologies should be developed for nonspherical particle shapes to accommodate efficiently engineering needs and to quantify results for process modelling.
At an early stage, only simple multiphase configurations have been investigated [9,2] that were extended to more and more complex engineering applications such as conveyor belt, cyclone and fluidised bed as demonstrated by Chu and Yu [10]. Chu et al. [11] described the complex flow of water, air, magnetite and coal particles of various sizes in a dense medium cyclone (DMC). Similarly, Zhou et al. [12] were able to model rich/lean combustion of pulverised coal in a complex burner with CCDM. Both predictions agreed remarkedly well with experimental data. CCDM is also well suited to describe the chaotic motion of a particulate phase in a gas-fluidised bed as shown by Feng and Yu [13] and Kafuia et al. [14].
A similar development has been seen for modelling of blast furnaces. Computational fluid dynamics (CFD) as a tool for continuous flow modelling has been applied with success to a large extent as reviewed by Chattopahyay et al. [15,16]. However, experimental observations [17][18][19][20][21] show that a pure continuous approach to the blast furnace is inaccurate. Therefore, the flow of the solid phase consisting of particles is to be modelled by a discrete approach as suggested by Dong et al. [22] and reviewed by Yu [23] for several engineering applications.
Simsek et al. [24] predicted grate ring systems by the CCDM method but obtained only qualitatively reasonable results highlighting the fact that more research is required that should point into polydisperse particle systems [7].
Although CCDM has seen significant progress during the past decade [2,25], applications including heat transfer are rarely attempted. A gas-fluidised bed with heat exchange for polymerisation reactions was undertaken by Kaneko et al. [26] employing the Ranz-Marshall correlation [27]. Only convective transfer was predicted for a two-dimensional spouted bed by Swasdisevi et al. [28]. Additional heat transfer due to conduction between particles during gas-solid transport in horizontal pipes was taken into account by Li and Mason [29][30][31]. Both convective and conductive heat transfer in a gas-fluidised bed for coal combustion were considered by Zhou et al. [32,33]. Similarly, Malone and Xu [34] estimated heat transfer for liquid-fluidised beds with CCDM and concluded that further investigations concerning heat transfer have to be carried out.
A newly developed extended discrete element method (XDEM) introduced by Peters [35][36][37][38] is a coupled continuum and discrete method and a novel approach to answer most challenging questions. The XDEM is an extension to the classical discrete element method with considering thermo-physical state of the particles as well while being in combination with the CFD. It is well described in the mathematical section and applied to different parts of a blast furnace with the following features: • Microscale: Generating a more concrete basis of the particle scale in both dynamic and conversion modules to develop a more comprehensive theory to study the interactions between particles and between particle and fluid under various conditions.
• Macroscale: To provide a general information on the continuum-based process modelling in terms of solving governing equations and to develop a general methodology to link the discrete and continuum approach modelling.
• Application: To present a robust and reliable model capable enough to be applied to the different parts of the blast furnace with transferring data between the multiphase system and particles.

The mathematical model for the continuous phases
The widely used Eulerian volumetric average [39] in the porous media concept was applied to the governing equations of the continuum phases denoted as phasei. In this regard, the conservation of mass, momentum and energy can be summarised as where E, α, r, v, Γ and S are void space, volume fraction, density, velocity, diffusion term and source term, respectively. In order to derive each equation, the variable, ϕ, should be replaced by the values reported in Table 1.
The m 0 , F and H on the right hand side of each equation are the source terms representing mass, momentum and energy transfer between phases. The very well-known Ergun and Ranz-Marshall correlations are used to consider momentum transfer and energy exchange between phases.
The void space and particle diameter must also be known to solve the governing equations for the fluid phases. To this purpose, the XDEM calculates the volumetric average value of void space for each CFD cell by knowing the number of particles in each CFD cell based on a kernelbased interpolation procedure proposed by Xiao et al. [40].

The mathematical model for discrete entities
The XDEM predicts both dynamics and thermochemical conversion of a particulate system. The thermochemical conversion of a particle includes not only the internal temperature distribution but also transport of species through diffusion and/or convection and chemical conversion via reactions in a porous media. The velocity, position and acceleration are calculated in the dynamics module via the DEM method while the temperature and species distribution in the conversion module of the XDEM via the DPM method.

Dynamics module
The proved DEM is based on the soft sphere model applied in the dynamic module. In this method, the particles are assumed to be deformable and allowed to overlap. The motion of each separate particle is tracked using the equations of classical mechanics. Newton's and Euler's second law for translation and rotation of each particle are integrated over time, and the particle positions are updated accordingly during time integration.
• Equation of motion: Here, F ext i ! is the summation of external forces acting on a particle such as drag and buoyancy forces from surrounding fluids. The contact force F c i ! of a particle is the sum of all normal F c, n Hertz-Mindlin model [41,42] is used to calculate the contact forces.

Conversion module
The current approach accounts for thermal conversion of a single particle via heat and mass transfer and chemical reactions. A discrete particle may be considered to consist of a solid, inert, liquid and gaseous phase, whereby local thermal equilibrium between the phases is assumed. The distribution of temperature and species within a single particle is accounted for by a system of one-dimensional and transient conservation equation as shown in Table 2.

Shaft
A blast furnace is a vertical shaft furnace that reduces iron ores to the so-called pig iron for subsequent processing into steels of various qualities. Metallic ore, coke and flux, containing elements such as limestone, are charged in layers into the top of a blast furnace, while usually air is introduced under pressure through the tuyeres (nozzles) at the bottom of the blast furnace. The air is preheated to temperatures between 900 and 1250 C and reacts vigorously with preheated coke to form reducing gas (carbon monoxide). It streams upwards through the void space between the solid materials at high temperatures of app. 1650 C. Hence, a significant reduction of the iron ores takes place in the upper shaft or stack between the top and the cohesive zone of the blast furnace. During indirect reduction with carbon monoxide as opposed to direct reduction with carbon, the iron-bearing materials are stepwise reduced from haematite (Fe 2 O 3 ) to magnetite (Fe 3 O 4 ) and wustite (FeO) to produce liquid iron following three equilibrium reactions: These reactions were experimentally investigated by Murayama et al. [43] under different CO-CO 2 gas mixtures and were employed to predict and validate the individual reaction rate of the abovementioned reduction scheme for a single ore particle. The reduction degree was measured and is defined as the fractional removal of oxygen related to the initially available oxygen as where the total mass of removable oxygen is defined as the oxygen mass to reach the next lower oxide in the sequence of haematite (Fe 2 O 3 ), magnetite (Fe 3 O 4 ), wustite (FeO) and iron (Fe) so that the reduction degree always varies between 1 and 0 for each individual reduction reaction. The reaction rate is expressed by an Arrhenius equation for which the relevant parameters are listed in Table 3.
The predicted results were compared to experimental data at different temperatures and are shown in the Figures 1-3. In general, a very satisfactory agreement between predicted and experimental data is achieved taking into account the temperature spread and experimental uncertainties such as material data. Hence, deviation occurrence is attributed to morphological changes such as porosity or crystal structure. In particular the latter was excluded from the study since crystal and grain structures address a length scale outside the scope of these investigations. Therefore, deviations of this scale are in general acceptable as pointed out by Peters [44].
The validated reaction mechanism was applied to a non-isothermal packed bed of iron ore particles [45] as encountered in the stack of a blast furnace, and the predicted reduction degree was compared to measurements as depicted in Figure 4. Both experiments and predictions were carried out for a packed bed consisting of pellets and pellets/nut coke. For both setups, the predicted integral reduction degree agrees well with experimental data and thus proves the superiority of the current DEM-CFD approach for predicting thermal conversion of packed beds.  Table 3. Kinetic parameters used in the XDEM simulation for reduction reactions.

Cohesive zone
Cohesive zone (CZ) is formed between shaft and dripping zone. In this region, the reduced iron-bearing particles start softening because of the weight of burden above and the high temperature in the CZ. As particles soften, the porosity of the ore layer decreases causing a high pressure drop. As the temperature increases further, the softened particles start melting  and generate two different liquids, molten iron and slag. The softening and melting process makes the CZ to form impermeable ferrous layers forcing gas to flow horizontally through the permeable coke slits. This means that the CZ has remarkable effects on the operation of the blast furnace and the description of the characteristics of the CZ is extremely important. Since it is not possible to interrupt the BF to investigate details of the phenomena occurred inside, the numerical simulation becomes more practical.
The iron-bearing materials start softening at a temperature depending on the chemical composition, reduction degree and the load of burden. In this study, the softening range is considered to be from 1200 K to 1400 K according to Yang et al. [46,47]. To describe the softening process in the CZ, a packed bed is randomly filled with particles assuming pellet which is shown in Figure 5. The radius of pellets varies in the range of 5.4-6.6 mm in a cylindrical vessel with 50 mm in diameter and height. The particles are coloured by their sizes. The value of load acting on the top of particles is determined 100 KPa according to the standard conditions in the softening test [48,49]. For this test case, the conversion and dynamic modulus of the XDEM are used. The conversion module allows to heat up the particles from 1200 K to 1400 K, and the other module is responsible for the movement and the deformation of particle with the assistance of the relationship of Young's modulus (E) and temperature. It was assumed that this relationship is linear within the softening range of pellets. Figure 6 shows the influence of the changes of Young's modulus (E) within XDEM simulation on the rearrangement of pellets in the packed bed. It is observed that to reach %40 in bed shrinkage, it is needed that the value of Young's modulus (E) is decreased approximately by two orders.  Due to the lack of appropriate experimental data for the melting of iron-bearing material by a hot gas in a real condition occurring in a CZ, the melting of a packed bed of ice particles in a solidliquid system is selected. The simulation setup is shown in Figure 7. Figure 8 compares the predicted results with experimental data [51] of the mass loss history for a packed bed exposed to an water mass flow of a velocity of 0.1 m/s. The comparison shows that the predicted mass   variations with time agree well with the experimental results. Figure 9 shows that the ability of the XDEM to fully describe solid and liquid flows in a four-way CFD-DEM coupling considering particle-particle, particle-wall and particle-fluid interaction in terms of dynamics and thermal conversion. As it can be observed, the particles are risen and pushed downwards due to the buoyancy and drag force. Moreover, the particles which are faced the warm liquid earlier are heated up and melted faster than the others. It can be also seen that liquid flows downwards gradually because of the flow resistance from the particles and its temperature decreases by exchanging convection heat with the ice particles and mixing with produced melt. Moreover, The XDEM has been validated for one single ice particle that undergoes melting [52].
Although we modelled solid-liquid flows under the melting process, the present technique is adequately robust and efficient to be applied to the complex melting process of granular packed bed by a hot gas which is a solid-gas-liquid system appropriate for many industrial applications particularly blast furnace modelling.

Dripping zone and hearth
The liquid iron and slag produced in the cohesive zone trickle down through the packed bed of coke particles towards the hearth in a zone known as the dripping zone. It is located deep inside the blast furnace where measurements are not easy to perform. It highly affects the production rate, the quality of hot metal and process efficiency [53]. In this region, the two liquid phases descend slowly and hot gas introduced through the tuyeres ascends upwards while exchanging heat and mass with the particles. The presence of solid particles besides the fluid phases makes the coupled discrete-continuous methods suitable for this application. In our previous study [54], the XDEM approach has been validated by comparing the numerical results with the experimental data [55] for a case representing the dripping zone. The experimental setup and the simulation domain are shown in Figure 10. A liquid distributor is located on the upper part, which is indicated as a mass source for some cells in the simulation domain in order to produce the same amount of liquid as the experimental study. The gas phase is introduced through a nozzle in the experimental study, which is specified as fixed uniform inlet velocity boundary condition in the modelling, while no slip boundary condition was used for the walls and the total pressure for the outlet boundary condition. Results obtained under these conditions are shown in Figure 11. The total absolute mean error of 3.27% was achieved, which shows the robustness of the XDEM for modelling of such cases. This model has also been applied to the trickle bed reactors [56].
In this regard, the model is applied to the lower part of the blast furnace including both dripping zone and hearth shown in Figure 12. The domain is randomly filled with non-uniform particle sizes of d p ¼1AE0:1 cm since the coke particles do not have the same sizes in this particular part. The molten iron and slag produced in the cohesive zone have completely different properties, thus showing different flow behaviours. Most of the previous studies considered only one liquid phase with either slag properties or iron, whereas in this study, both liquid phases are considered, and their mutual effects are studied. Therefore, the simultaneous downward flow of molten iron and slag through a packed bed of coke particles with the upward flow of gas phase is investigated. Despite the existence of several forces, drag force is the only force considered between each pair of phases. In this figure, the calculated porosity field is also illustrated. A predefined cavity is considered right in front of the gas inlet in order to avoid numerical instabilities as shown in Figure 12.
The liquid saturation, which is defined as the volume of the fluid per volume of the void space for both liquid iron and slag, is illustrated in Figure 13 over time. An incline mass source representing    the cohesive zone was specified to investigate its effect on the movement of liquid phases as well as the gas phase. The figures on the left are showing slag saturations and on the right the molten iron saturations. Higher saturations are calculated for the slag, which is due to the higher viscosity and lower density. The slag properties create higher resistance due to the solid particles and in consequence higher saturations. Therefore, the liquid iron flows with higher velocities and lower hold ups through the packed bed. The velocity difference of the two liquid phases causes a slip velocity, which accelerates the slag movement and decelerates the iron downward flow.
Within this concept, the stratification of liquid iron and slag can also be studied. The liquidliquid interface affects the mass and heat transfer rates between liquid iron and slag and thus, determines the final product quality. The molten iron with higher density is settling at the bottom where the slag floats on the top of it.
The effect of the gas phase due to the drag force on the deviation of the liquid phases from the raceway zone, which leads to the formation of the dry zone, is shown in Figure 14. The gas phase tends to replace the liquid phase at the entrance of the gas phase. This behaviour creates a dry zone, whose size is proportional to the gas inlet magnitude [54]. It can also be concluded that the liquid iron is pushed more towards the centre than the slag phase due to the larger drag force between liquid iron and gas phase. Although constant properties are considered for all phases, it should be noted that the fluid properties are not constant while moving through the blast furnace. They are mainly function of temperature, pressure and composition. The gas phase stream lines for gas inlet velocity of 30 m=s are illustrated in Figure 15. The gas velocity loses its energy dramatically as it touches the particles. The recirculation of the gas phase in the raceway besides the redirection of the stream lines in the cohesive zone is observed.

Conclusion
In this study, a coupled discrete-continuous method known as XDEM is introduced, which is able to model different zones of a blast furnace. The governing equations for both continuous and discrete phases are solved through a coupling procedure for mass, momentum and energy equations. The multi-scale and multi-physics characteristic of the XDEM method make it suitable for the complex phenomena modelling such as blast furnace process. The proposed model was applied to different parts of the blast furnace. The results for the shaft, cohesive zone and dripping zone have been validated and compared to experimental data, which shows the reliability of the XDEM method to take into account the whole blast furnace process in the near future.