DEM-

CFDMODELLI
NGOFTHE
I
RONMAKI
NGBLASTFURNACE

Al
l
ertAde
ma
DEM‐CFD Modelling of the
ironmaking blast furnace

PhD Thesis
March 2014

Allert Adema
This research described in this thesis was performed in the Department of
Materials Science and Engineering of Delft University of Technology.

This research was performed under project number MC5.06255 in the frame of the
Research Program of the Materials innovation institute M2i in the Netherlands
(www.m2i.nl)
 DEM‐CFD Modelling of the
ironmaking blast furnace

Proefschrift

Proefschrift ter verkrijging van de graad van doctor

aan de Technische Universiteit van Delft;

op gezag van de Rector Magnificus prof. ir. K.Ch.A.M. Luyben;

voorzitter van het College van Promoties

in het openbaar te verdedigen op dinsdag 4 maart 2014 om 12:30 uur

door

Allert Tjipke ADEMA

ingenieur in de Technische Aardwetenschappen

geboren te Surhuizum
Dit proefschrift is goedgekeurd door de promotor:
Prof.dr. R. Boom

Co‐promotor:
Dr. Y. Yang

Samenstelling promotiecomissie:
Rector Magnificus voorzitter
Prof.dr. R. Boom Technische Universiteit Delft, promotor
Dr. Y. Yang Technische Universiteit Delft, co‐promotor
Prof.dr. H. Saxen Åbo Akademi University, Finland
Prof.dr. ir. J.A.M. Kuipers Technische Universiteit Eindhoven
Dr.ir. T. Peeters Tata Steel R&D, IJmuiden
Prof.dr. J.H.W. de Wit Technische Universiteit Delft
Prof.dr. ir. J. Sietsma Technische Universiteit Delft
Prof.dr. I.M. Richardson Technische Universiteit Delft, reservelid

Keywords: blast furnace, solid flow, gas flow, modelling, coupled DEM‐CFD,
cohesive zone

ISBN: 978‐94‐91909‐06‐1

Copyright ©2014, by A.T. Adema

All rights reserved. No part of the material protected by this copyright notice may
be reproduced or utilized in any form or by any means, electronic or mechanical,
including photocopying, recording or by any information storage and retrieval
system, without written permission from the author.

Printed by CPI‐Wöhrmann Print Service – Zutphen, the Netherlands

Table of contents

Chapter 1 Introduction
1.1. Blast furnace ironmaking 1
1.2. Research question of this thesis 7
1.3. Structure of this thesis 13
References 14

Chapter 2 Blast furnace ironmaking

2.1. Introduction 15
2.2. Overview of integrated steelmaking process 15
2.3. Blast furnace process and equipment 18
2.4. Chemical reactions 19
2.5. Burden materials 21
2.6. Gas and solid flow 24
2.7. Cohesive zone 25
2.8. Process disruptions 26
References 27

Chapter 3 Review of blast furnace modelling

3.1. Introduction 29
3.2. Continuum modelling 30
3.3. DEM Modelling 34
3.3.1. Burden flow 34
3.3.2. Raceway modelling 53
3.3.3. Burden charging modelling 56
3.4. Conclusion 59
References 61

Chapter 4 Development of a DEM‐CFD blast furnace model

4.1. Introduction 65
4.2. Modelling methods 66
4.2.1. CFD Modelling 66
4.2.2. DEM Modelling 71
4.2.3. DEM‐CFD Coupling 73
4.3. Modelling approach 75
4.3.1. Lab scale slot model 75
4.3.2. Pie‐slice model of scaled‐down blast furnace 83
4.3.3. Blast furnace model with different geometries 87
4.3.4. Experimental blast furnace model 91
 4.3.5. Charging and burden descent models 104
4.4. Concluding remarks 112
References 113

Chapter 5 Conclusions and recommendations

5.1. Conclusions 115
5.2. Recommendations 118

Appendix A Parallel development of sub‐models 121

References 128

Summary 129
Samenvatting 131
List of publications 133
Acknowledgements 135
Curriculum Vitae 136
Chapter 1
Introduction

1.1. Blast furnace ironmaking

In 2012 the world production of iron was 1100 million tonnes [1] equivalent to
approximately 90 million tons per month; about 75 % of this was created from pig iron
produced by blast furnace ironmaking. Every second there is a staggering 26 000 kg of
liquid iron (called pig or iron hot metal) produced by this process worldwide. The pig iron is
further processed into steel for use in e.g. the construction, packaging or automotive
industry.
The blast furnace itself is a very large reactor of approximately 40 m high and 15
m wide, producing over 10 000 t/d of pig iron. The process is counter‐current and takes
place at high temperatures. The solids, iron oxides (in the form of pellets, pellets or lump
ore) and coke, are charged in layers at the top of the furnace as shown in Figure 1.1. Oxygen
enriched hot air is blown into the furnace through so called tuyeres at the bottom. Coke
reacts with oxygen to form the reduction gas, carbon monoxide, and part of the heat
required in the process. As the reduction gas rises through the furnace it reacts to reduce
the oxides creating liquid iron and carbon dioxide. Molten iron is tapped from the bottom
of the furnace as well as molten slag which contains the remaining nonferrous materials.

Figure 1.1. Zones in the blast furnace [2]

During the descent through the furnace the solid material charged at the top slowly heats
up and at a temperature of 1100 °C the ore starts softening and melting. The zone in the
furnace where the softening and melting takes place is called the cohesive zone. Below the
cohesive zone only coke is present in solid form with liquid iron and slag trickling down.
Where the blast enters the furnace the Raceway is created; a large cavity created by the
high gas velocity and consumption of coke, with temperatures of approximately 2250 °C. In
the centre of the furnace hearth is the so‐called ”deadman”, a slowly moving pile of coke.

1
Blast furnace ironmaking has been used for centuries and has evolved into an extremely
efficient process. Revolutionary improvements such as the use of coke instead of charcoal
th
in the 18 century are combined with a continuous fine‐tuning of the process. The use of
coke per tonne of hot metal, the amount of reductant required, has nearly halved in the last
50 years as can be seen in Figure 1.2.

Figure 1.2. Decrease in reductant rate

due to technological improvement [3]

A blast furnace today can operate continuously for more than 15 years before it requires
the inner lining to be replaced. In Table 1.1 the evolution in blast furnace size and
production at Tata Steel, IJmuiden (former Hoogovens and Corus) is shown. From the small
Blast Furnace No. 1 to the much larger No. 6 and 7 currently in use, the production has
vastly increased. The final or current productivity of each furnace is more than double of
what is was initially.

The rapid development of China in recent years has greatly increased the global iron and
steel production. Figure 1.3 shows that blast furnace production in China has more than
tripled between 1980 and 2000. Growth then really took off and production grew another 5
times from 2000 to 2012 level, a total of 17 times the 1980 level. During the same period
iron production in the rest of the world grew to a maximum in 2007 with 23%, after which
the global financial crisis caused a severe drop in production. This increase in production
went together with a shift to larger volume blast furnaces. Tables 1.2 and 1.3 show that for
the initial growth years between 2001 and 2006 the distribution of furnace capacity over
the volume ranges still remained the same, but the total number of blast furnaces
increased from 196 to 475. After 2006 the Chinese iron production capacity in blast
3
furnaces with an inner volume of less than 1000 m went from 41% to the current level of
19%. This decrease of share in capacity shifted entirely to furnaces with an inner volume

2
>2000 m3. The number of <1000 m3 furnaces has also significantly decreased from 364 to
206. Compared to the rest of the world China still has a relatively large percentage of the
total production capacity in small furnaces, but is actively restructuring production to larger
and more efficient furnaces [4].

Table 1.1. Development of the blast furnaces at Tata IJmuiden,

The Netherlands [2]
1 2 3 4 5 6 7
Hearth diameter (m) 4.8/5.6 4.8/5.6 5.2/5.9 8.5 8/9 10/11 13/13.8
3
Working volume (m ) 519 519 598 1413 1492 2328 3790
Built 1924 1926 1930 1958 1961 1967 1972
Initial production (t/d) 280 280 360 1380 1700 3000 5000
Most recent production
1000 1000 1100 3500 3700 8000 11000
(t/d)
Most recent production
3 1.9 1.9 1.8 2.5 2.5 3.4 2.9
(t/d∙m )
Last renovation 2002 2006
Demolished 1974 1974 1991 1997 1997

Table 1.2. Percentage of total production capacity

per inner volume range in China [4, 5]
3
Inner volume (m ) >2000 1000‐2000 <1000
2001 38% 21% 40%
2006 38% 21% 41%
2013 59% 21% 19%
Rest of world 2013 72% 23% 4%

Table 1.3. Total number and percentages of blast furnaces

per inner volume range in China [4, 5]
3
Inner volume (m ) >2000 1000‐2000 <1000
# % # % # %
2001 21 11 29 15 146 75
2006 51 11 60 13 364 77
2013 109 27 90 22 206 51
Rest of world 2013 178 43 144 35 89 22

3
 Figure 1.3. Blast furnace iron production [1]

The world’s largest blast furnace is currently POSCO’s Gwangyang No. 1 with an inner
3
volume of 6000 m after finishing relining in June 2013. By using a thinner lining in the blast
furnace the volume of an older furnace can be increased, and in this case the inner volume
3
in 1987 was only 3800 m . Table 1.4 lists the top 15 largest furnaces in the world, and
noticeable is that 13 of these are located in Asia. The shift to larger blast furnaces in China
can be seen from the fact that the three furnaces from China are all built after 2009. There
are currently 27 blast furnaces under construction or scheduled for blow‐in in 2013, 25 of
which are in Asia, 10 in China, 7 in India. Of the 27 furnaces, 17 are big furnaces larger than
3
4000 m . These include what will become India’s largest blast furnace: NMDC’s Nagarnar
No.1 built by Danieli‐Corus with an inner volume of 4506 m3. Tata Steel is building a
3
greenfield steel plant in Odisha: the Kalinganagar No. 1 with an inner volume of 4300 m .

A bigger blast furnace has a larger production capacity for that single unit; it is however, a
trade‐off for productivity. Table 1.5 shows the top 5 highest productivity for blast furnaces
3
>2000 m and the productivity of the top 5 largest furnaces. For productivity the working
volume of a blast furnace is used; this is the volume from stock level to the raceway and
does not include the hearth. The data in the table is limited to furnaces for which the
working volume was available.
The highest current productivity for blast furnaces larger than 2000 m3 is
achieved by Tata Steel Europe’s Blast Furnace No. 6 in IJmuiden with 3.35 tonnes per
3 3
day/m . All top 5 furnaces have productivity above 3 tonnes per day/m , and this is
considerably more than the 5 large furnaces which have an average of 2.5 tonnes per
3
day/m .

4
 Table 1.4. Top 15 largest blast furnaces [5]

Hearth diameter
Inner volume

Last reline
Company

(Mt/year)
Nominal
capacity
Country

Plant

Built
(m )
3

(m)
POSCO S. Korea Gwangyang No. 1 6000 16.1 1987 2013 5.48
Shagang China Zhangjiagang II No. 4 5800 15.7 2009 5.00
NSSMC Japan Oita No. 1 5775 15.6 1972 2009 4.80
NSSMC Japan Oita No. 2 5775 15.6 1976 2004 4.80
POSCO S. Korea Pohang No. 4 5600 15.6 1981 2010 5.31
Severstal Russia Cherepovets No. 5 5580 15.1 1986 2006 3.90
Shougang China Caofeidian No.1 5576 15.5 2009 4.50
Shougang China Caofeidian No.2 5576 15.5 2010 4.50
NSSMC Japan Kimitsu No. 4 5555 15.2 1975 2003 4.53
ThyssenKrupp Germany Schwelgern No.2 5513 14.9 1993 4.30
POSCO S. Korea Gwangyang No. 4 5500 15.6 1992 2009 5.00
JFE Steel Japan Fukuyama No. 5 5500 15.6 1973 2005 4.18
NSSMC Japan Nagoya No.1 5443 15.2 1979 2007 4.25
Kobe Steel Japan Kakogawa No. 2 5400 15.3 1973 2007 3.89
NSSMC Japan Kashima No.1 5370 15.0 2004 4.00

Table 1.5. Blast furnace productivity [5]

Working volume (m)

Nominal production
Inner volume (m )
3

Hearth diameter

(tonnes/day)

Productivity
Company

(tpd/m )
Country

3
Plant

(m)

3
Top 5 productivity for blast furnaces >2000 m inner volume
Tata Steel EU NL IJmuiden No.6 2678 2328 11 7800 3.35
Severstal USA Dearborn No.3 2130 1797 9.2 5890 3.28
Ternium Siderar Argentina Ramallo No. 2 2610 2340 10.4 7200 3.08
Wuhan I&S Co. China Hubei No.5,6,7 3200 2700 12.4 8300 3.07
Tata Steel EU UK Port Talbot No.5 2560 2134 10.8 6500 3.05

Productivity of Top 5 inner volume blast furnaces

Severstal Russia Cherepovets 5580 5200 15.1 11000 2.12
No.5
Shougang China Caofeidian 5576 4670 15.5 12650 2.71
No.1,2
ThyssenKrupp Germany Schwelgern No.2 5513 4769 14.9 12000 2.52
Hyundai Steel S. Korea Dangjin No.1,2,3 5250 4425 14.8 11650 2.63
NSSMC Japan Kimitsu No. 3 4822 4043 14.5 10600 2.62

5
A recent topic which has a large influence on global industry is the reduction of CO2
emissions. This has a large impact on blast furnace ironmaking due to the high amount of
CO2 it produces. Therefore several large projects have been started in a number of
countries. In Japan the governments Cool Earth 2050 initiative has a goal of a 50%
reduction of greenhouse gas emissions by 2050. Approximately 7% of CO2 emissions in
Japan are from the iron making process, so it is an important sector to reduce greenhouse
gas emissions. Part of the response from the Japanese steel industry has been to set‐up the
COURSE 50, standing for: CO2 Ultimate Reduction in Steelmaking process by innovative
technology for Cool Earth 50. The idea behind COURSE 50 is to reduce CO2 emissions by
utilising hydrogen from the coke plant for iron ore reduction in the blast furnace and by
applying CO2 separation to the blast furnace top gas and reusing the CO. [6, 7]
In Europe the ULCOS programme, Ultra Low CO2 Steelmaking, has been
initiated, aiming to reduce CO2 emissions with 50%. The programme is a consortium of 10
companies, of which Tata Steel Europe is one. Four technology options are investigated,
one of which is based on the blast furnace, ULCOS‐BF. The other three are based on
smelting reduction (HISARNA), direct reduction (ULCORED) and electrolysis
(ULCOWIN/ULCOLYSIS). The HISARNA plant is located on the Tata Steel site in IJmuiden
and has had several successful runs.

Figure 1.4. ULCOS‐BF Flowsheet example [8]

The ULCOS‐BF principle is based on top gas recycling, in which the top gas containing
approximately as much CO as CO2 is recycled into the furnace after removing the CO2. The
second part of CO2 reduction comes from Carbon Capture and Storage (CCS); therefore the
efficient removal of CO2 from the top gas is an important part. Experimental campaigns
have been done at the LKAB experimental blast furnace (EBF) in which different process
set‐ups where used, one of which is shown in Figure 1.4. The recycled top gas can be
injected in the raceway with the normal blast or/and at a higher level in the furnace. The
test campaigns have shown that the EBF can be successfully operated with top gas

6
recycling. CO2 emissions in the top gas were reduced by 24% and with CCS a total
reduction of 60% should be possible for an integrated steel plant. [8, 9]

The above methods describe novel technologies to reduce emissions; the first step should
be to increase process and energy efficiency. A highly productive blast furnace using a low
reductant rate will reduce emissions and will also be more cost effective. In Japan NSSMC
(Nippon Steel & Sumitomo Metal Corporation) is increasing the productivity of its blast
furnaces; the same amount of iron can be produced by fewer blast furnaces. This increases
energy efficiency and cost effectiveness per tonne of hot metal produced. Techniques used
to achieve this are DEM simulation and 3D real time monitoring systems which assist in
process optimisation[10].
For a blast furnace to have a high productivity and low reductant rate it firstly
requires high quality input materials. The iron oxide ore, in the form of pellets and sinter,
and coke need to be of an optimal quality, strength and size. When these demands are
fulfilled the blast furnace has to be operated as efficiently and smoothly as possible. This
requires control over the burden distribution and knowledge of the processes inside the
blast furnace, which can be realised through process modelling and simulation, including
CFD‐DEM based approaches.

1.2. Research question of this thesis

Inside the blast furnace a large volume of gas is blown through a slowly descending packed
bed. The balance between those two is critical and has a large influence on the stability and
efficiency of the process. In the cohesive zone the layers of ore slowly soften and deform
after which the melting starts. This greatly reduces the permeability forcing the gas flow
through the untouched coke layers, generating an upward force working against the
particle descent. The cohesive zone also distributes the reduction gas to the upper part of
the furnace where most of the reduction of the iron oxides takes place. Due to the melting
of ore a large solid volume is removed from the furnace, and this is one of the drivers for
the solid flow above the cohesive zone. The softening and melting of the ore influences the
solid flow in the furnace and could potentially create uneven flow or hanging of the
material. Because it is connected to many phenomena critical to the blast furnace process,
the cohesive zone determines the process efficiency and stability of operation to a large
extent [2, 11].
Conditions in the blast furnace make prediction of the behaviour of the cohesive
zone extremely complex; temperatures are very high and there is a large number of
interacting parameters. Observation of the process can only be done indirectly by
measuring the properties of the outgoing material flows e.g. hot metal temperature and
composition. What exactly is happening inside the furnace is impossible to measure
because besides of its sheer size, it is completely filled with a packed bed of solid material
and temperatures go up to 2300C. One method of investigating the inside conditions of
the blast furnace is to quench it and dig it out, like an archaeological dig. In most cases [12]

7
the furnace was quenched by spraying water on the top of the burden; however, during the
resulting slow quenching reactions continue and the burden also reacts with the quenching
water. An example can be seen in Figure 1.5 which shows the results of the dig‐out of the
Sumitomo Metals Corporation Kokura No. 2 Blast Furnace in Japan in 1974. The furnace has
an 8.4 m hearth diameter and was quenched with water for nearly 5 hours to “freeze” the
internal conditions.

Figure 1.5. Dig‐out of the Kokura 2 Blast Furnace [12]

A better and faster way is to quench by injection of liquid nitrogen, which has been done in
only a few cases. Mannesmannröhren‐Werke quenched their Blast Furnace No. 5 (hearth
diameter 8 m) in December 1981 using liquid nitrogen [13]. After quenching the furnace
was excavated and samples were taken in a grid pattern and with vertical probes. Figure 1.6
was constructed using the data from the excavation, both (a) and (b) show the location of
the cohesive zone during the excavation. It is the zone between the line where
temperatures were high enough for melting to start and the line where the ore is molten.
Figure 1.6(b) shows the location of softened and partially molten cohesive layers. By
comparing the reduction degree of the ore at the different grid points an idea of the
intensity of the gas flow can be obtained, this is illustrated in Figure 1.6(a) by the bars in the
furnace.

8
 (a) (b)
Figure 1.6. Dig‐out of Mannesmannrohren‐Werke BF 5
adapted from [13, 14]

Another example of a dig out which also shows the influence of the cohesive zone on the
process was from two trials of a British Steel experimental furnace [14]. The furnace is small
3
with a working volume of only 1 m , but this gives the advantage that after quenching with
nitrogen the furnace can be used again and dig‐outs can actually be compared. The
difference between the two trials was mainly the use of oxygen enrichment in trial B.
Results show a poor cohesive zone permeability in trial A and good in trial B. In trial A some
cohesive bridges had formed shown in Figure 1.7, which were not present in trial B where
melting took place lower in the furnace. The extent of the reduction of the ore above the
cohesive zone was very poor in trial A but good in trial B. This is determined primarily by
the permeability of the cohesive zone. The poor cohesive zone permeability and formation
of cohesive bridges caused operational difficulties during trial A, which were not present
during trial B. There was severe hanging, when the burden does not descend. The cohesive
masses force the gas flow predominantly along the walls of the furnace.

9
 Figure 1.7. Dig‐out of British Steel Pilot blast furnace [14]

These dig‐outs can give a general idea of what the process conditions were in the blast
furnace before shut down. The layer structure of the burden can be seen in Figure 1.5 where
alternating layers of coke and of ore are descending from the top. In the centre of the
furnace the ore and coke are mixed due to the high gas velocity through the centre of the
furnace causing fluidization and local mixing. The high gas velocity can also be seen in
Figure 1.6(a) where the arrows indicating gas flow intensity are much larger in the centre of
the blast furnace. When the ascending gas has heated the burden sufficiently the
descending ore layers starts softening; cohesive masses are formed and the layer
permeability starts decreasing. This is the top of the cohesive zone as shown in Figure 1.5
where the softening zone starts and in Figure 1.6 as the softening line. The cohesive zone
ends when all the ore is molten and only a coke bed remains. This is indicated in Figure 1.5
by the dashed line underneath which only coke is present and in Figure 1.6 by the melting
line.
The cohesive zones from both dig‐outs have a similar shape, high in the centre of
the furnace and sloping towards the walls. This shape is created by the gas flow which is
higher in the centre in the furnace. The cohesive zone in Figure 1.6 curves upward nearer to
the wall, this is an indication of increased gas flow along the walls.
In the cohesive zone the ore layers become cohesive masses which are shown in
both figure 1.5 and 1.6 as black areas, the same can also be seen in Figure 1.7 from the
much smaller furnace. Because the cohesive masses have a low permeability, the
ascending gas flow is forced through the cohesive zone via the coke layers. The influence of
the tuyeres, through which the gas is blown into the furnace, can be seen in Figure 1.5. The

10
two cross‐sections show 4 tuyeres: 140 mm, 120 mm and 70 mm in diameter and one
unused. By decreasing the diameter of the tuyere with an equal pressure drop a lower gas
volume enters the furnace. This results in a lower heat supply and less reductant being
formed in the raceway, consuming less coke. The effect is an increasingly thicker and lower
cohesive zone, with the cohesive zone layers still present below the tuyere.
In the case of the British Steel pilot blast furnace from Figure 1.7, the cohesive
masses block the centre of the furnace. This causes severe operational problems; the
burden descent becomes very irregular or can even stop altogether. Because the gas can
only flow along the wall, there is no reduction taking place in the centre. Since the furnace
is small the effect is rather extreme, but it does show the influence of the cohesive zone on
the process.

These dig‐outs supply valuable information, but have some significant limitations. They are
rare because they can only be done on blast furnaces taken out of production, they are very
expensive and labour‐intensive. They can never give an accurate picture of the conditions
during operation and also do not allow the study of the effect on internal conditions when
changing process parameters. Another method which can be used to study the internal
conditions is by horizontal core drilling from the walls. It has been applied by Tata Steel on
IJmuiden, but this method is limited by the length of samples which can be taken and it can
only be done when the blast furnace is out of operation.

The previous discussion clearly shows the need for a tool which can be used to study the
inner workings of the blast furnace. Mathematical modelling can deliver such a tool and
several models have been successfully used. However, such tools generally require a certain
degree of simplification. One major simplification applied is the simulation of solid particles
by a continuum method. In such models all the individual, discrete particles are replaced by
a single continuous phase with generalized flow characteristics. Simulations are not based
on the theoretical flow characteristics of individual particles but on empirical volume
averaged parameters. Discrete particle flow is rather complicated and cannot be accurately
simulated by these models. What the present research aims to improve over previous
models is the use of a Discrete Element Model to simulate the flow of the individual
particles.

Discrete Element Modelling (DEM) was developed in the late 70’s by Cundall and Strack
[15] but was never extensively applied because of its very heavy computational demands.
Due to the increase in computing power DEM came into practical use in recent years and is
rapidly gaining popularity. In the DEM method all the individual particles are tracked and
for each time‐step the location, velocity and acceleration are calculated. For the present
research the DEM method is coupled with Computational Fluid Dynamics (CFD) which can
calculate the gas flow in the furnace. A coupling module is used to combine both methods
so that e.g. the particles will experience drag force from the gas flow. Two additional
models are required to calculate the chemical reactions and to determine the softening and

11
melting of the ore. Figure 1.8 presents the framework of the full project which is divided
into a PhD and Post‐Doctoral part. This thesis presents the PhD work on the burden flow
model.

Blast Furnace Cohesive Zone Model

DEM + CFD + Thermodynamic + Kinetic models

Burden Flow Model Cohesive Softening and Melt Formation Model

Zone Thermodynamic + Kinetic

Solid Flow Gas Flow Thermodynamic Kinetic Model

Model DEM Model CFD Melt Formation Iron Ore
Model Reduction
Figure 1.8. Project framework

The full project “Prediction of physical and chemical properties of the burden materials in the
cohesive zone”, MC5.06255, was initiated by Materials innovation institute, M2i (formerly
Netherlands Institute for Metals Research, NIMR), and Delft University of Technology (TU
Delft) in cooperation with Tata Steel (then Corus) as industrial partner. The research work
was conducted in the group of Metals Production, Refining and Recycling (MPRR) in the
Department of Materials Science and Engineering of the TU Delft supervised by Prof. dr.
Rob Boom and Dr. Yongxiang Yang. Cooperation with the industrial partner Tata Steel was
with Jan van der Stel and Mark Hattink from the Ironmaking group of the Research and
Development department.

In order to accurately simulate the cohesive zone additional models are required besides
the solid and gas flow models:
 softening and melting can be simulated using a thermodynamic model based on
the chemical composition and changes of the burden
 Chemical reactions taking place in the blast furnace require a kinetic model

Using the combination of the four sub‐models illustrated in Figure 1.8 the final model
should be able to simulate the burden descent in the furnace, how the iron ore is reduced
by the CO in the ascending gas, when the ore starts softening and melting, and what the
thermal profile of the blast furnace is. This results in a model with the ability to simulate the
location and properties of the cohesive zone, and its interaction with the blast furnace
process.

12
For the development of the thermodynamic and kinetic models two post‐doc researchers,
Yuko Enqvist and Vilas Tatavadkhar have worked on both modelling and experiments. The
modelling was focussed on the thermodynamic and kinetic models but also on the
challenging integration of the four models. Experimental work focussed on studying the
softening and melting behaviour of the ore materials.

For both the DEM and CFD models commercial software is used [16, 17], which has several
advantages as well as disadvantages. The software is ready to use and does not require
development from scratch, which also allows for a much easier knowledge transfer within
the research group and to the industrial partner over time. Technical support is easily
available and the software is continuously improved. Because of a larger user base the
accuracy and quality of the results are better determined. Besides high costs the major
disadvantage, also encountered in this research, is the lack of adaptability of the software.
Users are tied to the software as supplied and only a very limited amount of user
adaptation can be made.

1.3. Structure of this thesis

This thesis is built up from 5 chapters. The first gives a short overview of the subject matter
as well as the background and structure of the research project. It explains the drive for the
research and the aims, and where this thesis fits into the larger research project.

Chapter two goes into more detail on ironmaking, explaining the process of producing pig
iron from the raw materials, the process this research project aims to model.

The third chapter gives an overview of the work which has been presented in literature on
the subject of the research project. It shows what has been done on Discrete Element
Modelling of the blast furnace, and where this work fits in and where it differs from
published work.

In the fourth chapter the development of our model and the results are presented and
discussed. Every gradually improved model is applied to a case and also the parallel work in
the research project is discussed. Chapter five contains the conclusions and
recommendations.

13
References
1. Worldsteel Association, Statistics archive. [Accessed 2013, 17‐8]; Available from:
www.worldsteel.org/statistics/statistics‐archive.html.
2. Geerdes, M., Toxopeus, H., and van der Vliet, C., Modern blast furnace ironmaking.
(2009), IOS Press.
3. Ameling, D. and Endemann, G., Ressourcenefficienz: gute argumenten für stahl.
Stahl Eisen 127 (2007), p. 85‐93.
4. Zhang, S.‐r. and Yin, H. The trends of ironmaking industry and challenges to
Chinese blast furnace ironmaking in the 21st century. in Proceedings of the 5th
International congress on the science and technology of ironmaking. (2009),
Shanghai.
5. VDEh, Blast furnaces worldwide. VDEh PLANTFACTS, Date published: 28‐06‐2013
6. Miwa, T. Development of iron‐making technologies. in Proceedings of the 5th
International congress on the science and technology of ironmaking. (2009).
7. Ariyama, T., Ueda, S., Natsui, S., Inoue, R., and Sato, M. Current technology and
future aspect on CO2 mitigation in Japanese steel industry. in Proceedings of the
5th Internatinal congress on the science and technology of ironmaking. (2009).
8. Stel, J.v.d., Louwerse, G., Sert, D., Hirsch, A., Eklund, N., and Pettersson, M., Top
gas recycling blast furnace developments for 'green' and sustainable ironmaking.
Ironmaking and Steelmaking 40 (2013), p. 483‐489.
9. Zou, G. and Hirsch, A., The Trial of the Top Gas Recycling Blast Furnace at LKAB's
EBF and Scale‐up. La Revue de Metallurgie (2009), p. 387‐392.
10. Ueshima, Y., Higuchi, Y., and Saito, K. Recent topics of iron‐ and steelmaking
technology in NSSMC. in 2013 Annual review of NSSMC. (2013), McGill Metals
Processing Centre (MMPC).
11. Biswas, A.K., Principles of ironmaking. (1981), Cootha Publishing House.
12. Omori, Y., ed. Blast furnace phenomena and modelling. (1987), Elsevier Applied
Science.
13. Engel, K., Fix, W., Grebe, K., de Haas, H., and Winzer, G., Ergibnisse von
Untersuchungen an dem mit Stickstoff abgekülten Hochofen 5 der
Mannesmannröhren‐Werke AG. Stahl und Eisen 106 (1986), p. 18‐23.
14. Gathergood, D., Jones, J., Juckes, L., and Golding, D., Progressive reduction of
burden in the blast furnace, in Technical Steel Research. (1996), British steel,
Teeside Technology centre.
15. Cundall, P.A. and Strack, O.D.L., A discrete numerical model for granular
assemblies. Geotechnique 29 (1979), p. 47‐65.
16. DEM Solutions, EDEM. Available from: www.dem‐solutions.com.
17. ANSYS, ANSYS Fluent. Available from: www.ansys.com.

14
Chapter 2
Blast furnace ironmaking

2.1. Introduction
The first iron objects were produced from meteoric iron and the oldest known examples are
iron beads from Gerzeh (3500 BC) and a dagger from Ur (3000 BC). Smelted iron objects
appeared around 2000 BC and from 1000 BC iron began to replace bronze as the most
important material for tools and weapons. Early iron smelting was done in a bloomery
furnace, a batch process in which the ore does not become liquid but remains solid while
being reduced to iron. The end product was a large block of iron called the bloom which
was then further processed by hammering out the impurities. The earliest blast furnaces
producing liquid iron were built in China as early as 200 BC. In Europe the process of
ironmaking in a blast furnace was first used in the middle ages, since then the process has
undergone large changes but is essentially the same. It is the first part of a two stage
process to produce steel, the first stage produces crude pig iron which is then refined to
steel in the second stage. The name pig iron comes from an old casting process in which
the crude iron was cast into a sand bed with a centre runner and a number of side runners.
When the iron is solidified it somewhat resembled piglets drinking from their mother pig.
In the medieval process the blast furnace was fed with iron ore and charcoal.
Water wheels were used to drive the bellows which pumped air into the furnace as well as
to drive hammers for processing of the steel. During the Industrial Revolution the demand
for steel increased dramatically with the development of steam engines, railways, bridges
and machinery. Blast furnace iron production rapidly expanded and innovated to meet this
th
demand. In the 18 century coke replaced charcoal and steam engines replaced water
th
wheels, the 19 century brought hot air blast. Continuous research and development
greatly increased the size, efficiency and productivity of the blast furnace. Table 2.1 shows
the development of blast furnaces from the industrial revolution to the modern age.

The current process of ironmaking is described in this chapter starting with the place of the
blast furnace in the steelmaking process. Focus will then shift to the blast furnace itself and
its operating parameters.

2.2. Overview of integrated steelmaking process

Blast furnaces are either located on what is called an integrated steel plant or feed different
steel plants at some distance from the site. On an integrated steel plant the whole process
from iron ore to steel products takes place as shown in Figure 2.1. The process starts with
the unloading of the raw materials from ships or trains. These materials are stored, mixed
and pre‐processed for use in the blast furnace. The raw materials consist of two types: the
iron containing ferrous oxides and the carbon containing reductants. Ferrous oxides can be
iron ore or processed pellets or sinter. Production of sinter generally takes place on‐site;

15
pellets are also commonly processed at the mine site. Only two ironmakers produce pellets
at the integrated site: Tata Steel IJmuiden (Netherlands) and Kobe Steel Kakogawa
(Japan). This has some big advantages: better control over pellet composition, direct
quality control on‐site, pellets are warm and dry, and less handling results in less pellet
degradation. The disadvantages are stockpiling of pellet ores, good quality control is
required and the fact that an extra operation is required on‐site. The large majority of
carbon containing material is coal but oil or gas can also be used in small quantities. Coal is
processed to coke by removing the volatile matter and carbonizing in coke batteries.
Unprocessed coal fines are directly blown into the furnace with the blast as Pulverized Coal
Injection (PCI). To this end large coal grinding facilities are needed on the site.

Table 2.1. Development of blast furnaces

Year and location Hearth Reductant Iron output,
diameter, rate, tonne per day
m2 kg/tonne of
hot metal
1796 (Gleiwitz, Germany) 0.6 3500 1‐2 (later 4)
1801‐1815 (UK) ‐ 2500 5‐7
1856 (Hasslinghausen, 2.1 1600 20‐23
Germany)
1880 (USA) 3.4 1510 120
1902 (USA) 4.4 1000 464
1929 (Bruckhausen, Germany) 6.5 740 1100
1993 (Schwelgern 2, Germany) 14.9 480 10600
2013 (Pohang 4, Korea) 15.6 ‐ 14600

Figure 2.1. Steelmaking process [1]

16
The ferrous oxides and coke are then charged into the blast furnace where the ferrous
oxides are reduced, thereby producing liquid slag and hot metal. Slag contains the
unreduced components from the burden and consists mainly of SiO2, CaO, MgO and Al2O3.
It is much more viscous than the hot metal and its composition is very important for the
process. Slag and hot metal are almost continuously tapped from the furnace from a
number of tapholes at the bottom of the furnace. The slag is solidified and when it meets
certain requirements, can be used in e.g. road construction or as raw material for cement
making for building purposes. The liquid hot metal is poured into large torpedo‐shaped
ladles, elongated firebrick insulated rail cars called torpedo cars, for transport. Because the
liquid iron is in constant contact with coke inside the blast furnace it contains a large
amount of dissolved carbon, much higher than required for steel. High carbon content
makes the iron very hard but also very brittle and it has to be reduced from approximately
4.5 wt% to 2 wt% ‐ 0.05 wt% depending on the type of steel. This is done in the Basic
Oxygen Steelmaking (BOS) process, where the hot metal is poured in a converter into
which a water‐cooled lance is used for oxygen blowing. The lance blows oxygen with a
velocity higher than Mach 1 onto the hot metal burning off the dissolved carbon. Steel
scrap is charged into the furnace before the liquid iron to recycle the steel and use the heat
created by the burning of the carbon for melting the scrap. Besides reducing the carbon
content the BOS process also removes silicon, titanium and phosphorus. Further impurity
removal can be done in the steps before or after BOS e.g. desulphurization in the torpedo
or charging ladle. Typical batch sizes of the converter process vary between 100 and 350
tonnes.
After finishing the oxygen blow the steel is tapped from the converter into a steel
ladle. A second step of refining follows in which alloying elements can be added and
impurities removed to fine‐tune the composition to match specifications. Alloying
elements can be added by feeding a wire into the melt or by injecting a powder. Excess
dissolved oxygen is removed by deoxidants such as aluminium or silicon. Before casting it is
ensured that the liquid metal has a homogenous temperature and composition by stirring
the metal using argon injection or Electromagnetic Stirring (EMS). The temperature can be
adjusted by scrap cooling or heating in a ladle furnace.
Steel is continuously cast into slabs, blooms or billets with different dimensions
depending on the final product, e.g. slabs 1600 mm wide, 250 mm thick and 12 m long. The
hot metal is cast into a large tundish, a holding reservoir and distributor which can typically
contain 30 to 80 tonnes of liquid steel. From the tundish the hot metal is fed into several
vertical water‐cooled copper moulds which solidifies the outer layer of metal. The metal
exits the mould supported by rollers and still has a molten core; water is sprayed on the
surface for secondary cooling. The next step in the process is the hot or cold rolling. The
former is sometimes placed directly after the continuous caster to take advantage of the
heat remaining in the metal. During hot rolling, the metal passes between rolls reducing
the thickness and width. By carefully managing temperature and cooling rates the physical
properties of the steel can be controlled. The deformation by successive cold rolling causes
hardening and increases the strength of the steel and the roughness of the steel surface.

17
2.3. Blast furnace process and equipment
The blast furnace is a counter‐current reactor, the ferrous oxides are charged at the top and
the reducing gas is blown in at the bottom. Solids are charged in alternating layers of coke
and ferrous oxides to ensure good permeability when the oxides start melting. At the
bottom of the furnace oxygen enriched and pre‐heated air is blown in at high velocity
through the tuyeres. Gas pre‐heat temperatures are approximately 1200 °C and the flame
where the oxygen reacts with the coke has temperatures in excess of 2100 °C. Additional
reductants such as pulverized coal, gas or oil are also injected with the hot blast, reducing
the amount of more expensive coke required. The Tata blast furnaces at IJmuiden can
currently achieve coke rates of 270 kg per tonne of hot metal with pulverized coal injection
(PCI) rates of 230 kg/tHM. The lowest sustainable rate achieved in IJmuiden for a week was
250 kg/tHM at which more than half of the reductants are injected. Yearly averages of 260
kg/tHM have been reached. The gasification of the coke by the blast causes a large cavity
called the raceway.
As the solids slowly descend they are heated by the ascending reductant gas. A
large part of the oxides are reduced in the upper part of the furnace with the remainder
being reduced in the cohesive zone. As the burden reaches temperatures of around 1100 °C
the ore starts softening and melting. In the cohesive zone the permeability of the ore layers
decreases and they become nearly impermeable. All of the ore melts and only coke
remains solid below the cohesive zone. During melting two liquid phases are created, the
liquid iron and the slag, trickling down through the coke bed into the hearth. The liquid iron
contains dissolved carbon as well as some impurities which are reduced besides the iron
such as silicon, manganese, titanium, sulphur or phosphorus. All the unreduced oxides from
the ore and fluxes form a silicate slag.

Around the blast furnace there are several pieces of equipment required for the process
shown in Figure 2.2. The oxygen enriched air is pre‐heated in the hot blast stoves: these are
large cylinders filled with fire bricks. A burner inside the hot blast stoves burns furnace gas
and/or natural gas to heat up the bricks. When the required temperature is reached the
burner is turned off and air is blown through, whereby the heat stored in the bricks is
transferred to the air. Using multiple stoves the gas is heated continuously. The hot blast is
blown into the bustle pipe which runs around the furnace and then injected into the blast
furnace via the tuyeres. At the top of the furnace the top gas is removed via the uptakes
and a down‐comer into a dust catcher for the removal of fine particles.

The solid materials are stored in the stock house where they are weighed and screened
before being charged into the blast furnace. Material is charged into the furnace by skips or
conveyor belts which bring the material up and by a charging mechanism which deposits
the particles into the furnace while keeping the furnace under pressure. There are two main
methods: a bell top or a bell‐less top. In the former the particles are charged into a series of
bells; the charge can be distributed over the material surface in the blast furnace by

18
movable armour plates which deflect the material. A bell‐less top uses a rotating chute to
distribute the material and is far more accurate in the distribution of the burden.

Figure 2.2. Blast furnace and equipment [2]

2.4. Chemical reactions

In the blast furnace the iron oxides present in the ore are reduced to metallic iron by the
carbon present in the coke and additional reductants. In the raceway the carbon from the
coke and coal reacts with oxygen to form the main reductant gas CO.

O2 + 2C → 2CO (2.1)

The iron oxide is present in the form of Hematite (Fe2O3) and the reduction takes place with
a decreasing ratio of O to Fe, forming Magnetite (Fe3O4) and Wϋstite (FeO). If we follow the
descending burden through the furnace the following reactions take place:

Hematite 3Fe2O3 + CO → 2Fe3O4 + CO2 (2.2)

Magnetite Fe3O4 + CO → 3FeO + CO2 (2.3)
Wϋstite 2FeO + CO → 2FeO0.5 + CO2 (2.4)

These reactions take place in the upper part of the furnace at increasing temperatures.
Hematite starts reducing at about 500 °C, magnetite at 600 °C to 900 °C and the reduction
of wϋstite at 900 °C to 1100 °C. At the start of the melting temperature of 1100‐1150 °C
FeO is generally reduced to FeO0.5.

19
Two types of reduction reactions take place in the blast furnace: indirect and direct
reduction. The reactions above where CO reacts to form CO2 take place by indirect
reduction; the coke is not directly involved in the reaction. Final reduction of FeO0.5 to Fe
takes place by direct reduction in a solid‐liquid reaction. In direct reduction the iron is
reduced directly by the coke and can be written as:

2FeO0.5 + C → 2Fe + CO (2.5)

However, the reaction actually takes place in two parts:

2FeO0.5 + CO → 2Fe + CO2 (2.6)

CO2 + C → 2CO (2.7)

The iron oxide reacts with CO gas forming CO2 which then immediately reacts with carbon
creating CO. This last reaction is called the Boudouard reaction and describes the balance
between CO on the one hand and CO2 and C on the other hand. Figure 2.3 shows both the
temperature dependent equilibrium relationship of the Boudouard reaction and the Fe‐O‐C
equilibrium diagram. Above 1100 °C all CO2 reacts to form CO, at lower temperatures an
increasing amount of CO2 is present. If CO is decomposed very fine carbon is formed with
CO2, however, below 400 °C the reaction is very slow and only a negligible amount of
carbon forms. Figure 2.3 shows the reduction of magnetite to wϋstite requires only 40% to
20% CO/CO+CO2, depending on the temperature, compared to 60% and higher for the
further reduction of wϋstite to metallic iron.

Figure 2.3. The Fe‐O‐C equilibrium diagram combined

with the Boudouard curve [3]

The furnace can be divided into three zones: upper, middle and lower. In the lower zone
molten material temperatures reach around 1500 °C and the approximately 2100 °C gas
cools down to 900 °C – 1000 °C. All the oxygen reacts with coke to form CO and any formed

20
CO2 reacts back according to the Boudouard reaction. The reaction of oxygen with the coke
generates a very large amount of energy and supplies the heat to the furnace. Direct
reduction of iron oxide consumes a large amount of energy and thus significantly cools the
lower furnace.
In the middle zone the temperature remains relatively constant at about 900‐
1000 °C and is called the thermal reserve zone. The effects of endothermal and exothermal
reactions, burden heating and heat losses cancel each other. In the thermal reserve zone
the wϋstite is slowly reduced in equilibrium with CO as shown in Figure 2.3. By probing it
has been shown that inside the thermal reserve zone there also is also a chemical reserve
zone where no reactions take place and Fe/FeO and CO/CO2 are in equilibrium. The thermal
reserve zone will shrink as the furnace is pushed to higher productivity. The upper zone is
where the gas temperature drops from 900‐1000 °C to 100‐250 °C and the solids heat up to
800 °C. Hematite and magnetite are reduced to wϋstite and any moisture in the burden is
vaporized.
Besides the CO, hydrogen is present as a second reductant in the blast furnace. At
temperatures above 821 °C it has a higher reduction efficiency than CO. Hydrogen is
formed from the moisture in the blast and in the pulverized coal and can then reduce the
iron oxides:

H2O + C → CO + H2 (2.8)
H2 + FeO → Fe + H2O (2.9)

Water generated in the furnace can react with CO to form hydrogen:

H2O + CO → CO2 + H2 (2.10)

Reaction 2.8 only takes place at temperatures above 1000 °C. At lower temperatures the
reaction 2.10, called the water‐gas shift reaction, shifts strongly to the right. Because it
requires iron as a catalyser the water‐gas shift reaction does not take place in the upper
part of the furnace.

2.5. Burden materials

A blast furnace has two main components in its top charged burden: the raw materials that
contain iron oxide and the materials that contain carbon. Three types of oxide materials,
shown in Figure 2.4, are used: lump ore, sinter and pellets. Lump ore is used directly from
the mine where it is screened to get the required size fraction of approximately 6‐25 mm.
Lump ore is cheaper than pellets but it has poorer properties for the blast furnace. The ore
is weaker than pellets during reduction under the high temperature and pressure in the
blast furnace. Because it is naturally occurring it has a wider range of compositions and
physical properties. Lump ore is generally only used in small percentages of the total oxide
feed.

21
 Sinter Pellets Lump ore
90% < 25 mm 11 mm (±2 mm) 6‐25 mm
Figure 2.4. Blast furnace ferrous burden materials [2]

Sinter is a blocky material made by fusing iron ore fines into larger particles. Originally
applied to use revert materials from the whole steel making site in the blast furnace, it is
currently used much more widely as one of the predominant sources of oxides beside
pellets. The main component of the sinter is iron ore fines, either raw material or a return
flow from the process. Return sinter from the sintering process is also added. They are
mixed with limestone as a flux to ensure a final slag in the blast furnace with the desired
properties for the ironmaking process and for the use as raw material for cement
production. The final component is coke fines added as a fuel for the sintering process.
After mixing in a rotating drum where 5 % water is added for primary binding the mixture is
charged on a sinter strand in a 35‐65 cm thick layer. A sinter strand, shown in Figure 2.5 is a
moving grate through which air is being sucked from the bottom. At the start of the strand
the coke in the bed is ignited using heated air. The heat generated partially melts and
sinters the material. As the sinter moves over the strand the flame front moves down
through the bed until it reaches the bottom and all of the material is sintered. At the end of
the strand the material breaks off and is further reduced in size.

Pellets are created from very fine iron ore (much finer than used for sinter), which has been
pre‐processed to separate the iron oxides from gangue materials. These fines are too small
for sintering because the bed would be impermeable. Similar to the sintering process the
iron ore fines are mixed with fluxes and coke fines in a rotating drum. Then water and a
binder are added, the rotation causes snowballing and creates spherical particles. The
particles, called green pellets, are then charged onto a grate moving through a furnace and
heated to 300‐350 °C. This removes the water and the binder creates a chemical bond to
hold the particle together. Burners and the coke fines then increase the temperature to
1250‐1350 °C fusing the particle.

22
 Figure 2.5. Sinter strand [1]

The second component in the burden is the coke which is supplying carbon for the
reduction and the heat to melt the burden. Coke is produced by mixing and grinding coal to
the required specifications which is then charged into a coke oven. In the oven the coke is
heated to over 1200 °C in an air‐free environment. The volatile matter in the coke escapes
and a solid carbon matrix is formed by the carbonization of the coal. After 16 to 24 hours all
the material is converted and the hot coke is pushed out and immediately quenched with
water to ensure it does not react with oxygen.
With a size of 45‐55 mm coke is larger than the other burden materials. Coke
quality is very important for the process; the whole burden is supported by the coke and it
has to ensure a good permeability. A large amount of coke can be replaced by injecting coal
directly in the raceways, by PCI. The quality demands on pulverized coal are much lower
and it is therefore cheaper. There is a limit on the amount of coke which can be replaced by
PCI due to the permeability requirements of the burden. This increases the quality demand
of the burden materials on strength and degradation resistance.

Blast Furnace No. 6 in IJmuiden has achieved very low coke rates of 260 kg/t hot metal year
average, globally a more usual rate is above 300 kg/tHM. To achieve such low coke rate the
burden permeability has to be as high as possible. For this the pellets and sinter need to be
of good quality; both need a narrow size distribution to maximise burden porosity and the
particles should not break or degrade during descent. The amount of slag generated should
be low, as slag in the inter‐particle pores will decrease the permeability. Because the gas
flows through the packed bed of ore and coke layers, the structure of the layers is very
important to control the permeability. This also influences the cohesive zone shape and
location, which has a very large influence on the pressure drop in the blast furnace. Good
control of the burden charging and the resulting layer structure is crucial.

23
The choice of raw materials for the blast furnace is highly dependent on their costs and
availability as well as their properties. As the burden descends through the furnace the load
onto the particles becomes very high. The strength of the particles determines the amount
of breakage and attrition occurring. Burden composition determines the reduction
parameters as well as the softening and melting temperatures and behaviour. The size
distribution greatly influences the permeability of the bed; blocky irregular sinter has a
higher permeability than spherical pellets. Every blast furnace has its process specifically
balanced to the raw materials it uses.

2.6. Gas and solid flow

In the research presented in this thesis the main focus is on the solid and gas flow. Because
of the large size of the blast furnace and the complex character of the process, the
behaviour of the burden is hard to predict. A very important feature of the blast furnace
solid flow is the way the material is charged. The main characteristic is that the material is
charged in alternating layers of coke and ferrous oxides. A completely mixed burden would
become impermeable to gas flow if the ore starts melting. When charged in layers the gas
can flow through the coke slits if the oxide layers become impermeable. The composition
of the oxide layers can be varied by pellets, sinter, lump ore or mixtures to get a repeating
series of oxide layers in the furnace. Small sized coke is sometimes mixed in the ore layers
and has several advantages; it can increase layer permeability when the ore starts softening
and melting , and it can improve reduction of the iron oxides by supplying reductant inside
the layer. [4]
The amount of material charged can also be distributed in the radial direction.
This allows a change in thickness of the layers; a common practice is to charge more coke in
the centre of the furnace creating a highly permeable core. Because more of the hot gas
flows now through the centre, it determines the inner temperature profile and thus the
location of the cohesive zone and reductions.

Figure 2.6 shows the difference in cohesive zone shape between a furnace with a coke
centre and one where more coke is charged at the walls. Both methods have their
advantages and disadvantages. Clearly the distribution of the burden determines the gas
flow in the furnace, which in turn determines where the reduction takes place and how the
cohesive zone is formed. The permeability of the individual layers is determined by the size
distribution and shape of the particles. A layer of spherical pellets has a lower permeability
than a layer of irregular, blocky coke. If the size distribution increases smaller particles tend
to fill the gaps between larger ones and the permeability decreases. A high permeability is
preferred to allow good distribution of the reductant gas through the furnace.

24
 Figure 2.6. Two types of cohesive zones: central working (left)
and wall working (right) [2]

The burden and solid flow profiles in the furnace are for a great part determined by the
consumption of solids; they flow to where there is space created. There are a number of
processes which consume solids in the furnace of which the largest is the melting of the ore
in the cohesive zone. This accounts for over 55 % of the solid volume consumed in the blast
furnace. The second and third are coke gasification in the raceway and in direct reduction
via the Boudouard reaction around the cohesive zone, both at approximately 15‐20 %.
Carburization of the liquid iron is the fourth solids consumer at around 5 %.

2.7. Cohesive zone

Focus of the larger project of which this thesis is a part, is the cohesive zone. It starts when
the oxide layer begins softening and ends when all solid iron oxides have melted. The
cohesive zone is highly complex due to the large number of interacting processes taking
place. Beside softening and melting there are a number of reactions, most importantly the
final reduction of the iron oxides combined with the Boudouard reaction.
Softening and melting temperatures are highly dependent on the chemical
composition of the ore, sinter or pellets as well as the reduction degree. Softening
temperatures range from 1000 °C to 1100 °C and melting temperatures from 1250 °C to
1450 °C. When the particle starts softening, its shape changes due to the pressure from the
load above it. The layer height reduces and the pores between the particles shrink. Melt
formation fills in the remaining pore space. Both softening and melting cause cohesive
forces between the particles and the layers more or less become a sticky mass.
Pellet softening and melting is complicated due to its round and homogenous
shape. The reduction takes place from the outside of the pellet inwards as shown in Figure
2.7. When melting temperatures are reached the outside of the particle is an iron shell with
a molten core due to the lower melting point of the oxides inside. Pure iron has a melting
point of 1539 °C and wϋstite of about 1380 °C. Melt formation from the pellet occurs when

25
the iron shell is sufficiently weakened by deformation and cracks. Because sinter and lump
ore have a much more heterogeneous composition there is no iron shell formed.

Figure 2.7. Pellet reduction [1]

Due to the low permeability of the ore layers gas can only flow through the coke slits. The
cohesive zone distributes the reductant gas to the upper furnace where most of the
reduction takes place. This makes the cohesive zone very important for an efficient and
stable process. If the ore layers touch each other in a section of a furnace, that area
becomes impermeable to the reductant gas and the oxides above it will be reduced much
more slowly. This causes an inhomogeneous process in the furnace and thus reduced
productivity. Another risk is that low permeability causes a gas pressure build up under the
cohesive zone. The upward force can become larger than the weight of the burden and the
burden above starts hanging. Below the cohesive zone the flow continues and a void is
created, at some point the balance will be disturbed and the burden will drop down into the
void. This is a very large furnace disruption and could have large consequences.

2.8. Process disruptions

An efficient and stable process can only take place when the gas can flow through the
furnace with minimum resistance. Increasing resistance to gas flow can be caused by
process disruptions or can in turn cause process disruptions. If the resistance becomes too
high the pressure builds up causing hanging and slipping of the burden and this can
severely disrupt the burden layer profile. Recovery from such events generally requires
decreasing production and increasing coke rates, as well as sufficient time to replace the
burden.

26
References
1. Worldsteel association, Steeluniversity.org. [Accessed 2013, 05‐06]; Available
from: www.steeluniversity.org.
2. Geerdes, M., Toxopeus, H., and van der Vliet, C., Modern blast furnace ironmaking.
(2009), IOS Press.
3. Biswas, A.K., Principles of ironmaking. (1981), Cootha Publishing House.
4. Song, Q., Effect of nut coke on the performance of the ironmaking blast furnace,
(2013), PhD Thesis at Department of Materials Science and Engineering, Delft
University of Technology.

27
28
Chapter 3
Review of blast furnace modelling

3.1. Introduction
The blast furnace described in chapters 1 and 2 is one of the most widely used packed bed
reactors in the world; it has been for many years and will remain so for a long time. Inside
the furnace several types of flow are present: a packed bed of solids descending, liquid
dripping in the lower zone and gas with powder (dust, pulverized coal) ascending through
the packed bed. A stable and efficient operation of the furnace depends for a great part on
the interaction of these flows. Because the flows are counter current they influence each
other’s flow pattern as well as react with each other. Because of the high complexity of the
processes taking place inside the furnace and the impossibility of ‘looking inside’ due to
high pressure and temperature, computer simulation offers a way to gain knowledge about
the process. Therefore, the modelling of the material flows has received considerable
attention.

The models used to simulate the blast furnace flows can be divided into two types based on
the method applied to describe the solids flow: continuum models and discrete element
methods. Continuum models are based on the finite element or finite volume method in
which a geometry is divided into a finite number of cells over which a set of equations are
defined describing flow of scalars such as velocity, mass and temperature. These are solved
numerically to achieve a solution for each cell and time step if transient. In each cell the
mass, velocity, temperature, etc. of that specific cell is stored. This allows a graphic
representation of the solution to be made.

Discrete element or discrete particle models are based on the movement of separate
particles. The method stores the location, velocity vector and acceleration vector of each
particle and uses these to calculate where each particle will be the next time step. When
two or more particles collide it uses a contact model to calculate the forces working on the
particles and how the collision changes the velocity and acceleration vectors. The discrete
element method is by definition a transient model. Because of the high computational
requirements of the discrete element method it has only become widely used in recent
years in contrast to the continuum method which has been in use for the last decades.

In this study a combination of both methods is used; a continuum model for gas flow and a
Discrete Element Method (DEM) for solid flow. Though the continuum method is not used
for solid flow modelling, its use in that area is briefly discussed here for comparison and
completeness. After that the use of the discrete element method for blast furnace
modelling is described.

29
3.2. Continuum modelling
Models of the blast furnace have been developed using the continuum method since the
80’s. Yagi [1] has written an extensive review on the work done since that period, in which
general simulations of packed bed flow phenomena are described by the flow of four fluids
representing gas, liquid, solid particles and fine powder. The flow is mathematically
modelled by the generalized equations for continuity and the Navier‐Stokes equations for
motion. As this is a continuum approach based on local averaging, the four fluids are
defined by the volume fraction of the continuum they occupy as shown in Figure 3.1(a). The
gas phase G fills the voids of the packed bed, and the flowing powders are regarded as the
powder phase F. The liquid phase L contains the flowing liquids as well as the amount of
powder entrained in the liquid εf(l,d). The solid phase S consists of the solid packed particles
as well as static hold‐up consisting of powder εf(l,s), liquid εl(s) and powder entrained in the
liquid εf(p,s). The four fluids affect each other through exchange of momentum, mass and
heat as illustrated in Figure 3.1(b). This is defined by the interaction parameters, , for each
of which there are underlying equations describing the interaction, e.g. drag force by gas
flow on the powder.

The source term in the continuity equation includes sources due to chemical reactions. It
expresses the amount of gas phase increase due to phase changes and absorption or
desorption at the surface of liquid, powder and packed interface. The source terms for fine
particles define powder generation due to physical disintegration and chemical reactions as
well as powder entrainment by the liquid and solid phases. For the liquid phase the source
terms express liquid generation due to melting, chemical reactions and powder
entrainment. The packed particle source terms describe mass generation and loss due to
chemical reactions, melting and physical disintegration as well as static hold up of liquids
and powders.

(a) (b)

Figure 3.1. (a) Four fluids in the packed bed,

(b) Interaction parameters among four fluids, [1]

30
 Figure 3.2. Simulation results of timelines for packed particles descent by gravity [2]

The paper from Yagi describes two models developed to simulate solid flow by the
continuum method: the kinematic and the potential flow model. These methods can, up to
a certain extent, predict simple solids flow but rely on empirical data of the bulk. An
example of the potential flow model is shown in Figure 3.2 where the solid flow is defined
by a solid viscosity µs and compared to experimental results.

Austin, Nogami and Yagi have presented a model of an industrial scale blast furnace based
on the four fluid model including heat transfer and reactions, [3, 4], some results of which
are shown in Figure 3.3. The method was then applied to several cases on static liquid hold‐
up [5], scrap charging [6] and top gas recycling [7]. The model was expanded by De Castro
[8] to include separate hot metal and slag phases and used to do transient simulations. It
further included separate phases for pulverized coal and pre‐reduced ore [9]. The model
was developed for 3D simulations with five phases [10].
The model has since been refined and applied to several cases: carbon composite
agglomerates charging and top gas recycling, [11, 12]; hydrogen bearing materials, [13];
novel feed materials [14, 15]; char and fine coke [16]

Zaïmi presented a model based on Yagi’s four fluid model but using the hypo‐plasticity
model to simulate the solids flow in the blast furnace [17]. This model has the great
advantage of being able to predict the deadman instead of it being predefined, which is
required in the original four‐fluid model. In this model the deadman is not static but a zone

31
of very low moving solid and not a boundary for solid flow. Further work and experimental
validation were done [18].

(a)

(b)
Figure 3.3. (a) CO and CO2 gas phase mass fraction [3],
(b) Velocity vectors [4]

Zhang published a blast furnace model including solid flow based on the kinematic model
and gas flow [19]. A model to investigate the dripping zone of the blast furnace was
presented by Wang based on a discrete fluid flow model to simulate the interaction

32
between gas and liquid flow including an impermeable cohesive zone [20]. The model
shows that the shape of the impermeable cohesive zone has considerable influence on the
liquid flow. Chew presented a model including gas, liquid and solid flow; heat and mass
transfer; and chemical reactions [21]. The main focus of this model is the liquid flow which
is shown to be highly influenced by the packed bed profile and the gas flow.

Zhang proposed a model to simulate solids flow in a small scale blast furnace to improve
the earlier kinematic and potential flow models as an extension of the earlier work [22]. The
work includes solids consumption and a layered burden, an example of the work is shown in
Figure 3.4. An investigation on gas‐powder flow in a packed bed was published by Dong
[23, 24]. The continuum method used was able to predict powder accumulation and the
effect of the presence of an impermeable cohesive zone.

A more detailed in‐depth review of these and other blast furnace models was published by
Dong [25].

Figure 3.4. Layered flow pattern with solids consumption under two initial burden conditions:
(a) uniform and (b) V‐shaped [22]

33
3.3. DEM Modelling
The basis of discrete particle modelling was laid down in 1979 by Cundall and Strack who
describe a discrete numerical model to simulate discrete particles [26]. Theoretical
development and some applied modelling took place during the 1980’s [27, 28] and ever
increasing applied modelling in the 1990’s e.g. references [29‐33]. An overview of the work
done in Japan on discrete particle modelling was published by Tsuji [34]. With ever
increasing computing power the demanding DEM models became much more widely used
after 2000 in fields such as minerals processing, bulk solids handling, pharmaceuticals and
geology. The focus in the present project is restricted to the application of DEM simulation
on the ironmaking blast furnaces.

Due to the computational demand of discrete element modelling and complexity of the
processes in the blast furnace, DEM models are generally simplified or only a section of the
furnace is simulated. The work published in literature can be divided into three general
areas:

 Burden flow: simulation of the solid particle flow through the whole or part of the
blast furnace. In some cases combined with gas flow or the influence of liquids in
the hearth.
 Raceway: simulation of the raceway zone coupling a packed bed with gas
injection.
 Charging: modelling of the charging of burden materials into the blast furnace;
simulation of bins, chutes and the top layers.

These different areas of simulation all have their benefits, however, only the last one can
accurately model reality. The first two attempt to simulate more complex processes with
more variables and are simplified. In reality they include chemical reactions, high
temperatures, powder flow and softening and melting; all of which are absent in charging
simulations.

3.3.1. Burden flow

The area of the blast furnace where DEM has been most widely applied is the general
simulation of the solid burden flow. These simulations derive from the desire to understand
the development of the burden structure as it descends through the furnace as well as the
formation of characteristic phenomena such as the deadman. The focus of these models
can be on different areas of the furnace such as the top, cohesive zone or hearth.

Nouchi presented the application of the DEM method for the simulation of solid flow in the
hearth and studied how it is influenced by the liquid levels, and the results were
experimentally validated [35]. Both the experimental and DEM models consist of a box
filled with particles which are removed through an outlet; the vessel is partly filled with
water to investigate buoyancy effects on the solid flow at different liquid levels and outlet
locations. The DEM approach was able to reproduce the experimental results quite well.

34
 Figure 3.5. Simulation of solid flow in the hearth: experimental (a) and DEM results (b),
residence time; (c), particle velocity vectors [35]

A larger model is presented simulating the solids flow in the presence of liquids in the
hearth [36]. The paper presents two parts: a 10° slice of the total blast furnace and of the
hearth. In the simulations the solid flow as well as the force networks are investigated. A
novel feature is the addition of a coal abrasion model in which the coal particles are
reduced in size under pressure. Results show the formation of a force network which
supports the packed bed as illustrated in Figure 3.6. The stress in the network is much
larger than the average stress. The coke abrasion due to the force networks causes the
continuing reconfiguration of those networks. The force network was more locally
simulated around the raceway in combination with buoyancy forces due to the liquids in
the hearth.

Additional results were published where the influence of the shaft angle and tuyere length
of the furnace were investigated [37]. The geometry of this model is a 10 pie‐slice of an
industrial size blast furnace containing 20 cm diameter particles. Such a large particle
diameter is required to reduce the number of particles to what is computationally
managable. The results focus on the two clear indicators of solid flow present; the layer
structure and the deadman size and shape. Figure 3.7 shows the effect of the angle of the
furnace wall, the shaft angle, on the packed bed layer structure; particles are coloured to
indicate visualise the flow lines.

35
 (I) (II)

Figure 3.6. Force networks in: (I) blast furnace and (II) hearth where (a) base case, (b) shallow
hearth, (c) increased load and (d) moved consumption area [36]

In 2003 Zhou presented a small scale 2D blast furnace model combining solid particle flow
and gas flow [38]. The model uses DEM to simulate particle flow and computational fluid
dynamics (CFD) for gas flow. The CFD method is a modelling technique extensively used
for the simulation of liquid and gas flows. It divides the computational domain into small
cells over which a number of conservation equations are defined. Using a number of
defined boundary conditions and an iterative process the flow pattern, temperatures,
pressures and other parameters are solved for the domain.

36
 Figure 3.7. Effect of the shaft angle on solid flow; (a) 75, (b) 80 and (c) 85 [37]

DEM and CFD are combined by calculating the drag force between the particles and gas
flow with a drag force model. The simulations of solid flow are compared to experimental
results and show comparable layer structures. Gas flow is shown to strongly influence the
solid flow and size of the deadman. The boundary conditions of the front and back walls of
the geometry were found to have a large influence on the size of the deadman, most
significantly the sliding friction between the particles and wall. Further work on the 2D
model shows a decrease on deadman size with increased solid flow velocity [39]. The
model was compared to 3D simulations of a blast furnace slot geometry in which the
influence of the front and back walls on the particle flow is much lower, this results in the
formation of a smaller deadman.

37
Figure 3.8. Normal force network at different gas flow rates; (a) 1.21e‐2 kg/s, (b) 3.62e‐2 kg/s,
(c) 6.05e‐2 kg/s;
(d) probability density distribution ((a): , (b):, (c):)[40]

The small scale 3D model (1 m height) is coupled with a CFD model to investigate the
influence of gas flow [40]. For this, the DEM model is coupled with a 2D continuum gas flow
model. Results confirm what was observed in the 2D model: increased gas flow velocity
increases the deadman size. Additional study was made here of the microscopic properties
of the burden such as porosity, force‐structure and particle coordination number. From
Figure 3.8 it can be seen that an increase in gas flow rate decreases the normal force
network between the particles, shown here in red for high and blue for low forces. The
coordination number also decreases and the porosity increases.

This work was combined with an averaging method to extract time averaged particle
properties [41], allowing the discrete DEM results to be studied as approaching steady‐
state. The method is used to study the velocity field and stress distribution. It is shown that
a high solid flow rate influences the stress distribution in the burden. The authors mention

38
that in an industrial blast furnace changes in descend velocity are so small that they are not
likely to have any effect on the stress distribution.

Zhou et al. present further developments of the model in 2010 [42]. The geometry is a 4 m
high slot model with a thickness of 4 particle diameters; this is alternatingly filled with
layers of coke and ore particles with different diameters and densities. Static and rolling
friction coefficients between the particles are the same for both coke and ore. In a pre‐
defined cohesive zone the ore particles shrink and disappear, and in the raceway the coke
particles are removed. The ore layers in the cohesive zone are identified by the model and
made impermeable in the CFD gas flow model. Liquid is introduced in the hearth to
simulate the presence of molten iron and slag. It interacts with the solid particles by
buoyancy and drag force. Interaction with the gas is considered negligible and is not taken
into account.

The packed bed structure is very similar to the previous simulations; the addition of
separate particle removal zones for ore and coke result in a difference in descent velocity
above and below the cohesive zone. Coke particle removal in the raceway is mainly fed by
particles directly from above and below, but there is also a gradual renewal of coke in the
deadman zone in the centre of the hearth.

The effect of the cohesive zone shape on the gas flow is demonstrated by comparing two
very different cohesive zone shapes: V and W. As is shown in Figure 3.9(a) these result in
opposite gas flow directions through the cohesive zone. For the V shape the cohesive zone
distributes the gas inward and for the W shape outwards.

Results for the stress state in the simulation in Figure 3.9(b) show that the highest normal
forces on the particles are in the centre hearth due to the weight of the burden above. In
the raceway normal stress is low due to the removal and fast movement of particles.
The influence of liquid in the hearth was studied as shown in Figure 3.10, when
the liquid reaches a certain level the burden starts floating due to the buoyancy forces.
Because of the coke consumption in the raceway the coke‐free zone is larger along the
wall. It is also shown that an increasing coke‐free zone along the wall can form due to
increasing coke removal rates.

The authors conclude with a remark that current work is progressing on including thermo‐
chemical behaviour into the model.

39
Figure 3.9. (a) solid flow (top) and gas flow (bottom) for V (b) Normal contact forces
and W shapes cohesive zones [42]

Figure 3.10. Solid flow pattern with increasing (top) and decreasing (bottom) liquid levels [42]

40
In 2008 a paper was presented by Ueda describing a DEM blast furnace model [43]. The
3
geometry of the model is based on a large scale blast furnace, with a 5000 m inner volume
and 40 tuyeres. An 18 pie‐slice of this furnace was simulated as illustrated in Figure 3.11,
using periodic boundary conditions the particles exiting the model on one side will re‐enter
on the other. In order to make such a large scale simulation possible, the particles have
significantly increased diameters of 0.12 m for ore and 0.24 m for coke reducing the total
amount of particles in the simulation. Particles are alternatingly charged at the top of the
furnace, ore particles shrink and are removed in a pre‐defined cohesive zone and coke is
removed in three pre‐defined raceway
zones. Additional coke is removed in the
hearth at a certain interval to account for
dissolution loss of carbon in the hot metal.

Using this model the influence of a particle

hardness and rolling friction on solid flow
was investigated. The particle hardness
defined by the Young’s modulus has a large
influence on simulation time; a lower Figure 3.11. Cross‐section of DEM model [43]
hardness greatly accelerates simulations.
9 6
Figure 3.12 shows that decreasing the Young’s modulus from 2x10 Pa to 2x10 Pa does not
have any clearly visible effect on the layer structure. Rolling friction is studied here as a
method of simulating non‐spherical particle shapes while using spherical particles. When a
particle is perfectly rounded it rolls very easily, the more non‐spherical a particle becomes
the more difficult it becomes to rotate. The rolling friction is determined by the particle
weight and the rolling friction factor, a property for the resistance to rolling between two
materials. By making it artificially high it is attempted to take into account the increased
rolling resistance caused by the particle shape. Simulations using rolling friction factors of
2, 5, 10 and 20 are done; the latter two are unrealistically high and result in very steep
burden levels. With increased rolling friction the deadman size slightly increases and there
is a change in the periodicity in which the packed bed descends, however the influence
seems small.

The burden profile shows a decrease of the layer angle closer to the centre and as layers
descend the furnace. Further it is shown that the influence of the tuyeres does not extend
to more than 2 m above the raceway.

A more detailed description of the work was published including a liquid in the hearth [44].
Results of the particle hardness again show that with decreasing hardness there is a larger
fluctuation of the descent velocity. By tracing 5 random particles it is shown that
decreasing hardness results in smoother path‐lines of the particles in the lower part of the
furnace. The hardness has an influence on the normal stress distribution in the furnace:
lower hardness results in more evenly distributed stresses, with a higher hardness there are
particles with much higher stresses exerted on them.

41
 Figure 3.12. Influence of Young’s modulus on packed bed structure [43]

To visualize the influence of the rolling friction, again 5 particles are tracked while
descending in the furnace. This shows that increasing the rolling friction creates smoother
path‐lines in the lower part of the furnace. The rolling friction has no effect on the normal
stress distribution.
The residence time of a particle in a DEM simulated blast furnace is much shorter than is
the case in reality, e.g. from 6 h to 60 s, this shortens simulation time to a manageable
level. Figure 3.13 shows the influence of descent velocity acceleration of 0, 90, 180, 360 and
720 times the actual velocity in an industrial blast furnace. The total stress distribution is

42
very similar, but there is an influence on the stresses in the raceway due to the increased
particle velocity in this region. Between (a), (b) and (c) the difference is not significant and it
is concluded that the stress distribution in a blast furnace can be simulated in accelerated
conditions up to some extent.

Figure 3.13. Influence of descend velocity on stress distribution [44]

The model was applied to investigate asymmetric phenomena in a blast furnace with a
3
2000 m inner volume and 16 tuyeres [45]. The modelled geometry consists of half the
furnace with a frictionless front wall and 8 spherical raceways. Ore particles disappear
below a pre‐defined cohesive zone and coke in the raceway at a defined rate. Coke particles
have a diameter of 0.3 m and ore of 0.15 m, based on the angle of repose of real coke and
ore the rolling friction factors are 2.5 and 1.25, respectively. The descent velocity was
accelerated 80 times.

The base case in Figure 3.14(a) shows that layers descend the furnace with equal velocity
across the width of the furnace. When a new layer of particles is charged on the burden the
descent velocity peaks, most significantly when the much heavier ore is charged. Besides
these regular peaks there are also irregular slips in the raceway causing increases in descent
velocity. However this affects only the area around the raceway. Between the peaks the
blast furnace reaches a steady state condition.

43
 (a) (b)
Figure 3.14. Layer structure, path‐lines and tuyere configuration (a) base case, (b)
asymmetrical case [45]

To study the asymmetric behaviour of the burden half of the tuyeres were made inactive as
illustrated in Figure 3.14(b), resulting in coke removal from half of the modelled geometry.
The layer structure seems very much the same as the base case, however the path‐lines
show that the particles flow to the right hand side of the furnace. From the horizontal time‐
lines it can be seen that the descent velocity in the shaft is only slightly influenced, but is
higher in the lower part of the furnace on the right hand side where the tuyeres are active.
The peaks in descent velocity generated by charging are equal to the base case above the
cohesive zone; below this zone the flow is strongly influenced by the inactive tuyeres. Slips
occurring around the raceway have a more expanded area of influence on descend velocity.

Figure 3.15 shows a comparison between the normal stresses for both cases, the inactive
tuyeres now allow the burden to build up a stress network from the left side. The area of
low stress around the raceway is increased indicating higher particle velocity. The authors
introduce a slip index parameter to quantify the stability of solid flow based on the ratio of
shear stress to normal stress and friction coefficient. When the shear stress is higher than
the normal stress and contact friction, and thus the slip index is higher than one, slip occurs.
Results show more high slip index particles in the raceway zone at asymmetrical conditions
due to a decrease in normal stress. This causes more slips to occur with an increased
influence into the upper zone, decreasing the stability of the blast furnace.

44
 Figure 3.15. Normal stress distribution [45]

The influence of the cohesive zone shape on burden flow was investigated by Fan et al. in a
3
geometry based on a 5000 m blast furnace with 40 tuyeres, half of this is simulated [46].
Coke and ore have diameters of 0.4 m and 0.2 m, respectively. Properties are the same as
in the previously described article of Natsui et al. [45]. The presence of molten iron and slag
is taken into account by adding a buoyancy force in the hearth. Three different types of
cohesive zones are simulated: in the base case (A) it starts above the tuyere and has a
maximum height of 8.6 m in the centre, the second case (B) offsets the top by 5 m to the
left and in the third case (C) the cohesive zone has a height of 15 m.

The results in Figure 3.16 show that due to the buoyancy forces exerted by the liquids in the
hearth, the burden only rests on the bottom in the centre. For the biased melting or
cohesive zone (B) the layers become asymmetric. The higher cohesive zone does not
influence the layer angle and for all cases there is no influence on the burden level.

Figure 3.17 shows the normal stress distributions for the three cases. There are small forces
mostly below 0.5 MPa in the shaft of the furnace and very high forces up to 2 MPa between
the raceways in the coke zone. The biased cohesive zone (B) causes asymmetrical stress
and the elevated cohesive zone (C) results in an increased high normal stress zone. When
the stress distribution in the centre of the furnace is plotted in a graph, Figure 3.18, it can be
seen that the maximum normal stress shifts from 7 m above the tuyere level to 2.5 m and
becomes lower with a higher cohesive zone.

45
The forces exerted by the burden on the wall of the furnace gradually increase in the shaft
with a maximum in the belly. In the hearth the maximum lies approximately 2.5 m below
the tuyeres for all cases. The wall stress also shows an increase when ore is being charged
into the furnace. There is no influence of the cohesive zone shape on the normal force on
the walls.

Figure 3.16. Layer distribution with three cohesive zone shapes [46]

Figure 3.17. Normal stress distribution for three cohesive zone shapes [46]

46
 Figure 3.18. Normal stress distribution in the centre [46]

The influence of the blast furnace size and shape on the solid flow is investigated by Fan
[47]; in which cases (A) 5775 m3, (B) 3250 m3 and (C) 1444 m3 are compared. Cases (B) and
(C) are generated by reducing the diameters of case (A) with a factor of 0.75 for each
subsequent step, the height is unchanged. For each furnace the coke removal rate is
adjusted so each furnace has the same descent velocity. Ore has a diameter of 0.30 m and
coke 0.15 m.

The results of the influence of inner volume on solid flow and stress distribution are
summarised in Figure 3.19. With an increased blast furnace diameter the relative deadman
volume increases. This results in high particle velocities between the wall and deadman in
large blast furnaces. In smaller furnaces the descent velocity is more uniform in the lower
part. The stress distribution is dependent on the inner volume. In large blast furnaces
normal stress is reduced near the walls due to the higher particle velocity there. Both the
increased descent velocity and decreased normal stress cause a less stable burden descent
due to local slipping.

47
 Figure 3.19. Illustration of the influence of inner volume on solid flow
and stress distribution [47]

Figure 3.20. Solid flow pattern Figure 3.21. Simulated particle flow
[48] [48]

48
DEM solid flow modelling in the blast furnace was also undertaken by Kawai and
Takahashi, comparing results from an experimental model with 2D DEM simulation
focussing on deadman behaviour [48]. The experimental model consists of a 1 m high, 6 cm
deep blast furnace geometry filled with polymer particles. To simulate the influence of
hearth level on the deadman, the hearth of the experimental model is filled with water.
Figure 3.20 shows the solid flow pattern, plug flow in the upper part becomes funnel flow in
the lower half due to the presence of the deadman. The deadman behaviour was studied
when repeatedly increasing and lowering the water level in the hearth. Results show the
following flow patterns:

‐ Particles along the centre line and in the lower hearth display a vertical up‐
and‐down motion
‐ Particles in the deadman surface layer are pushed into the funnel flow when
the hearth level rises
‐ Particles below the outlet level move toward the outlet with a zigzag up‐
and‐down motion

A 2D DEM model was used to simulate particle flow and the deadman when repeatedly
raising and lowering the hearth liquid level. Figure 3.21 shows part of the cycle with a rising
hearth level. Particles in the deadman show the same behaviour as the experimental
model, with the addition that only the particles charged along the centre line of the furnace
attribute to the renewal of the deadman. When the hearth level rises the normal forces on
the wall increase and decrease when lowering.

The DEM model was extended to 3D and the influence of particle density on the deadman
was investigated by Kawai et al. [49]. The geometry as shown in Figure 3.22 was filled with
only alumina spheres (white) or a combination of alumina spheres and higher density glass
spheres (blue). Figure 3.22(a) shows the result of a simulation where glass particles are
charged in the centre of the furnace and alumina at the wall. The timelines in Figure 3.22 (b)
represent the alumina‐only simulation in red and two simulations containing a combination
of alumina and glass with different friction factors in white and grey. It shows that charging
heavier particles in the centre of the furnace causes horizontal flow lines and thus a smaller
deadman.

49
 (a) (b)
Figure 3.22. DEM model with timelines: (a) packed bed structure,
(b) comparison of timelines [49]

The flow of solids in the blast furnace was also studied by Mio et al. [50] [51] where
experimental results are compared to an identical DEM simulation. The validated model is
then used for a simulation of a 1:3 scale blast furnace.

Both the experimental and DEM model are 780 mm high and filled with 31098 coke
particles (9.3 mm diameter) and 450095 ore particles (3 mm diameter). Results of both
simulations show that the descent velocity of tracer particles is very similar and the
parameters used are suitable for the blast furnace simulation.

The model of the 1:3 scale blast furnace is based on a bell top charging system and
approximately 10 m high. Coke particles have a median diameter of 0.0694 m and ore of
0.025 m, the total number of particles is approximately 330000. Particles are removed by
shrinking when they reach the cohesive zone for ore and the raceway for the coke particles.
The shrinkage rate decreases inversely proportional to the surface area and is balanced
with the charging rate to ensure a constant burden level. Figure 3.23(a) shows the DEM
simulation after 2 million time‐steps; the furnace is filled mostly with the initial packing. In
Figure 3.23(b) at 10 million time‐steps the packed bed is fully replaced creating the final
layer structure. The structure seems to indicate the presence of funnel flow in the area
between the wall and cohesive zone. Particles are observed to flow with pulses as shown in
Figure 3.23(c), where the particles have a higher flow rate on the right hand side of the
furnace.

50
 (a) (a) (c)
Figure 3.23. DEM model: (a) 2M steps; (b) 10M steps;
(c) descent velocity at 8M steps [50]

A very large scale CFD‐DEM simulation of a blast furnace is presented by Yuu et al. in 2010
[52]. Using the Earth Simulator at the Japan Agency for Marine‐Earth Science and
Technology a quarter section of a full scale blast furnace was simulated. The furnace is
approximately 28 m high with a hearth diameter of 14.7 m and contains up to 16 million
particles. Coke has a mean diameter of 0.057 m and ore of 0.045 m, particles are shrunk in
the cohesive zone and are removed in the raceway. Gas is injected through 10 tuyeres with
a velocity of 250 m/s. Due to the large amount of particles the simulated time is relatively
small with 3.05 s, this represents approximately 50 minutes of the real process because of a
much faster particle removal rate of 600 ‐ 1000 times faster than reality. The model takes
into account the volume of melt created in the cohesive zone and the cohesive forces the
melt creates between particles. Figure 3.24(a) shows the simulation at the end time of
3.05 s. Below 11 m the initial particles have settled into a packed bed but above it particles
are still falling. Because the vertical velocity of the particles decreases approaching the
furnace wall on the right hand side, the layers become nearly vertical there. The raceway
size shown in Figure 3.24(b) is in good agreement with those measured experimentally,
even though the conditions are not identical. Figure 3.24(c) shows the gas flow in radial
direction through the furnace. Below the cohesive zone gas flows to the centre, but the
reduced permeability of the ore layers in the cohesive zone than changes the gas flow
direction diagonally to the wall through the coke layers.

51
 (b)

Tuyere level 7m .

(a) (c)

Figure 3.24. Large scale CFD‐DEM model: (a) layer structure at t=3.05s and cross‐sections; (b)
vertical particle velocity; (c) radial gas velocity [52]

3.3.2. Raceway modelling

The DEM method has also been used to model the blast furnace raceway, where the hot
blast is injected at high velocity into the furnace. Due to the high velocities of the gas and
its interaction with the particles the raceway is a challenging section to model.

An early combined DEM‐CFD model is presented by Xu et al., in which gas is blown into a
2D particle bed containing 10000 particles [53]. During the simulation the gas velocity was

52
steadily increased into a packed bed. Figure 3.25 shows the development of the pressure
drop over the bed. From point A to B gas velocity is increased and the pressure drop rises.
At point B a cavity or raceway starts forming, and the pressure drops as the cavity size
increases. When point C is reached bubbles start forming and the bed becomes locally
fluidised. If the gas velocity is decreased at point C the cavity collapses and the pressure
drop decreases to point F. The difference in paths between A‐C and C‐F shows a strong
hysteretic effect.

Figure 3.25. Bed pressure drop vs. gas velocity [53]

Nogami et al. compared a 3D DEM‐CFD raceway model to an experimental model to study

the effects of temperature, gas composition and top gas recycling [54]. Besides gas and
solids the simulation also contains pulverized coal injection.

53
 Figure 3.26. Raceway shapes [54]

Figure 3.26 shows a comparison between the experimental and computational models.
Taking into account the gas composition and blast temperature in the raceway, it was
confirmed that these can control raceway shape and temperature. The model was used to
study the influence of top gas recycling on the raceway. Top‐gas recycling is a way to
reduce CO2 emissions of the blast furnace by re‐injecting CO from the top‐gas. This can be
done in both the raceway and/or shaft of the furnace. The blast composition will be very
different, containing CO, higher O2 and H2, and much less N2; the blast is at room
temperature. This can significantly influence the raceway conditions.

A study of hysteresis of the raceway was reported by Singh et al. in which 2D DEM
simulations were compared to experimental results [55]. A reasonable agreement was
found between the results. The cause of the hysteresis was attributed to the change in the
inter‐particle interactions with gas flow.

Umekage et al. simulated the raceway region using a combination of DEM, Hard Sphere
Method (with a simplified contact model) and CFD; for particles, pulverized coal injection
and gas, respectively [56]. The model was applied to study the influence of the tuyere angle
and the presence of scaffolding on the furnace wall. At increasing tuyere angle the
difference between raceway depth and height become smaller, at a too high tuyere angle
flow becomes unstable. Scaffolding, such as shown in Figure 3.27, affects gas and solid
flow, the raceway shape becomes elongated and is unstable. The effects increase when the
scaffolding is closer to the tuyere.

54
 Figure 3.27. Influence of scaffolding on raceway size [56]

A comparison between DEM‐CFD model and experimental results was presented by Yuu et
al., showing fairly good agreement [57]. Figure 3.28 shows a comparison of the raceway
height (Hr), depth (Dr) and the length of the particle inflow area into the blast (Ds).

Figure 3.28. Comparison of calculated and experimental results

of raceway simulation [57]

55
3.3.3. Burden charging modelling
The Discrete Element Method is very well suited for simulation of blast furnace charging,
from the stock house to the burden, due to the relative simplicity of the physical
phenomena. Main areas of interest are particle mixing and segregation due to size and
weight, and flow behaviour when charged into the furnace. These phenomena have a large
influence on the burden distribution in the furnace (which determines gas flow) and thus on
the process stability.

In 2007 Mio et al. [58] presented DEM simulations of particle charging in the blast furnace,
focussing on hopper, chute and burden flow. Figure 3.29 shows the simulation of particles
charging in a chute to investigate particles segregation. A more detailed study combined
with experimental results was further presented [59]. Through a static chute with
adaptable angle sinter particles are discharged into a sampling box. The same is repeated
in DEM to determine the rolling coefficient and then used to study the effect of a damper
on the chute. The DEM proved to be very good in reproducing the experimental results.
Size segregation takes place with the smaller particles on the bottom of the chute and the
larger on top, the dampers reduces segregation by remixing the particles.

Figure 3.29. DEM simulation of blast furnace charging [58]

Ho et al. investigated the behaviour of the material from the chute impacting on the
burden level by experiments and DEM simulation [60]. When heavier material is charged on
a layer of lighter material it gouges into the burden level. Figure 3.30 shows the DEM results
of charging wood (red), glass (yellow) and steel (green) balls on a wood (blue) burden.

56
 Figure 3.30. Burden level after charging wood, glass and steel particles
(left to right) [60]

The effect of chute angle on particle flow behaviour and segregation was investigated by
Mio et al. [61]. A hopper and rotating chute are simulated for charging the materials on the
burden level, similar to the set‐up shown in Figure 3.29. The results indicate that the
material is centrifuged in the chute due to the rotation, shown in Figure 3.31, causing
particle segregation. When the chute angle is decreased more particles slide downwards on
the burden leading to size segregation.

Figure 3.31. Particle motion in chute [61]

An extensive simulation of the behaviour of nut‐coke (small size coke) using a bell‐type top
was presented [62]. The simulation contains a conveyor, chute, quad‐hopper and small and
large bells. This allows the segregation of nut‐coke from the sinter during the full charging
route to be studied. Figure 3.32 shows the particles flowing through the chute into the
quad‐hopper and the material in the quad‐hopper, small bell and large bell. In these figures
the nut‐coke is coloured blue and is charged first into the hopper. Two other cases were
simulated: one where the nut‐coke was charged last and one where the nut‐coke is fully
mixed with the sinter. The results show particle segregation at each transfer where nut‐

57
coke moves upwards in the sinter. After charging most of the nut‐coke ends up near the
wall due to segregation.

Size segregation in top‐bunkers was investigated by Yu and Saxén using DEM and
experimental methods [63]. Results show that segregation is influenced by a number of
particle properties. Figure 3.33 shows the effect of particle‐wall static friction on the
segregation. Also rolling friction, the amount of small particles, the filling method and the
size ratio have an influence on segregation.

Figure 3.32. Charging set‐up and size segregation [62]

58
 Figure 3.33. Influence of wall friction (top: high, bottom: low) on segregation; large (red),
intermediate (grey) and small (green) particles. [63]

3.4. Conclusions
As this chapter has shown, there has been a significant increase in the use of Discrete
Element Modelling for the ironmaking blast furnace in recent years. The literature shows
the DEM method has great potential but still has some significant challenges. Currently
quite a lot of simplifications are required to be able to do any simulation.

For DEM simulation of the blast furnace there are three fields of interest: the burden flow,
the raceway and charging. The last is well within the current capabilities of the DEM
method; there are no high temperatures, no high gas velocities and no chemical reactions.
The DEM method can be used very well to study blast furnace charging from the stock
house to the burden level. However, problems remain because of limitations on the
number of particles in the simulation due to the high computational load and because of
the characterisation of the particles. The DEM method uses spherical particles, but in the
blast furnace burden materials there is a large range of shapes and sizes. Particle properties
such as the rolling friction factor are not well known for burden materials.

Simulation of the raceway requires DEM to be coupled with a CFD model in order to
calculate the gas flow. The process in the raceway is much more complicated to model due
to the high gas velocities, very high temperatures, chemical reactions and shrinkage of the
particles. Simulations generally use some simplifications, but good results have been
achieved by modelling the raceway with DEM. Also here the simulations are limited in the
number of particles.

Large part of the work has been done on modelling the solid flow of the blast furnace.
Some simulations include gas flow and some include a liquid level in the hearth. These
simulations have focussed on the solid flow distribution and the resulting burden layer
structure, on deadman shape, on particle stress networks and on gas flow through the
packed bed. The results give a picture of how the layer structure of the burden develops

59
when descending through the furnace. It shows the location of high normal stresses in the
hearth and deadman.

All these models are significantly simplified, either by scaling up particles for large
(industrial) scale furnaces, or scaled‐down furnace size for real size particles. This
simplification lowers the computational load by reducing the number of particles. The
geometry is further reduced in size by using a slot or pie‐slice instead of the full cylindrical
blast furnace shape, again reducing the number of particles in the simulation. Besides these
geometry adaptions, the model of the process is also simplified. One aim is in order to
increase the calculation speed; the particle material’s Young’s modulus is usually reduced
allowing a larger time‐step to be used, and the descent velocity is much higher requiring a
smaller time which is needed to be simulated.

The other aim is the simplifications due to the inability of the DEM method to capture the
total process as it is taking place inside the blast furnace. Therefore the process is simplified
on the most complicated part, the cohesive zone. It is pre‐defined and no softening or
melting takes place. Further simplifications have included the absence of heat and
chemical reactions.

As can be seen in literature, with these simplifications reasonable simulations of the burden
flow can be made. But in order to really predict and understand the process in the blast
furnace the model should fully include a simulated cohesive zone and be able to determine
its shape and location. This requires the coupling of the DEM model to a CFD model as well
as reaction kinetics and softening and melting models, and is one of the main focus points
of the current project where this thesis research is part of.

60
References
1. Yagi, J.‐i., Mathematical modeling of the flow of four fluids in a packed bed. ISIJ
International 33 (1993), p. 619‐639.
2. Chen, J., Akiyama, T., Nogami, H., Yagi, J.‐i., and Takahashi, H., Modeling of solid
flow in moving beds. ISIJ International 33 (1993), p. 664‐671.
3. Austin, P.R., Nogami, H., and Yagi, J.‐i., A mathematical model for blast furnace
reaction analysis based on the four fluid model. ISIJ International 37 (1997), p. 748‐
755.
4. Austin, P.R., Nogami, H., and Yagi, J.‐i., A mathematical model of four phase
motion and heat transfer in the blast furnace. ISIJ International 37 (1997), p. 458‐
467.
5. Austin, P.R., Nogami, H., and Yagi, J.‐i., Analysis of actual blast furnace operations
and evaluation of static liquid holdup effects by the four fluid model. ISIJ
International 38 (1998), p. 246‐255.
6. Austin, P.R., Nogami, H., and Yagi, J.‐i., Computational investigation of scrap
charging to the blast furnace. ISIJ International 38 (1998), p. 697‐703.
7. Austin, P.R., Nogami, H., and Yagi, J.‐i., Prediction of blast furnace performance
with top gas recycling. ISIJ International 38 (1998), p. 239‐245.
8. De Castro, J.A., Nogami, H., and Yagi, J.‐i., Transient mathematical model of blast
furnace based on multi‐fluid concept, with application to high PCI operation. ISIJ
International 637‐646 (2000), p. 7.
9. De Castro, J.A., Nogami, H., and Yagi, J.‐i., Numerical analysis on multiple injection
of pulverized coal, prereduced iron ore and flux with oxygen enrichment to the blast
furnace. ISIJ International 41 (2001), p.
10. De Castro, J.A., Nogami, H., and Yagi, J.‐i., Three‐dimensional multiphase
mathemathical modeling of the blast furnace based on the multifluid model. ISIJ
International 42 (2002), p. 44‐52.
11. Chu, M., Nogami, H., and Yagi, J.‐i., Numerical analysis on charging carbon
composite agglomerates into blast furnace. ISIJ International 44 (2004), p.
12. Chu, M., Nogami, H., and Yagi, J.‐i., Numerical analysis on blast furnace
performance under operation with top gas recycling and carbon composite
agglomerates charging. ISIJ International 44 (2004), p. 2159‐2167.
13. Chu, M., Nogami, H., and Yagi, J.‐i., Numerical analysis on injection of hydrogen
bearing materials into blast furnace. ISIJ International 44 (2004), p. 801‐808.
14. Nogami, H., Chu, M., and Yagi, J.‐i., Multi‐dimensional transient mathemathical
simulator of blast furnace process based on multi‐fluid and kinetic theories.
Computers and Chemical Engineering 29 (2005), p. 2438‐2448.
15. Nogami, H., Chu, M., and Yagi, J.‐i., Numerical analysis on blast furnace
performance with novel feed material by multi‐dimensional simulator based on
multi fluid theory. Applied Mathematical Modelling 30 (2006), p. 1212‐1228.
16. Nogami, H., Pintowantoro, S., and Yagi, J.‐i., Numerical analysis on behavior of
unburned char and fine coke in blast furnace. ISIJ International 45 (2005), p. 1489‐
1495.
17. Zaïmi, S.A., Akiyama, T., Guillot, J.‐B., and Yagi, J.‐i., Sophisticated multi‐phase
multi‐flow modeling of the blast furnace. ISIJ International 40 (2000), p. 322‐331.
18. Zaïmi, S.A., Guillot, J.‐B., and Biausser, H., A solid flow model for blast furnace
based on hypo‐plasticity theory. Ironmaking and Steelmaking 30 (2003), p. 475‐
482.

61
19. Zhang, S.J., Yu, A.B., Zulli, P., Wright, B., and Tüzün, U., Modelling of the solids
flow in a blast furnace ISIJ International 38 (1998), p. 1311‐1319.
20. Wang, G.X., Litster, J.D., and Yu, A.B., Simulation of gas‐liquid flow in dripping
zone of blast furnace involving impermeable fused layers. ISIJ International 40
(2000), p. 627‐636.
21. Chew, S.J., Zulli, P., and Yu, A., Modelling of liquid flow in the blast furnace.
Application in a comprehensive blast furnace model. ISIJ International 41 (2001), p.
1122‐1130.
22. Zhang, S.J., Yu, A.B., Zulli, P., Wright, B., and Austin, P., Numerical simulation of
solids flow in a blast furnace. Applied Mathematical Modelling 26 (2002), p. 141‐
154.
23. Dong, X.F., Pinson, D., Zhang, S.J., Yu, A.B., and Zulli, P., Gas‐powder flow and
powder accumulation in a packed bed II: Numerical study. Powder Technology 149
(2004), p. 10‐22.
24. Dong, X.F., Pinson, D., Zhang, S.J., Yu, A.B., and Zulli, P., Gas‐powder flow in
blast furnace with different shapes of cohesive zone. Applied Mathematical
Modelling 30 (2006), p. 1293‐1309.
25. Dong, X.F., Yu, A.B., Yagi, J., and Zulli, P., Modelling of multiphase flow in a blast
furnace: recent developments and future work. ISIJ International 47 (2007), p. 1553‐
1570.
26. Cundall, P.A. and Strack, O.D.L., A discrete numerical model for granular
assemblies. Geotechnique 29 (1979), p. 47‐65.
27. Campbell, C.S. and Brennen, C.E., Computer simulations of granular shear flows.
Journal of Fluid Mechanics 151 (1985), p. 167‐188.
28. Walton, O.R. and Braun, R.L., Viscosity, granular‐temperature, and stress
calculations for shearing assemblies of inelastic, frictional discs. Journal of
Rheology 30 (1986), p. 949‐980.
29. Zhou, Y.C., Wright, B., Yang, R.Y., Xu, B.H., and Yu, A.B., Rolling friction in the
dynamic simulation of sandpile formation. Physica A 269 (1999), p. 536‐553.
30. Tsuji, Y., Tanaka, T., and Ishida, T., Lagrangian numerical simulation of plug flow of
cohesionless particles in a horizontal pipe. Powder Technology 71 (1992), p. 239‐
250.
31. Gu, Z.H., Arnold, P.C., and McLean, A.G., A simplified model for predicting the
particle flow‐rate from mass flow bins. Powder Technology 74 (1993), p. 153‐158.
32. Langston, P.A., Tuzun, U., and Heyes, D.M., Discrete element simulation of
granular flow in 2D and 3D hoppers: dependence of discharge rate and wall stress
on particle interactions. Chemical Engineering Science 50 (1995), p. 967‐987.
33. Xu, B.H. and Yu, A.B., Numerical simulation of the gas‐solid flow in a fluidized bed
by combining discrete particle method with computational fluid dynamics. Chemical
Engineering Science 52 (1997), p. 2785‐2809.
34. Tsuji, Y., Activities in discrete particle simulation in Japan. Powder Technology 113
(2000), p. 278‐286.
35. Nouchi, T., Takeda, K., and Yu, A.B., Solid flow caused by bouyancy force of heavy
liquid. ISIJ International 43 (2003), p. 187‐191.
36. Nouchi, T., Sato, T., Sato, M., Takeda, K., and Ariyama, T., Stress field and solid
flow analysis of coke packed bed in blast furnace based on DEM. ISIJ International
45 (2005), p. 1426‐1431.

62
37. Nouchi, T., Sato, M., and Takeda, K. Studies of ironmaking by discrete element
method and examples of application to analysis of phenomena in blast furnace. in
Recent progress on mathematical modeling in ironmaking. (2008), Tokyo, Japan.
38. Zhou, Z.Y., Pinson, D., Zhu, H.P., Yu, A.B., Wright, B., and Zulli, P. Discrete
particle simulation of gas‐solid flow in a blast furnace. in Third International
Conference on CFD in the Minerals and Process Industries. (2003), Melbourne,
Australia.
39. Zhou, Z., Zhu, H., Yu, A., Wright, B., Pinson, D., and Zulli, P., Discrete particle
simulation of solid flow in a model blast furnace ISIJ International 45 (2005), p.
1828‐1837.
40. Zhou, Z.Y., Zhu, H.P., Yu, A.B., Wright, B., and Zulli, P., Discrete particle
simulation of gas‐solid flow in a blast furnace. Computers and Chemical
Engineering 32 (2007), p. 1760‐1772.
41. Zhu, H.P., Zhou, Z.Y., Yu, A.B., and Zulli, P. Stress fields of the solid flow in a
model blast furnace based on CFD‐DEM simulation. in Discrete Element Methods
'07. (2007), Brisbane, Australia.
42. Zhou, Z., Zhu, H., Yu, A., and Zulli, P., Numerical investigation of the transient
multiphase flow in an ironmaking blast furnace. ISIJ International 50 (2010), p. 515‐
523.
43. Ueda, S., Nogami, H., Kano, J., and Ariyama, T. Development of blast furnace
operation simulator including discrete element concept for solid motion. in Recent
progress on mathematical modeling in ironmaking. (2008), Tokyo, Japan.
44. Ueda, S., Natsui, S., Fan, Z., Nogami, H., Soda, R., Kano, J., Inoue, R., and
Ariyama, T., Influences of physical properties of particle in discrete element method
on descending phenomena and stress distribution in blast furnace. ISIJ International
50 (2010), p. 981‐986.
45. Natsui, S., Ueda, S., Fan, Z., Andersson, N., Kano, J., Inoue, R., and Ariyama, T.,
Characteristics of solid flow and stress distribution including asymetric phenomena
in blast furnace analyzed by discrete element method. ISIJ International 50 (2010),
p. 207‐214.
46. Fan, Z., Natsui, S., Ueda, S., Yang, T., Kano, J., Inoue, R., and Ariyama, T.,
Transient behavior of burden descending and influence cohesive zone shape on solid
flow and stress distribution in blast furnace by discrete element method. ISIJ
International 50 (2010), p. 946‐953.
47. Fan, Z., Igarashi, S., Natsui, S., Ueda, S., Yang, T., Inoue, R., and Ariyama, T.,
Influence of blast furnace inner volume on solid flow and stress distribution by three
dimensional discrete element method. ISIJ International 50 (2010), p. 1406‐1412.
48. Kawai, H. and Takahashi, H., Solid behavior in shaft and deadman in a cold model
of blast furnace with floating‐sinking motion of hearth packed bed studied by
experimental and numerical DEM analyses. ISIJ International 44 (2004), p. 1140‐
1149.
49. Kawai, H. and Takahashi, H. Computer simulation of the effect of center charging
solid density on the solid flow and deadman behavior in blast furnace. in Recent
progress on mathematical modeling in ironmaking. (2008), Tokyo, Japan.
50. Mio, H., Yamamoto, K., Shimosaka, A., Shirakawa, Y., and Hidaka, J., Modeling of
solid particle flow in blast furnace considering actual operation by large‐scale
discrete element method. ISIJ International 47 (2007), p. 1745 ‐ 1752.

63
51. Mio, H., Yamamoto, K., Shimosaka, A., Shirakawa, Y., and Hidaka, J. Analysis of
granular flow in blast furnace by discrete element method. in Recent progress on
mathematical modeling in ironmaking. (2008), Tokyo, Japan.
52. Yuu, S., Umekage, T., Matsuzaki, S., Kadowaki, M., and Kunitomo, K., Large scale
simulation of coke and iron ore particle motions and air flow in actual blast furnace.
ISIJ International 50 (2010), p. 962‐971.
53. Xu, B.H., Yu, A.B., Chew, S.J., and Zulli, P., Numerical simulation of the gas‐solid
flow in a bed with lateral gas blasting. Powder Technology 109 (2000), p. 13‐26.
54. Nogami, H., Yamaoka, H., and Takatani, K., Raceway design for the innovative
blast furnace. ISIJ International 44 (2004), p. 2150‐2158.
55. Singh, V., Gupta, G.S., and Sarkar, S., Study of gas cavity size hysteresis in a
packed bed using DEM. Chemical Engineering Science 62 (2007), p. 6102‐6111.
56. Umekage, T., Kadowaki, M., and Yuu, S., Numerical simulation of effect of tuyere
angle and wall scaffolding on unsteady gas and particle flows including raceway in
blast furnace. ISIJ International 47 (2007), p. 659‐668.
57. Yuu, S., Umekage, T., and Kadowaki, M., Numerical simulation of particle and air
velocity fields in raceway in model blast furnace and comparison with experimental
data (cold model). ISIJ International 50 (2010), p. 1107‐1116.
58. Mio, H., Miyazaki, S., Seki, T., Shimosaka, A., Shirakawa, Y., Hidaka, J.,
Matsuzaki, S., and Kunitomo, K. Modeling of solid particle flow in ironmaking
process by large‐scale discrete element method. in Discrete Element Methods '07.
(2007), Brisbane, Australia.
59. Mio, H., Komatsuki, S., Akashi, M., Shimosaka, A., Shirakawa, Y., Hidaka, J.,
Kadowaki, M., and Matsuzaki, S., Validation of particle size segregation of sintered
ore during flowing through laboratory‐scale chute by discrete element method. ISIJ
International 48 (2008), p. 1696‐1703.
60. Ho, C.K., Wu, S.M., Zhu, H.P., Yu, A.B., and Tsai, S.T., Experimental and numerical
investigations of gouge formation related to blast furnace burden distribution.
Minerals Engineering 22 (2009), p. 986‐994.
61. Mio, H., Komatsuki, S., Akashi, M., Shimosaka, A., Shirakawa, Y., Hidaka, J.,
Kadowaki, M., Matsuzaki, S., and Kunitomo, K., Effect of Chute Angle on Charging
Behaviour of Sintered Ore Particles at Bell‐less Type Charging System of Blast
Furnace be Discrete Element Method. ISIJ International 49 (2009), p. 479‐486.
62. Mio, H., Komatsuki, S., Akashi, M., Shimosaka, A., Shirakawa, Y., Hidaka, J.,
Kadowaki, M., Yokoyama, H., Matsuzaki, S., and Kunitomo, K., Analysis of
traveling behavior of nut coke particles in bell‐type charging process of blast furnace
by using discrete element method. ISIJ International 50 (2010), p. 1000‐1009.
63. Yu, Y. and Saxén, H., Experimental and DEM study of segregation of ternary size
particles in a blast furnace top bunker model. Chemical Engineering Science 65
(2010), p. 5237‐5250.

64
Chapter 4
Development of a DEM‐CFD blast furnace model

4.1. Introduction
In the previous chapters it has been shown that blast furnace ironmaking is a highly
complex process in which the cohesive zone plays a very important role. It is connected to
many phenomena critical to the blast furnace process and determines process efficiency
and stability of operation to a large extent. The cohesive zone strongly influences the gas
resistance of the descending burden and distributes the ascending reductant gas to the
upper furnace part. The burden flow resistance has a direct influence on the burden descent
rate and can cause operational problems when a burden hangs and slips. Gas distribution
determines the amount of reduction taking place in the upper blast furnace and influences
the process in a number of ways.

It is very difficult to know what exactly takes place where inside the blast furnace during
production. Measurement of blast furnace conditions is limited to the in‐ and outgoing
material flows, and conditions along the side walls. Valuable information has been gained
in the past by blast furnace dig‐outs, but these show the process in a frozen state and do
not give information on changing conditions during operation.

There is a clear demand for better knowledge of the process inside the blast furnace during
production and how it reacts to changing process parameters. This demand is driven by
economic and environmental developments. The current financial crisis asks for optimal
use of costly resources and the drive for a reduction in greenhouse gases requires high
efficiency and productivity as well as new blast furnace technologies.

Mathematical modelling can produce a method to study the conditions inside the blast
furnace and how it reacts to a changes in operational parameters. The blast furnace has
been modelled extensively using continuum methods in recent decades, but these cannot
accurately simulate the discrete nature of the burden. Individual particles can be modelled
using the Discrete Element Method (DEM) which can accurately simulate solids flow.
Recent increases in computing power have seen a fast rise in its use. The research project of
which this thesis is a part aims to improve the previous continuum models by combining
the DEM model with Computational Fluid Dynamics (CFD) for gas flow simulation, and a
number of sub‐models for reaction kinetics and softening and melting.

The set‐up and results of all of the simulations done in this study are presented in this
chapter. Six different cases are presented showing the development of the model from
small test conditions to much more realistic simulations. Increasing realism was introduced
by changing geometry size and shape, by using non‐spherical particles, by including gas
temperature and by using more realistic operating parameters. Because this PhD thesis is

65
part of a larger project, work is being continued on introducing increased realism to reach
our goal of predicting the properties of the cohesive zone.

4.2. Modelling methods

All the Discrete Element Modelling was done using the EDEM software from DEM
Solutions [1]. For the Computational Fluid Dynamics part ANSYS Fluent from ANSYS was
used [2], and for the coupling between DEM and CFD the EDEM Coupling module [3].

4.2.1. CFD Modelling

Partial differential equations

Computational Fluid Dynamics (CFD) is a numerical tool to simulate fluid flow related
problems. In this work the general purpose CFD code FLUENT was used with GAMBIT as
pre‐processor [2]. A problem is solved by dividing the geometry into a number of cells and
solving a set of conservation equations for them, constrained by pre‐defined boundary
conditions. The conservation equations are a set of partial differential equations (PDE’s)
representing simultaneous conservation of mass, momentum, thermal energy and
chemical species. These govern fluid flow and related transport phenomena for fluids and
solids. Because the simulations in this thesis do not consider thermal energy or chemical
species, these PDE’s can be disregarded. The PDE’s for continuity and momentum are
given by respectively equation (4.1) and (4.2).

∙ (4.1)

where  = density, = velocity factor, t = time and Sm = mass source. The continuity
equation governs mass conservation. The first term on the left hand side represents the
rate of increase of mass, the second represents rate of flow of mass into a fluid element.
The right hand term represents the mass added to the continuous phase from a dispersed
second phase (e.g. due to vaporization) and any user‐defined sources.

∙ ∙ ̿ (4.2)

In which the stress tensor ̿ is defined as:

2
̿ ∙ (4.3)
3

where p= static pressure, = gravitational acceleration, = external body forces,  =

molecular viscosity, I = unit tensor. The Navier‐Stokes momentum equation satisfies
Newton’s second law; the rate of change of momentum equals the sum of forces acting on

66
the element. The terms on the left hand of equation (4.2) respectively refer to rate of
momentum increase of the element and the net rate of momentum into the element. On
the right hand the terms refer to the surface force on the element due to viscous stress, due
to the pressure gradient and due to gravity and other body forces.

An extra set of equations is required to close the system of fluid dynamic equations; they
also supply relations between thermodynamic variables as well as between transport and
thermodynamic properties. These equations are called equations of state, an example of
which is the ideal gas law.

Numerical schemes
FLUENT can use two numerical schemes to solve the PDE’s, the segregated or the coupled
solver shown in Figure 4.1 [2]. Both are control volume based methods that consist of:

 Division of the domain into discrete control volumes using a computational grid
 Integration of the governing PDE’s on the control volumes to construct algebraic
equations
 Linearization of the discretized equations and solution of the resultant linear
equation system to yield updated values of the dependent variables

Segregated solver
The segregated solver solves the governing equations sequentially, i.e. segregated from
one another. To obtain a converged solution many iterations are necessary. Each iteration
consists of the following steps:

1. Fluid properties are updated based on the current solution (for the first iteration
initialized values are used)
2. The u, v and w momentum equations are each solved in turn using current values
of the pressure and face mass fluxes in order to update the velocity field
3. The velocities obtained in step 2 may not satisfy the continuity equation for the
volume, and an equation for a pressure correction is derived from the continuity
equation and the linearized momentum equations. This pressure correction
equation is then solved to obtain the necessary corrections to the pressure and
velocity fields and face mass fluxes to satisfy continuity.
4. The equations for the scalars such as turbulence, energy, species, radiation and
user‐defined scalars are solved using the previously updated values for the other
variables
5. A check for convergence of the equation set is made

These steps are repeated until convergence is reached.

67
Coupled solver
The governing equations for continuity, momentum, energy and species transport are
solved simultaneously, in other words coupled. Governing equations for additional scalars
will be solved sequentially using the following procedure:

1. Fluid properties are updated based on the current solution (for the first iteration
initialized values are used)
2. The continuity, momentum, energy and species equations are solved
simultaneously
3. Equations for additional scalars such as turbulence, radiation or user defined are
solved using the previously updated values of the other variables
4. The convergence is checked, the steps are repeated until convergence is reached

Update properties

Solve momentum equations Update properties

Solve pressure correction Solve continuity, momentum,

(Continuity equation) energy and species equations
Update pressure, face mass flow rate simultaneously

Solve energy, species, turbulence Solve turbulence and other

and other scalar equations scalar equations

Converged? Stop Converged? Stop

Figure 4.1. Segregated and coupled solution methods [2]

Linearization: implicit versus explicit

The discrete non‐linear governing equations are linearized to produce a system of
equations for the dependent variables in every computational cell [2]. This linear system is
solved to obtain an updated flow field solution.

The governing equations can be linearized in two ways, “implicit” and “explicit”:

 Implicit: The unknown variable in a cell is calculated using a relation which

includes existing and unknown values from neighbouring cells. Each unknown
appears in multiple equations, and these equations are solved simultaneously.
 Explicit: Each unknown variable is calculated using a relation that includes only
existing values. Each unknown will appear in only one equation and the equations
are solved one at a time.

68
In the segregated solution method implicit linearization is used. The segregated approach
solves for a single variable field (e.g., pressure) by considering all cells at the same time. It
then solves the next variable field, again considering all cells at the same time, and so on.
There is no explicit option for the segregated solver.

In the coupled solution method implicit or segregated linearization of the governing

equations are available. These, however, only refer to the simultaneously solved equations.
The additional variables such as turbulence, radiation or a user defined scalar are linearized
and solved implicitly using the segregated solving method. The coupled implicit approach
solves all variables in all cells at the same time. A point implicit (block Gauss‐Seidel) linear
equation solver is used in combination with the Algebraic Multigrid method (AMG). This is a
method to accelerate convergence by solving the problem on a coarser grid or resolution
and interpolating this data on a finer grid. The coupled explicit approach solves for all
variables one cell at a time, using a multi stage (Runge‐Kutta) solver and the AMG method.

Discretisation
To solve the PDE’s, first the flow domain is discretised by generating a large amount of cells
or control volumes. FLUENT then uses a control‐volume‐based technique (finite volume
method, FVM) to discretise the governing PDE’s to algebraic equations that can be solved
numerically for each cell [2]. This technique consists of integrating the governing equations
about each control volume, giving discrete equations that conserve each quantity on a
control‐volume basis. Using this technique the integral form for an arbitrary control volume
V can be written as:

∅ ∙ Γ ∅∙ ∅ (4.4)

where:
ρ density
 scalar quantity
velocity vector
surface area factor
Γ diffusion coefficient for 
∅ gradient of 
SΦ source of  per unit volume

Equation (4.4) is applied to each control volume, or cell, in the computational domain. The
discretisation of equation (4.4) on a given cells is:

69
 ∅ ∙ Γ ∅ ∙ ∅ (4.5)

where:
Nfaces number of faces enclosing cell
∅ value of  convected through face f
∅ ∙ mass flux through the face
area of face f
∅ magnitude of   normal to face f
V cell volume

The equations solved by FLUENT are similar to equation (4.5) and can be used for multi‐
dimensional, unstructured meshes composed of polyhedra.

First‐ and second‐order upwind

FLUENT stores discrete values of the scalar at the cell centres, however, face values for øf
are required in equation (4.5). These values therefore have to be extrapolated from the
values stored at the cell centres. This can be done using an upwind scheme; upwinding
means that the face value is derived from quantities in the cell upstream relative to the
direction of the normal velocity in equation (4.5). There are several upwind schemes
available in FLUENT: first‐order upwind, second‐order upwind, power law and QUICK. In
this thesis work only first‐ and second‐order upwind schemes were used.

When the first‐order upwind scheme is used quantities at the cell faces are determined by
assuming that the cell centre values of any field variable represent a cell average value and
hold throughout the entire cell; the face quantities are identical to the cell quantities. The
second‐order upwind scheme cell face quantities are computed using a multi‐dimensional
linear reconstruction method. Higher order accuracy is obtained at cell faces through a
Taylor series expansion of the cell centred solution about the cell centroid. The face value is
computed by:

∅ ∅ ∅∙∆ (4.6)
where

∅ ∑ ∅ (4.7)


where  and  are the cell centred value and its gradient in the upstream cell, and s is
the displacement vector from the upstream cell centroid to the face centroid. The gradient

70
is calculated using the divergence theorem. Face values f are computed by averaging 
from the two adjacent cells.

Under‐relaxation
Because of the non‐linearity of the equation set solved by FLUENT, it is necessary to
control the change of ø. Under‐relaxation reduces the change of ø during each iteration,
and can be seen as:

∅ ∅ ∆∅ (4.8)

where:
∅ old value of variable ø
∆∅ computed change in ø
α under‐relaxation factor

4.2.2. DEM Modelling

In the Discrete Element Method every individual particle is tracked and its motion is
calculated based on Newton’s second law of motion and the governing equation for the
translational motion can be written as:

F F F (4.9)

where particle i has mass mi and velocity vi. The right hand side contains terms for contact,
gravity and drag forces. This general governing equation for the solid motion is solved with
the general purpose DEM software package EDEM [1]. In equation (4.9), the collision forces
are calculated using the Herz‐Mindlin No Slip contact model based on the work of Mindlin
and Deresiewicz [4]. The two forces governing the contact model are the normal force and
the damping force. The former is a function of the equivalent Young’s modulus E* (stress‐
strain relation) according to equation (4.10) where R* and δn are the equivalent radius and
normal overlap. Damping force shown in equation (4.11) is a function of the particle
rel
properties; equivalent mass m* and the normal component of the relative velocity vn ; and
material properties; the normal stiffness Sn and the coefficient of restitution e.

4 ∗ ∗ ⁄
F (4.10)
3

F 2 56 ∗ (4.11)

71
with
ln
√ln
and
∗ ∗
2

Besides translational motion as governed by equation (4.9) particles also undergo

rotational motion, which is governed by equation (4.12).

T M (4.12)

where Ii is the moment of inertia, ωi the angular velocity, Ti the torque generated by the
tangential forces and Mi is the torque generated by the rolling friction. The tangential force
torque given in equation (4.13) depends on two components: the tangential force Ft and the
tangential damping force Ftd. The former is shown in equation (4.14), where δt is the
tangential overlap and St the tangential stiffness; and the latter in equation (4.15) where
rel
vt is the relative tangential velocity.

T R F F (4.13)
F (4.14)
with
∗ ∗
8

F 2 56 ∗ (4.15)

Coulomb friction, μsFn with μs the coefficient of static friction, limits the tangential force.
Rolling friction is included in the equations by applying the negative torque shown in
equation (4.16), where μr is the coefficient of rolling friction, Ri the distance from the centre
of mass to the contact point and ωi the unit angular velocity vector at the contact point.

M (4.16)
For an extensive background on the theory of discrete particle modelling we would like to
refer to Zhu et al. [5].

72
4.2.3. DEM – CFD Coupling
The EDEM – Fluent Coupling Module is used to couple the DEM simulation with CFD, and
uses the existing Eulerian – Eulerian multiphase model in Fluent. Equations (4.17) and (4.18)
show the continuity and momentum equations for the gas phase, where ε is the volume
fraction, ρ the density, u the velocity, μ the viscosity, p the pressure and S the momentum
sink. The momentum equation is based on the Model B as proposed by Gidaspow [6] where
the pressure drop is only in the gas phase and is not shared by the solid phase as is
described by Model A. For mono‐sized particles there is little difference between both
models, for the fluidization of binary mixtures Model B is preferred as shown by Feng [7].
Although in this simulation two particle sizes are used they are separated in mono‐sized
layers and even if mixed, Model B is preferred.

∙ u 0 (4.17)

u
∙ uu ∙ u g S (4.18)

The solid volume fraction is copied from the DEM to the CFD model; the coupling module
over‐rides the continuity equation for the solid phase such that it is not solved by FLUENT.
To determine the solid fraction in a cell, regular sample points are taken in the bounding
box of each particle. A point is kept if the point lies in a particle and it is then checked in
which CFD cell it lies. When this is done for all particles the solid fraction for each cell can be
determined by:

1 (4.19)

In which nc is the number of sample points of particle p in the mesh cell, N is the total
number of sample points for a particle and Vp the particle volume.

The momentum coupling causes an additional force on the DEM particles based on the
local drag force. In the CFD simulation a momentum sink is added to each of the mesh cells
to represent the effect of the momentum transfer to the DEM particles. Consider the
momentum sink, S, on a mesh cell:

∑ ∑
(4.20)

where F is the force on a particle in a particular iteration from the fluid. The sum is over the
number of DEM iterations carried out between CFD iterations which generally have a larger
time step than the DEM simulation. The drag force on the individual particles is calculated
using the Di Felice [8] drag model as shown in Equation (4.21).

73
 (4.21)

where:
0.5 v v v v

4.8
0.63 .

where ε is voidage and χ is given by:

1.5 log
3.7 0.65 (4.22)
2

74
4.3. Modelling approach
The work described here is based on several different geometries to describe the
ironmaking blast furnace. A full scale industrial blast furnace has a height of approximately
35 m and is up to 15 m in diameter. Due to high computational requirements for DEM
simulations, it is as yet impossible to model the full scale furnace with its millions of
particles. Therefore, a reduced scale geometry has to be chosen in order to decrease the
total amount of particles to be simulated. During the course of this work two main options
have been used: the slot model and the pie‐slice model. The former is a cross‐section of the
furnace with parallel front and back planes and is more or less a semi‐2D model as it does
not truly capture the three dimensional cylindrical blast furnace shape. The pie‐slice model
does and its front and back planes converge on the centre‐line of the furnace.

This chapter is divided into several sections each describing the particulars of a case; the
geometry, boundary conditions, results and discussion. Six cases are described:

 Lab scale slot model:

Small scale slot model with single spherical particle type, to compare with
experimental and simulated literature data.
 Pie‐slice model of scaled‐down blast furnace:
Pie‐slice model of a scaled down industrial blast furnace, simulations with a
single particle type of spherical or non‐spherical particles. To investigate
influence of particle shape on solid flow.
 Blast furnace models with different geometry:
Study on the influence of the geometry shape; pie‐slice, slot model and full
cylindrical furnace.
 Slot model of Experimental blast furnace:
Slot model of an experimental blast furnace with both spherical and non‐
spherical particle layers. Includes gas flow to investigate its influence on solid
flow.
 Modelling the burden charging and wall conditions:
Small supplementary study on the charging of the burden in the blast furnace

4.3.1. Lab scale slot model

The first simulations were performed to assess the performance of the software used in this
research with the results from the work by Zhou et al. [9] and to study the influence of non‐
spherical particles. In the experiments from the literature an approximately 250 mm high
blast furnace scale model was filled with glass beads. In the model there is no hearth
present and the particles are removed from the bottom while new particles are charged in
layers at the top. This is compared with a DEM model with the same scale and geometry. A
second simulation in the study of Zhou compares the results to a DEM model of a 1 m high
geometry including a hearth, where the particles are removed from the raceway zones.

75
These literature results were compared to two DEM simulations. The first has a very similar
geometry as the experimental model, though at the larger scale of the DEM simulation
from literature, which we shall call Case Ia. The second adds a hearth and non‐spherical
particles to the simulation to study their effect on particle flow, Case Ib. Our results are
compared to both the experimental results and the DEM simulation presented in the
journal paper of Zhou et al .

Geometry and particles

The geometries in this simulation are based on the experimental and simulation set‐ups
from the literature [1]. Two geometries are presented in which the particles at the start of
the simulation are coloured grey and the ones with which the model is refilled after the
start are red and white. The sole purpose of these colours is to create distinct layers
allowing the flow pattern to be studied; all particles are physically exactly equal. Both the
front and back planes of the model are made of Perspex to make the particle flow visible in
the experiments. Particles are discharged by gravity through two gaps in the bottom next
to side‐walls of the furnace. The first literature DEM simulation uses the same geometry as
is used in the experiment, and the second literature DEM simulation is four times larger in
height and width. These literature results are compared to a DEM simulation of a geometry
with an equal size to the larger DEM simulation. This geometry is shown in Figure 4.1; it has
a height of 0.925 m, width at the top is 0.275 m and the depth is 0.05 m. All particles are
spherical and have a diameter of 0.01 m, particle properties and model parameters are
shown in Table 4.1 and are determined based on the values from literature.

The second geometry used here is similar to the first but extends below the raceway zones
to contain a hearth zone. In the literature [9] particles are removed in two rectangular
zones representing the raceway zones rather than through the bottom for comparison with
the experimental results. Due to software restrictions particles in our model are removed
by gravity through 0.03 m high outlets in the side walls. Case Ib is a slot model 1.025 m
high, 0.275 m wide and 0.050 m thick. Expanding the literature research besides
comparison of the results, the influence of non‐spherical particles has also been
investigated. A simulation using only 0.01 m spherical particles has been compared to a
case including non‐spherical particles. The former represents only pellets and the latter
coke and pellets in the actual blast furnace process. By combining 8 spheres of 0.01 m
diameter in a cubic arrangement with 0.004 m overlap non‐spherical particles are created
with a height, width and depth of 0.016 m, shown in Figure 4.1(b).

76
 0.275 m

0.925 m
Ia

0.430 m

Ib 0.1 m

(a) (b)

Figure 4.1. (a) DEM model geometry of the present study, case Ia and Ib
with hearth extension; (b) particle shape

Table 4.1. Simulation parameters

Variable Case Ia Case Ib
Particle diameter 0.01 m 0.01 – 0.016 mm
Particle density 2500 kg/m3 2500 kg/m3
Poisson’s ratio 0.172 0.172
Shear modulus 3×108 Pa 3×108 Pa
Coefficient of restitution 0.0001 0.6
Coefficient of static friction (part.‐part.) 0.4 0.4
Coefficient of static friction (part.‐wall) 0.4 0.4
Coefficient of rolling friction 0.1 0.1
Time step 2.5×10‐5 s 3×10‐5 s

Boundary conditions and simulation procedure

Both front and back walls in the models are physical walls, representing the Perspex walls
used in the experiments of the literature [1]. The simulation parameters are taken from the
literature or are very similar to it due to a difference in the elastic constants used in the
model, in EDEM described with the Shear modulus and Poisson’s ratio, but in the literature
model with Young’s modulus. To be able to run the simulations within an acceptable time

77
the shear modulus, which determines for a large part the size of the time step, was
reduced. This is an unavoidable and general practice for larger size DEM simulations and
simulations show it does not have a significant influence on the packed bed flow [10]. In
case Ia the coefficient of restitution, which is the ratio of the velocity before and after
particles collision, is kept at a minimal value to ensure the stability of the simulation with
which there were some issues. At higher values the simulations would regularly become
unstable with particles generating very high velocities within a time‐step.

Each simulation starts by creating random particles in the geometry, and they settle into a
packed bed which is the initial packing for the simulation. In the first picture in Figure
4.2(a), which shows the results of case Ia, the initial particles are coloured grey. When the
furnace is filled up to approximately 100 mm under the top, the two outlets in the bottom
are opened and the particles flow out by gravity. When the surface of the packed bed drops
down to the level where the vertical throat ends and the sloping shaft starts , a new layer of
1000 particles is created. These layers are alternatingly coloured red and white to visualize
the emerging flow patterns.
In case Ib the initial packing consists of 17300 spherical beige particles. After
filling, the outlets in the sidewalls are opened and the particles flow out under gravity. As
the packed bed level drops layers of grey and green particles are created. For the
simulations containing only spherical particles, all particles are the same in their physical
properties. When non‐spherical particles are included in the simulation these are coloured
grey and alternate with green layers of spherical particles. Each spherical layer contains 750
particles and each non‐spherical layer 150. The simulation is stopped when 1230 spherical
particles are removed or the equivalent volume in spherical and non‐spherical particles.

Results and discussion

The results of case Ia and the results from the literature are presented in Figure 4.2, where
the first row shows the results from our DEM simulations, the second row the experimental
results from the literature and the third row the DEM results from the literature. Due to the
uncontrolled outflow of particles the simulations cannot be compared based on simulation
time. From the results in Figure 4.2(a) three different flow regions can be observed in the
model. In the top part of the furnace there is plug flow where particles have the same speed
across the width of the furnace. In the stack the flow changes and the layer structure starts
to resemble a W‐shape caused by slow flow in the centre of the furnace, rapid flow to the
outlets at ¼ and ¾ of the width and again slow flow close to the walls. The W‐shape of the
layer is generally associated with materials flow in the blast furnace. The third picture in
Figure 4.2(a) shows the formation of a low velocity zone in the lower area of the furnace. In
the real blast furnace this zone is present as the Deadman which is a stagnant zone in the
lower part of the furnace where the coke has a much longer residence time than the
particles in other parts of the furnace. In the last picture the zone becomes invisible due to
the gradual replacement of grey particles by red and white ones, though the outline can
still be seen in these layers as well.

78
 t = 5.5 s t = 17.6 s t = 23.3 s t = 40.8 s

(a) DEM Simulation from this work

(b) DEM Simulations, literature [9]

(c) Experimental results, literature [9]

Figure 4.2. Comparison experimental results and DEM simulations

79
Most of the results are in good agreement with the literature; the first three snapshots are
very similar. The fourth snapshot shows a different flow pattern in the lower area. Both the
literature experiment and the literature DEM model show particle mixing directly above the
Deadman. This is caused by layer break‐up in the literature model, the slightly thicker
layers in our model do not break up and therefore do not show this particle mixing. This
mixing is also concealed due to the uncontrolled outflow in our model, the rate of which
can only be partially influenced, e.g. reducing the outlet size to slow down the flow causes
hanging and all flow is stopped. Particle mixing does occur when the simulation is
continued a bit longer, however, at that point the Deadman has disappeared. The smaller
Deadman is the main difference between the literature and our simulation. Downward
pressure from the particles in the centre push the particles outwards to the outlets.

Due to the small and rapidly disappearing Deadman the geometry was changed to include
the hearth below the raceway level. The presence of particles between the raceways
creates an uneven level compared to the smooth bottom from Case Ia; this was expected to
assist in the creation of a Deadman extending above the raceways. A similar geometry was
also studied in the literature [9]. In addition to the changed geometry there are non‐
spherical particles introduced to the simulation to study the shape effect on the particle
flow.

The results presented in Figure 4.4(a) show similar results to the previous model in Figure
4.2(a). There is plug flow in the upper part, a W‐shaped layer structure in the lower zone
and the presence of a small Deadman extending above the raceway level. Comparison with
the results from literature in Figure 4.4(c) again shows a smaller Deadman and much less
particle mixing. In this case the ways in which the particles are removed are rather different
in our work compared to the literature. This is likely to have a significant effect on the flow
pattern, but software restrictions at this stage did not allow for a similar method of particle
removal. The literature DEM results seem to differ from the experimental results: the layers
in the simulation curve downwards at the walls where the experimental results in Figure
4.2(c) clearly show upwardly curving layers also seen in our DEM results presented in Figure
4.4(a) and (b).
The simulation was also performed using a reduced particle‐wall static friction
coefficient [11], shown in Figure 4.4(d) and is much more similar to our simulation results in
Figure 4.3(a). For this simulation the static friction coefficient was reduced from 0.4 to 0.1.
Layer mixing is still stronger, but the layer structure has become much more horizontal.
Most importantly, the Deadman is significantly smaller and much more in line to our
results.

Presence of non‐spherical particles has a large influence on the particle flow as can be seen
in Figure 4.4(b). Layer curvature at the walls and the W‐shape of the layers become much
less pronounced. The particles become more evenly distributed along the width of the
furnace. This can also be seen in the reduced size of the Deadman due to the increased flow
through the centre pushing it down. Cause of these effects is the increased resistance to

80
particle shear. Non‐spherical particles are interlocking in a larger degree, making it harder
for a cluster of particles (in this case a vertical band) to move past another cluster. A close‐
up view is shown in Figure 4.3, (a) and (b) of the mixed simulation of Figure 4.4(b), and
4.3(c) of the spherical simulation in Figure 4.4(a). From the angular velocity in Figure 4.3(a)
can be seen that the spherical particles rotate considerably more than the non‐spherical
particles. The non‐spherical particles resist rotation due to their shape. This causes a higher
torque to be applied to the non‐spherical particles compared to the spherical particles, as
shown in Figure 4.3(b). In the simulation where only spherical particles are present the
torque is also significantly lower. This causes the non‐spherical particle layers to be more
resistant to deformation than those consisting of only spherical particles.

(a) (b) (c)

Figure 4.3. (a) Angular velocity, scale blue to red: 0 – 15 rad/s,
(b) and (c) Torque, scale blue to red: 0 – 0.008 Nm

The results from literature have not been entirely comparable to our simulations. General
flow patterns are similar but the Deadman is consistently smaller than those present in the
literature experiments and simulations. These dissimilarities are likely in part due to the
differences in geometry and models. But the low‐friction results from Figure 4.3(d) show
that the front and back walls have a very large influence on the particle flow. This
demonstrates that these small models are very strongly influenced by boundary conditions
and simulation parameters. The use of non‐spherical particles has also proven to have a
large effect on particle flow and should be included in simulations of the blast furnace.
The experimental models used to validate the DEM simulations represent small
slots filled with plastic or glass beads, they do not represent a blast furnace. The validity of
the results for the DEM method when used on much larger problems remains an issue.
However, several papers have described good results comparing experimental to DEM
models [12‐19]. This sufficiently proves that the DEM method can accurately simulate
particle flow. It is the goal of this project to simulate a large scale blast furnace, making it

81
infeasible to do validation experiments on a meaningful scale. Therefore, no lab‐scale
experiments have been done as part of this work.

(a) (b)

(c) (d)
Figure 4.4. Slot DEM model simulation results: (a) spherical (beige, green, grey),
(b) spherical (beige, green) and non‐spherical (grey),
(c) spherical (blue, green, red) literature [9],
(d) spherical (blue, green, red) low particle‐wall friction [11]

82
4.3.2. Pie‐slice model of scaled‐down blast furnace
With the previous cases experience was gained with the DEM modelling technique and our
method was compared to the literature results. The two main points from the slot model
(Case I) were a smaller Deadman than what was to be expected from the literature results
and the influence of non‐spherical particles on the flow pattern. The large degree of
simplification and abstraction of our model compared to the actual blast furnace as well as
the limited availability of data for validation, have made us decide to focus on comparative
parameter studies rather than direct simulation of real problems. Care should therefore be
taken when drawing conclusions from these results on the actual blast furnace.

In this case the geometry is changed to a small size blast furnace rather than based on an
experimental set‐up. The solid burden in the blast furnace can contain pellets, sinter, lump
ore and coke. These have very different sizes and shapes ranging from small spherical
pellets to large blocky coke. Because the particle shape has proven to influence the particle
flow in this case several different non‐spherical particle types are compared.

Geometry and particles

A pie‐slice model (Case II) was created based on a scaled down industrial blast furnace.
Compared to the slot model from the previous case the pie‐slice model is a better
representation of the 3D blast furnace shape. The surface area of the horizontal cross‐
section of the pie‐slice model has cubic growth when the furnace diameter widens, but with
the slot model the cross‐section only increases linearly. The pie‐slice model has a height of
3.5 m and a maximum radius of 0.77 m, giving the total blast furnace a width of 1.54 m. The
angle of the pie‐slice at the centre‐line is 15 degrees, the sharp edged zone created here is
removed from the geometry because the depth is smaller than the particle diameter and
thus it would be open to gas flow if it were included. As in Case I the particles are removed
by gravity through an outlet in the side wall, and thus the size of the outlets determines the
removal rate of the particles.

Four different particle types are used in the simulations: one spherical and three non‐
spherical. The latter are a square bi‐pyramidal formation consisting of 6 spheres, a cubical
formation of 8 spheres and an elongated particle of three spheres. The particles all have a
maximum diameter of approximately 0.022 m and only a single type was used for each
simulation. The particles and the pie‐slice geometry are shown in Figure 4.5.

Boundary conditions and simulation procedure

The simplified geometry used in the simulation contains front and back (physical) walls not
present in the real blast furnace. To minimize the influence of the walls on the particle flow
specific wall parameters were chosen to mimic realistic conditions. In reality the bed
extends forwards and backwards and at the location of the front and back planes the
particle packing is unrealistically disturbed by the presence of the walls. Ideally, periodic
boundary conditions should be used where a particle exiting from the back plane will re‐
enter through the front plane. The software used, however, did not allow periodic

83
boundary conditions at that stage. Therefore, in order to realistically simulate the
continuation of the packing, the static wall friction was set at a very low value to account
for the absence of friction when the bed is moving down at an equal rate at both sides of
the boundary. In reality the particles are interlocking causing a large resistance to rolling as
the particles would have to roll over each other. To prevent the particles rolling along the
wall a high wall rolling friction was chosen. Boundary conditions for the simulations are
shown in Table 4.2, where the last two values for static and rolling friction are for the front
and back walls.

Spherical Non‐spherical 1

15°

3.5 m
1 particle, 22 mm 6 particles, 6 x 13 mm
square bi‐pyramidal
Non‐spherical 2 Non‐spherical 3

0.77 m

8 particles, 8 x 14 mm 3 particles, 3 x 16 mm
cubical elongated
Figure 4.5. Particle and geometry shapes

84
 Table 4.2. Particle properties
Variable Pie‐slice model
Particle size ±0.022 m
Particle density 1000 kg/m3
Poisson’s ratio 0.5
Shear modulus 1×107 Pa
Time step 5×10‐5 s
Coefficient of restitution 0.5
Coefficient of static friction (particle‐particle) 0.6
Coefficient of static friction (particle‐wall) 0.4
Coefficient of rolling friction 0.05
Coefficient of static friction (particle‐inner wall) 0.001
Coefficient of rolling friction (particle‐inner wall) 0.5

The shear modulus is decreased to speed up the simulations to workable levels. Density
and Poisson’s ratio are default values as the simulation focuses on comparing the influence
of particle shapes, not attempting to accurately simulate the real process. A fifth case was
done where the front and back walls are given the same friction coefficients as between the
particles themselves to simulate a real wall in order to compare with the previous
simulations.

Each simulation is started by filling the geometry with particles for the initial packing. In
each simulation only a single particle type is present both as initial packing as in the
charged layers during the simulation. When the furnace is filled particle removal starts from
the outlet in the sidewall. The initial packing consists of layers with equal thickness as
shown in Figure 4.6(a). As the burden level at the top descends new layers are added each
containing 1500 particles. These layers will have slightly different layer heights for each
particle type as the particle volumes and resulting bed porosities differ. The use of coloured
particles in layers is to create distinct layers allowing the flow pattern to be studied The
simulation is stopped when 1000 magenta particles have been removed from the furnace
and all initial layers have been removed. When the simulation would continue the layer
structure could change slightly due to flow patterns evolving which take longer to develop
than the current simulation time. This could also result in a slightly different deadman, but
this would be increasingly hard to see due to the increased particle mixing.

Results and discussion

The results from the pie‐slice simulations are shown in Figure 4.6; in (a) the initial packing
of equal height spherical particle layers is shown before the particle removal is started; this
initial packing is similar for the other cases but with the different non‐spherical particle
types. The spherical case (b) is compared to the three non‐spherical cases (c), (d) and (e).
The latter are filled with non‐spherical particles only and the particle colours serve to

85
visualize the flow pattern. Figure 4.6(f) shows the results when front and back walls are real
walls as was done in the previous cases.

(a) (b) (c) (d) (e) (f)

Figure 4.6. Pie‐slice model: (a) initial packing; and simulation results: (b) spherical,
(c) non‐spherical 1, (d) non‐spherical 2, (e) elongated, (f) physical inner walls

The pie‐slice model used in Case II for scaled down blast furnace shows similar results to
the slot model of Case I when looking at the overall layer structure, though the cases use
very different wall conditions. The simulation of the spherical particles in Figure 4.6(b)
shows strongly upward curving layers at the side wall, a distinct W‐shape cannot be seen
due to the much smaller Deadman. The non‐spherical particles of the first type show a
significantly reduced but still present layer curvature, but the non‐spherical 2 and 3
particles hardly show any layer curvature. These results are comparable to the ones from
the slot model and show again that the non‐spherical particles create a more even flow
across the width of the furnace.

86
A noticeable difference can be seen in the Deadman size which grows with increasingly
non‐spherical particles. This seems to contradict the results from the slot model where the
contrary happens. It can be explained by the fact that the Deadman in Case I is composed
of the spherical particles from the initial packing and in the pie‐slice model only a single
particle shape is used in each simulation. A Deadman formed by non‐spherical particles is
more stable than a spherical one due to its resistance to rolling and thereby resists the
increased pressure in the centre caused by the more even flow across the furnace width.
The more the particle shape deviates from a perfect sphere the larger the Deadman
becomes.

When the simulations containing only spherical particles are compared it can be seen that
the Deadman in Figure 4.3(b) is considerably larger than the one formed in Figure 4.6(b).
This can be attributed to the front and back wall friction coefficient, which for the slot
model simulation represent actual walls and for the pie‐slice are defined to represent the
conditions inside the blast furnace as mentioned in Section 4.3.2. Deadman formation in
slot models with physical walls can be attributed to the presence of these planes rather
than particle or solid flow properties. If the simulations of the pie‐slice model are made with
physical walls with friction coefficients equal to the particle‐particle coefficients the results
change drastically, as can be seen in Figure 4.6(f). It should be mentioned that the effects of
the front and back planes are strongly magnified by the decreasing thickness of the
geometry near the centre. It is, however, clear that even taking this into account the
influence of the planes on the formation of the Deadman is very large.

In cold solid flow models with a simplified geometry the presence of the Deadman can
largely be attributed to the front and back plane boundary conditions rather than particle
flow characteristics. The geometry used should be carefully chosen as it has a very large
influence on the solid flow. A small increase in the size of the Deadman can be seen if it is
build‐up from non‐spherical rather than spherical particles.

The presence of non‐spherical particles creates a more evenly distributed flow across the
width of the geometry. This is caused by the increased resistance to shearing of the packed
bed due to the reduced freedom of particles to rotate in the bed.

4.3.3. Blast furnace models with different geometries

The previous simulations have shown a very large effect of the shape and boundary
conditions of the geometry on the solid flow pattern. A very important consideration in the
choice of geometry is the amount of particles present in the simulation. An increasing
number of particles quickly decreases computation speed. The geometry influence clearly
shows the need for a more detailed investigation, specifically on the formation of a
Deadman in cold single phase solid flow. Five simulations were done to compare the main
geometry choices available to find the optimal shape. The simulations investigate three
aspects: slot model thickness, slot model width and the pie‐slice geometry.

87
Geometry and particles
The first geometry is a half diameter, 0.77 m, slot model of 3.5 m high filled with 0.022 m
spherical particles. Thickness is 0.4 m and the geometry is filled with 9000 particles. This
case is compared to two other slot models: one with double thickness and one with the full
diameter of 1.54 m, both of these contain 18000 particles. The fourth geometry is a pie‐
slice model, as was used in previous simulations, with an 8° angle and containing 9000
particles.

For the simulations in the next chapter (4.5. Experimental blast furnace model) a 3D
geometry was modelled to investigate the deadman size and shape. Due to the time
required to run this simulation only one 3D simulation was done. The results are presented
here as a comparison, even though the geometries are slightly different. The 3D model
contains 105000 spherical particles with a diameter of 0.04 m. If the 0.022 m diameter of
the other simulations would have been used the number of particles becomes too large to
solve within a reasonable time. In all the models the particles are removed at a pre‐defined
rate from a volume representing the raceway. For the 3D case this is a ring‐shaped volume,
not separate raceway zones, to allow accurate comparison with the other models

Table 4.3. Geometry comparison

(a) (b) (c) (d) (e)
Half Double Full Pie‐slice Full 3D
diameter Thickness diameter
Thickness (m) 0.04 0.075 0.04 ‐ ‐
Pie angle ‐ ‐ ‐ 8° ‐
No. of Particles 9 000 18 000 18 000 9 000 105 000

Boundary conditions and simulation procedure

New software developments allowed the use of periodic boundary conditions for the front
and back planes. Using these conditions particles moving through the back wall will re‐
enter the simulation at the same position at the front wall. This creates much more realistic
boundary conditions at the front and back walls compared to the previous cases. The
volumetric particle removal has a pre‐defined rate for each case ensuring the same draw‐
down rate. The cases are compared after 345 seconds of real time and each case has the
same amount of layers removed. Particle properties are the same as have been used in
Case II and are presented in Table 4.3.

Results and discussion

The earlier simulations of the pie‐slice model have shown the influence of the presence of
front and back planes in the simulation. Using a slot model allows us to use periodic
boundary conditions but is less representative of the actual 3D blast furnace shape. When
compared to a full diameter slot model, simulating only half of the diameter is a convenient
way of halving the amount of particles in the simulations, assuming that the results can be
mirrored along the centre line. This is shown in Figure 4.7(a) and (b) where the left wall, the

88
centre line, is a frictionless wall. To investigate the influence of the thickness of the slot
model, it was doubled from 0.04 m to 0.075 m, doubling the amount of particles from 9000
to 18000. This has little influence on the Deadman size as shown in Figure 4.7(a) and (b),
and on the solid flow lines (not shown here). However, a thicker slot is always preferred to
lower the wall effects on the particle flow. A full diameter slot model geometry with 0.04 m
thickness again doubles the amount of particles to 18000. The result in Figure 4.7(c) shows
that restriction of flow in horizontal direction by the frictionless bounding wall at the
furnace centreline does significantly influence the flow. There is a clearly distinguishable
deadman present in Figure 4.7(a) and (b), which has mostly disappeared in Figure 4.7(c).
The 3D case in Figure 4.7(e) confirms the absence of a deadman when modelling cold solid
flow. The pie‐slice geometry in Figure 4.7(d) with realistic front and back plane conditions,
as described for Case II, also shows good agreement to the results of Figure 4.7(c) and (e).
The geometry, however, has several problems; it has no periodic boundary conditions and
the sharp edge along the centre line means that the area remains empty or is removed
from the simulation. The former would be highly unrealistic if gas flow is included in the
simulation. Therefore, the full diameter slot geometry is considered most appropriate for
the simulations, confirmed by the full 3‐D model results.

89
 (a) (b) (d)

(c) (e)
Figure 4.7. Geometry comparison, coloured by residence time:
(a) half diameter, (b) double thickness, (c) full diameter, (d) pie‐slice, (e) full 3D

90
4.3.4. Experimental blast furnace model
In the previous subsection it has been determined what geometry approximation gives the
optimal balance between realistic results and computation speed within the constraints
given by the software. To increase the realism of the simulation the geometry for the
present case is based on an experimental blast furnace. The smaller size of the furnace
allows us to simulate the actual process rather than downscaling parameters as done when
simulating the industrial blast furnace. It also opens up the opportunity to validate the
results with data from the experimental furnace when the model has been developed
further.

In the previous simulations gas flow has not been included due to the instability of the
DEM‐CFD coupling module. All attempts to use the coupling module failed due to program
crashes. A new coupling module became available which was found to be much more
stable and gas flow was included in this case. The case focuses on the influence of gas flow
on the particle flow: a simulation without gas flow is compared to two cases with increasing
gas velocity. In the simulations both spherical pellets and non‐spherical coke are present.
Coke is removed from the raceways but pellets shrink and are removed in a pre‐defined
cohesive zone. Both the non‐spherical coke and the pellet shrinking and removal in the
cohesive zone increase the realism of the simulations.

Geometry and particles

The experimental blast furnace, EBF, on which this case is based is considerably smaller
than an industrial blast furnace, allowing full scale simulations of the process with
possibility for validation. It is based on the pellet producer LKAB’s experimental blast
furnace at Swerea MEFOS, shown in Figure 4.8, and is used to test blast furnace pellets and
operational conditions [20‐22]. The EBF has a hearth diameter of 1.2 m, a working volume
3
of approximately 8.2 m and produces approximately 35 tonnes of hot metal per day. On
average 2 campaigns of 50 days are done every year. Most research work has been done on
product development of pellets, and thus focussed on the mineralogy and chemistry of the
burden.

91
 Figure 4.8. LKAB Experimental blast furnace [23]

For our model 4 m of the total 7 m of the EBF is simulated: from the middle of the stack, 1.6
m below burden level, to the top of the hearth. This case is a slot model with a geometry
thickness of 0,06 m and periodic boundary conditions on the front and back walls. It uses
the full diameter geometry which the previous chapter has shown to be most realistic. Any
particles protruding or leaving the geometry through the back wall will enter through the
front wall. Particles in this geometry are removed from two pre‐defined volumetric zones:
coke from the raceway and pellets from the cohesive zone. Both the raceway and cohesive
zone geometries are based on the actual geometry of the experimental furnace [20, 22,
24]. After most of the campaigns of the EBF it has been quenched with nitrogen and
excavated. From these results the location and shape of the cohesive zone is known for
different operational conditions. Figure 4.9 shows one example of the cohesive zone
location and shape. In our model the cohesive zone is similar to those found in the EBF; the
top is 3.4 m under the burden level, the height is 0.95 m and the thickness is 0.3 m. The
position and shape of the cohesive zone is strongly dependent on the burden material
composition, e.g. melting point and reducibility. These are not considered in our
simulations. Ideally the parameters from the EBF data should match those from our
simulations.

92
 Figure 4.9. EBF Cohesive zone location [24]

Figure 4.10(a) shows the shape of the furnace and the cohesive zone and raceways. In the
blast furnace the cohesive zone is defined as the zone starting where the pellets begin to
soften and end where all pellets have molten. Melting in the cohesive zone is simulated by
a stepwise decrease of the particle diameter. Each particle passing one of the two first lines
defined in the cohesive zone will be replaced by a smaller particle. When passing the last
line the particle will be removed from the simulation. There is no actual melt formation
taking place in the simulation. Softening which causes shape change of the particles is not
taken into account here.

93
 4m Cohesive zone

Raceway

(a) (b)
Figure 4.10. Blast furnace geometry (a) and CFD mesh (b)

For the CFD simulation the geometry is divided into a grid of 67×22 cells and one cell deep
shown in Figure 4.10(b). Gas is injected into the raceways at 20 and 40 m/s through inlets in
the walls, indicated in Figure 4.10(b) by the thick black lines. In a real blast furnace the
velocity will be considerably higher, mainly due to the high temperatures at the raceway
greatly increasing the volume of gas. The gas velocity at the top of the furnace is, however,
close to realistic values.

As in the previous cases a combination of spherical and non‐spherical particles shown in

Figure 4.11 is used to represent pellets and coke. The non‐spherical coke consists of 6
particles with a 0.016 m diameter in a bi‐pyramidal shape; the widest diameter is 0.036 m.
To simulate the ore pellets melting during the simulation the spherical particles shrink in
two steps in the pre‐defined cohesive zone; from 0.02 m to 0.015 m and then to 0.01 m,
they then disappear. At the start of a simulation the packed bed is filled only with coke on
which the first layer of pellets is charged. The coke particles are generated across a large
part of the width of the furnace, the pellets in a small area close to the side wall. This
creates pellet layers which are becoming thicker close to the wall and correspondingly
thickening coke layers in the centre. The created layer structure is close to the structure
present in the actual blast furnace where in the centre‐line of the furnace only coke is
charged to ensure good gas flow in the furnace.

94
 Pellet, 0.02 m Coke, 0,036 m
Figure 4.11. Simulation particles

Boundary conditions and simulation procedure

The particle properties are shown in Table 4.5 and are based on generic values for coal from
the database of the EDEM software package used for these simulations. Coupling DEM and
CFD is done with the EDEM Coupling Module using the Eulerian coupling method with the
Di Felice drag mode [8]; no models for lift and heat transfer were used. After 50 DEM time‐
steps a Fluent CFD time‐step is calculated. For the CFD simulation the standard k‐ε model
for turbulence is used. As an approximation, gas properties for carbon monoxide are used.

Three simulations have been done to investigate the influence of the gas flow on the
particle flow. The first simulation does not contain any gas flow, in the second the gas
velocity is 20 m/s at the inlets which is increased for the third simulation to 40 m/s. The first
is a DEM model and the latter two are DEM‐CFD coupled models. Simulations are started
by filling the geometry with coke; when filled, particle removal is started from the raceways
and the first layer of pellets is charged.

Table 4.5. DEM Parameters

Coke Pellet
Poisson ratio 0.25 0.25
Shear modulus, Pa 1×107 1×107
Density, kg/m3 1000 4000
Coefficient of restitution 0.2 0.2
Coefficient static friction 0.5 0.2
Coefficient rolling friction 0.05 0.02
DEM Time step, s 5×10‐5

95
Results and discussion
The DEM particle simulation results are shown in Figures 4.12 and 4.14. Results from the
CFD simulations in Figure 4.13 naturally only include the gas flow cases.

Layer structure
From the top row of cross sections shown in Figure 4.12(a) the resulting layer structure can
be seen. The layers become more upward curving at the centre line of the furnace with
increasing gas velocity. From Figure 4.13(a) it can be seen that the gas velocity through the
centre is considerably higher than it is closer to the side‐walls. Figure 4.13(b), which shows
the solid fraction in the CFD cells, shows that the porosity of the coke layers above the
cohesive zone and the coke bed below it, is higher than that of the pellet layers. This is
caused by the non‐spherical particle shape of the coke particles. The higher gas velocity
increases the resistance on the solid flow and slows it down causing the curved layers. The
largest influence of the gas flow on the solid layers is the splitting up of the pellet layers
resulting in a coke filled centre. This is caused by both the increase in gas flow at the
tuyeres and its increase by the decrease in porosity of the pellet layers. The layer structure
of the 20 m/s gas velocity case is skewed downwards in the left part of the furnace
indicating a higher solid flow velocity. The addition of gas flow to the simulation decreases
the stability of the particle flow by loosening the structure. It should be remarked here that
in the real blast furnace the expansion and contraction of the gas due to the high
temperature difference will have a significant influence on the materials flow.
Comparison of Figure 4.12(a) to the data from the LKAB EBF in Figure 4.9 shows
a similar layer pattern above the cohesive zone, with slightly dipping layers towards the
centre. On the level of the cohesive zone the layering is different; in the EBF the layers start
to dip stronger, in our model they are curved. This is caused by the way the cohesive zone is
defined in our model. Its location is pre‐defined and particles are removed based on
location, not on local process conditions. The softening and melting rate differs in the
cohesive zone; this model does not take this into account as it requires thermodynamic and
kinetic models. With which it could also determine the location of the cohesive zone. The
parallel work which has been done in the total project has successfully managed to couple
the gas and solid flow model with a softening and melting, thermodynamic equilibrium and
reaction kinetic models. This full coupled model has the ability to simulate the cohesive
parameters. More details on this model are given in Appendix A.

Figures 4.14(b) and (c) shows the residence times of the particles, as shown in previous
chapters there is no deadman present in these conditions. Simulations including gas flow
require a longer simulation time with increasing gas flow due to the resistance which the
descending particle flow encounters from the ascending gas flow. As mentioned in the
previous subsection 4.4.3, a 3D simulation of this case was done and shown in Figure
4.14(a). It models half of the cylindrical blast furnace. The 3D case uses spherical particles
with a diameter of 0.04 m, similar in size to the non‐spherical particles in the slot model
simulation in Figures 4.14(b) and (c).

96
The use of spherical particles was necessary in the 3D simulation because the non‐spherical
particles consist of 6 particles, which would bring the number of particles up from 105 000
to 630 000. The 3D case is uncoupled with CFD and does not include any gas flow. In the
cohesive zone in Figures 4.14(b) and (c) the particle removal can be seen by the newly
created particles which have a much smaller residence time than the surrounding particles.

There is significant amount of hold‐up along the diagonal side‐walls under the cohesive
zone. This effect can also be seen in the 3D simulation, but in a much smaller degree. The
effect is caused by the diagonal side‐walls restricting flow and by the non‐spherical particle
shape which restricts particle rotation. The geometry thickness increases these effects; it is
only 0.06 m, which is nearly four times the smallest particle, but less than two times the
diameter of the largest combined particles. Even though the front and back have periodic
boundary conditions, the small thickness negatively influences the particle flow in this
region.
Besides this hold‐up along the side‐walls, the results are very similar to the 3D
simulation. The lines visible in the residence time distribution show the influence of the
burden layer structure caused by the charging of pellets closer to the wall and coke in the
centre. It is also influenced by the higher gas flow through the centre which can be seen in
Figure 4.13(a).

97
No gas flow Low gas flow (20m/s) High gas flow (40m/s)

(a) DEM Particles: Red ‐ pellets, Grey ‐ coke

350N

280N

210N

140N

70N

0N

(b) Compressive force

Figure 4.12. DEM Simulation results of the EBF slot model;
(a) layer structure, (b) particle compressive forces

98
(a) High gas flow (40 m/s) (b) Solid burden fraction

Figure 4.13. Gas velocity and porosity results of the EBF slot model

99
 50 s

40 s

30 s

20 s

10 s

0s

(b) Low gas flow

(a) 3D Model
(20 m/s)

75 s

60 s

45 s

30 s

15 s

0s

(c) High gas flow (40 m/s)

Figure 4.14. Residence time

100
Compressive forces
The second row of cross sections in Figure 4.12(b) shows the compressive forces on the
particles creating a dendritic network structure on which the material is resting. In all cases
the force network bridges the width of the furnace and is supported by the inward sloping
side walls. The network of high compressive forces in the simulation without gas flow is
considerably larger compared to those in both cases with the gas flow. The material in the
furnace is supported by a larger area of the inwards sloping side wall. Injection of gas in the
raceway causes the expansion of an area where the particle packing is loose. From Figure
4.12(b) it can be seen that this area increases upward when gas flow is added and then
increased.

Results from literature show a very similar network structure, as demonstrated in Figure
4.15(a), (b) and (c) [25‐27] . These models show relatively wide geometries compared to the
EBF presented in Figure 4.12(b). The build‐up of the force network from the side‐walls is
considerably less. Figure 4.15(a) shows the highest stresses in the centre of the hearth,
following the profile of the cohesive zone. In Figure 4.15(b) there some network arms
reaching to the wall, but the highest forces are also in the centre. The same is seen in the
results from Figure 4.15(c). Because the width of the EBF in Figure 4.12(b) is relatively
small, the raceway is relatively large and there is less room for a packed bed to build‐up in
the centre on which the burden rests. In the work of Fan et al. [25] it is also concluded that
the stress distribution is dependent on the inner volume of the furnace. In large blast
furnaces normal stress is reduced near the walls due to the higher particle velocity there.
Figure 4.15(c) shows the change in normal stress at different gas velocities, the
results seem to agree well with those in Figure 4.12(b). In both cases the increase in gas
flow causes a decrease in the normal forces or compressive forces on the particles.

CFD Results
Figure 4.13(a) shows the velocity profile of the gas flow from the CFD simulation. The
porosity in the pellet layers is lower compared to that of the coke layers due to the non‐
spherical shape of the coke particles which creates larger pores between the particles, as
shown in Figure 4.13(b). The higher porosity is also present along the centre line of the
furnace due to the larger amount of coke present there as a result from the particle
charging; and below the cohesive zone where all the pellets are removed. The higher
porosity of the coke bed in the hearth and along the centre line of the furnace combined
with the inward velocity of the gas result in a higher gas flow through the centre of the
furnace. In the part of the burden where pellet and coke layers alternate the gas velocity is
higher in the pellet layers due to the lower porosity. The average gas velocity in the packed
bed is between 2‐8 m/s; at the top this is close to reality but below the cohesive zone it is
far too low due to the isothermal conditions.

101
 (a) (b)

1.21×10‐2 kg/s 3.62×10‐2 kg/s 6.05×10‐2 kg/s

(c)
Figure 4.15. Normal stress distribution: (a) Fan et al. [25], (b) Nouchi et al. [26],
(c) Zhou et al., at different gas flow rates [27]

102
Discussion
This case investigates the interaction between the gas flow and the solid flow as a packed
bed. Pellet layers have a lower permeability than the coke layers due to the non‐spherical
coke shape. The method of charging pellets to the wall and coke across the width increases
the coke‐pellet ratio towards the centre‐line. Combined with the pre‐defined shape of the
cohesive zone, gas will flow preferentially through the centre of the furnace. At increasing
gas velocities this causes the pellet layers to curve upward at the centre‐line and eventually
the pellet layers separate in the centre. The created coke filled heart will cause even more
gas to flow through the centre.

103
4.3.5. Charging and burden descent models
All the previous cases were aimed at modelling the complete blast furnace; the DEM
technique is also very suitable for simulation of smaller problems. Two such problems are
related to the charging of the blast furnace and burden descent. Because the blast furnace
is completely filled with a solid burden the structure of this burden is very important to the
process inside the furnace. Material is charged into the furnace using a rotating chute with
an adaptable angle to allow control over the location of charging. This is used to build up a
burden with the required composition and structure. In general this is a layering of coke and
ore layers, with a decreasing ore/coke ratio to the centre. The behaviour of the charged
material as it is deposited on the burden is strongly dependent on the material properties,
burden structure and the chute configuration. Process irregularities can cause accretions on
or wear of the furnace lining. This influences the solid flow and layer structure potentially
resulting in process disturbances. In parallel with the other work, simulations have been
done to investigate the use of DEM for such problems.

Geometry and particles

Burden charging model

The geometry of the charging simulation is a slot model. It has a width of 4 m, a realistic
size to simulate the radius of an industrial blast furnace, depth of the slot is 0.120 m and the
height is 4 m. In the slot geometry there are four zones where the particles are created to
simulate the rotating chute moving inward. They are 0.2 m wide and are placed 2.5, 2.9, 3.3
and 3.7 m from the centre, as illustrated in Figure 4.16. If the rotating chute moves inward
the circle it describes becomes smaller, if the charging rate is constant more material will
be charged per square meter. Therefore, the number of particles charged from a zone per
time, increases from the wall to the centre; 1044, 1170, 1332 and 1545 for pellets and 351,
393, 448 and 519 for coke. Non‐spherical coke has the same 6 particle bi‐pyramidal shape
as used in the previous simulation approximately 0.036 m in diameter, and pellets are
0.02 m spheres.

Burden descent model

For the burden descent simulation the full experimental blast furnace was used as in the
previous chapter. For this case the total experimental blast furnace height of 6.8 m is
modelled, the thickness of the slot model is 0.12 m. The four geometries with four wall
obstructions and wear are shown in Figure 4.17; (b) and (c) show accretion in the upper
shaft and (d) in the lower, and (e) shows wear of the lining in the upper shaft. Particles are
spherical pellets and non‐spherical coke, the latter has a four particle pyramidal shape with
a widest diameter of 0.036 m. This is more than three times the particle diameter, and
twice that of the EBF simulations in the previous chapter. Pellets (0.02 m diameter) shrink
in four steps (0.0159 m, 0.0126 m, 0.01 m and 0.0079 m) each step halving its volume,
before disappearing in the cohesive zone as shown in Figure 4.17(a). Coke is removed from
the raceway as was done in the previous cases. Pellets are charged closer to the wall and

104
coke more across the furnace width, each pellet layer contain 1820 particles and each coke
layer 490 particles. The resulting layers are more or less of equal thickness.

Figure 4.16. Burden charging geometry and initial burden

(a) (b) (c) (d) (e)

Figure 4.17. Geometry, wall accretions and wear

105
Boundary conditions and simulation procedure

Burden charging model

The main difference in the boundary conditions between the two cases shown in Table 4.6
is the shear modulus. For the burden descent case the same shear modulus as the other
cases is used. This is an unrealistically low value which does not have significant influence
on results when simulating a packed bed, but greatly speeds up simulations. For the
charging simulation where the collision of particles is very important, a shear modulus of
7
1×10 is too low and gives unrealistic results. Due to the impact of the heavy pellets, the bed
deforms and then springs back in form. Therefore, the shear modulus for the charging
8
simulation is set to 1×10 . The charging models were simulated with and without gas flow,
burden descent is modelled without.

Table 4.6. DEM Parameters

Charging Burden descent
Coke Pellet Coke Pellet
Poisson ratio 0.5 0.5 0.35 0.35
Shear modulus, Pa 1×108 1×108 1×107 1×107
Density, kg/m3 1400 4000 1000 4000
Coefficient of restitution 0.5 0.5 0.5 0.5
Coefficient static friction 0.6 0.5 0.76 1
Coefficient rolling friction 0.05 0.01 0.05 0.05
DEM Time step, s 2.45×10‐5 5×10‐5

The influence of gas flow through the bed on charging was investigated first. Simulation
starts with the pre‐defined coke bed shown in Figure 4.18(a), the furnace wall is at the right
and the centre at the left. The slope of the coke bed is a rather steep and has a well‐defined
angle, a situation not very realistic but present here to investigate the principle. On the bed
16 charges of pellets are dropped at the four locations described before, starting at the wall
and moving inwards. Gas flows through the bed at 5 m/s from the bottom. On the model
including gas flow a coke layer is charged followed by another pellet layer.

Burden descent model

In the burden descent simulations the full furnace is filled with coke after which particle
removal from the raceways is started. At the start of simulation the first layer of pellets is
charged. As the top level drops below the original level alternating layers of coke and
pellets are charged. The simulation is considered to be finished after the volume above the
cohesive zone has been renewed twice.

106
Results and discussion

Burden charging model

The influence of gas flow on the charging simulation can be seen in Figure 4.18(b) and (c).
Material from the coke bed is being carried with the pellets to the left; this causes the
distinct angle present at the start to smoothen out. In the case with gas flow in Figure
4.18(c) much more of the angle is removed; besides the particles being dragged with the
flow of pellets, there is also a layer of material being pushed down. This is caused by a
loosening of the structure of the bed by the upward forces the gas flow exerts on the
particles. The round particle shape of the pellets translates into a low angle of repose which
can be seen in Figure 4.18(b) and (c), the angle is much lower than that of the non‐spherical
coke underneath.

On the previously charged layer of pellets another coke layer is charged shown in Figure
4.18(d). The case containing gas flow, Figures 4.18(c), (d), and (e), is the more realistic as it
is virtually always present in industry and therefore the simulation is continued with that
case. Figure 4.18(d) shows the layer has a very similar shape as the starting layer. The slope
of the layer is again steeper than the underlying layer of pellets due to its non‐spherical
shape causing a higher angle of repose. The slope levels off at the point where the last
amount of coke was charged. In Figure 4.18(e) another layer of pellets has been charged, as
in Figure 4.18(c) a layer of coke is being pushed down. The pellet layer on which the coke
layer is now laying has a much smaller angle of repose than the angle of internal friction or
repose of the coke. This causes the coke layer to slide down the pellet layer rather than
movement of particles in the coke layer itself. A four times higher density of the pellets
amplifies the effect of the above. The work of Ho et al. [18] shown in Figure 4.19
demonstrates the same effect, the impingement of dense particles is much greater than
less dense particles.

From the simulation results in the previous subsection we have seen that the layer porosity
in combination with the layer structure has a large influence on the gas flow in the blast
furnace, shown in Figures 4.12 and 4.13. The simulations here show that this layer structure
can be strongly influenced by the charging method. The burden level after the coke charge
in Figure 4.18(d) is very similar to the initial burden in Figure 4.18(a). The next pellet charge
in Figure 4.18(b) changes the underlying coke layer completely. In order to simulate a full
scale industrial blast furnace process, as is the ultimate aim of the project, the burden
charging should be taken into account in the model. For the EBF simulations the ratio
between the layer or particle size, and blast furnace throat diameter for these effects to
have a significant influence.

107
 (a)

(b) (c)

(d) (e)
Figure 4.18. Pellet charging;(a) initial conditions (b) no gas,
(c,d,e) with gas flow

Figure 4.19. Burden level after charging wood, glass and steel particles
(left to right)[18]

108
It is clear that the pellet properties such as shape and friction coefficients play an important
role in these simulations. When these are experimentally determined for the particles in the
simulation, the burden charging simulations can give accurate details on how the material
is distributed in the furnace. If this is combined with simulations of the stockhouse, skip,
top bunkers and charging chute for effects such as size segregation an accurate description
of exactly what ends up where in the blast furnace can be obtained. Using this model an
optimal burden distribution and layer structure can be achieved.

Burden descent model

Figure 4.20 presents the results of the burden descent simulations. The first accretion in
Figure 4.20(a) shows only a small influence on the layer structure. Below the accretion the
layers mix and dip towards the wall. Mixing the large coke particles with smaller pellets is
likely to decrease the porosity. This is, however, highly dependent on the size ratio
between the particles and the mixture ratio. A much larger effect can be seen with a larger
accretion as seen in Figure 4.20(b). The most visible effect is the presence of a large void
immediately under the accretion. Mixing is present as in (a) but is stronger and the mixed
layer extends further into the furnace. Next to the mixed layer there is a vertical layer with a
higher pellet/coke ratio. Because the pellet layers have a lower porosity this creates a
higher resistance to gas flow in this area. The formation of this layer starts above the
accretion where descent is blocked and material pushed inwards. Because the spherical
pellets have a low rolling friction and the non‐spherical coke a high one, the pellets are
pushed inward from between the coke layers. When the material passes the accretion coke
flows back to the wall mixed with pellets and the pellets which were pushed inwards remain
on their location. In Figure 4.20(a) the layer structure remains in a large part unaffected, in
(b) the structure in the whole furnace is influenced by the accretion.

If the accretion is placed lower in the furnace as in Figure 4.20(c) the structure in the upper
part remains unaffected. When the layers approach the accretion the material is held up
and the layers start curving. The accretion in (c) is lower and longer than that of (b), this
gives a gentler transition from the wall causing the layers to curve around the accretion
rather than breaking them up as in the previous simulation. Below the accretion there is
strong mixing of the materials. Because the area of this effect is quite large and it coincides
with the location of the softening and melting it will have significant effects on the cohesive
zone behaviour. This is highly important in an industrial blast furnace as it could block gas
flow and disrupt the process.
In Figure 4.20(d) the effect of wear of the furnace lining is shown. Widening the
furnace gives the particles additional room to flow into, which causes the layers to slope
downward to the left wall. At the end of the worn section the material is held up, slopes
steeply upward and gets mixed. This mixed layer next to the wall continues down to the
cohesive zone. The influence of the mixing is quite small in this case, but the leftwards
sloping layer structure is present everywhere in the furnace.

109
These simulations demonstrate the possibilities DEM gives to investigate solid flow in the
blast furnace. The full height experimental blast furnace is modelled here, though it is still a
slot model and particles are larger than in reality. Because gas flow was not included in
these simulations the calculation is faster and the geometry could be larger. However, each
simulation still cost a week to solve.

The formation of wall accretions is a possible process risk in the blast furnace which can
cause irregular burden descend and an inefficient process. They can be formed due to
build‐up of materials in the furnace which can form these accretions, such as alkalis or zinc.
As is also shown in the results in Figure 4.20, they can have a large effect on the layer
structure. These simulations can supply knowledge and an understanding on the effect of
wall accretions or wear, at least for the part above the cohesive zone. As mentioned before,
the results are hard to validate with real data. The open question is what the influence of
the cohesive zone on solid flow above it is, and how far it extends into the upper furnace.
This determines how far the burden can be simulated without taking into account details
on the cohesive zone, and thus the validity of the results of this model. This again shows
the need for the comprehensive blast furnace model which this project has been
developing.

110
(a) (b) (c) (d)
Figure 4.20. Burden descent simulations
with wall accretions and wear

111
4.4. Concluding remarks
This chapter has described the development of the coupled CFD‐DEM model, part of the
project to develop a comprehensive model which can simulate the blast furnace cohesive
zone. Both the DEM method and the coupling with CFD show promising results and
demonstrate the value of DEM simulations for the blast furnace.
The main issue is the simplification of the geometry required to do the
simulations in a reasonable time. However, improvements in the EDEM software have seen
the simulations increase in size, as can be seen in the simulations in this chapter. Continued
development of the software will make this less of an issue.
For the cases with gas flow, a thermal profile of the furnace is missing. Due to the
very high temperature drop over the furnace the gas density will drop significantly, greatly
influencing the gas flow through the furnace, and thus the solid flow. Simulations have
been done in which there was a profile predefining the energy loss over the furnace. A large
part of the heat is removed in the cohesive zone due to melting of the burden materials. An
example can be seen in Figure 4.21., where a 1200 K blast enters the furnace and an energy
sink is defined in the cohesive zone to simulate melting of the material. Adding
temperature to the simulations would greatly improve the realism.

Figure 4.21. DEM‐CFD Model with hot blast

and pre‐defined cohesive zone energy sink

112
References
1. EDEM 2.3 User Guide. Revision 2C ed. (2010), DEM Solutions.
2. FLUENT 6.3 User's Guide. (2006), Fluent Inc.
3. EDEM‐CFD Coupling for FLUENT User Guide. Revision 4B ed. (2010), DEM
Solutions.
4. Mindlin, R.D. and Deresiewicz, H., Elastic spheres in contact under varying oblique
forces. Journal of Applied Mechanics 20 (1953), p.
5. Zhu, H.P., Zhou, Z.Y., Yang, R.Y., and Yu, A.B., Discrete particle simulations of
particulate systems: Theoretical developments. Chemical Engineering Science 62
(2007), p. 3378‐3396.
6. Gidaspow, D., Multiphase flow and fluidization. (1994), Academic press.
7. Feng, Y.Q. and Yu, A.B., Assessment of model formulations in the discrete particle
simulation of gas‐solid flow. Industrial & Engineering Chemistry Research 43
(2004), p. 8378‐8390.
8. Di Felice, R., The voidage function for fluid‐particle interaction systems.
International Journal of Multiphase Flow 20 (1994), p. 153‐159.
9. Zhou, Z., Zhu, H., Yu, A., Wright, B., Pinson, D., and Zulli, P., Discrete particle
simulation of solid flow in a model blast furnace ISIJ International 45 (2005), p.
1828‐1837.
10. Ueda, S., Nogami, H., Kano, J., and Ariyama, T. Development of blast furnace
operation simulator including discrete element concept for solid motion. in Recent
progress on mathematical modeling in ironmaking. (2008), Tokyo, Japan.
11. Zhou, Z.Y., Pinson, D., Zhu, H.P., Yu, A.B., Wright, B., and Zulli, P. Discrete
particle simulation of gas‐solid flow in a blast furnace. in Third International
Conference on CFD in the Minerals and Process Industries. (2003), Melbourne,
Australia.
12. Nouchi, T., Takeda, K., and Yu, A.B., Solid flow caused by bouyancy force of heavy
liquid. ISIJ International 43 (2003), p. 187‐191.
13. Kawai, H. and Takahashi, H., Solid behavior in shaft and deadman in a cold model
of blast furnace with floating‐sinking motion of hearth packed bed studied by
experimental and numerical DEM analyses. ISIJ International 44 (2004), p. 1140‐
1149.
14. Mio, H., Yamamoto, K., Shimosaka, A., Shirakawa, Y., and Hidaka, J., Modeling of
solid particle flow in blast furnace considering actual operation by large‐scale
discrete element method. ISIJ International 47 (2007), p. 1745 ‐ 1752.
15. Mio, H., Komatsuki, S., Akashi, M., Shimosaka, A., Shirakawa, Y., Hidaka, J.,
Kadowaki, M., and Matsuzaki, S., Validation of particle size segregation of sintered
ore during flowing through laboratory‐scale chute by discrete element method. ISIJ
International 48 (2008), p. 1696‐1703.
16. Nogami, H., Toda, K., Pintowantoro, S., and Yagi, J.‐i., Cold‐model experiments on
deadman renewal rate due to sink‐float motion of hearth coke bed. ISIJ
International 44 (2004), p. 2127‐2133.
17. Yuu, S., Umekage, T., and Kadowaki, M., Numerical simulation of particle and air
velocity fields in raceway in model blast furnace and comparison with experimental
data (cold model). ISIJ International 50 (2010), p. 1107‐1116.
18. Ho, C.K., Wu, S.M., Zhu, H.P., Yu, A.B., and Tsai, S.T., Experimental and numerical
investigations of gouge formation related to blast furnace burden distribution.
Minerals Engineering 22 (2009), p. 986‐994.

113
19. Yu, Y. and Saxén, H., Experimental and DEM study of segregation of ternary size
particles in a blast furnace top bunker model. Chemical Engineering Science 65
(2010), p. 5237‐5250.
20. Hooey, L., Sköld, B.‐E., Sundqvist Ökvist, L., Seppänen, M., and Zuo, G. LKAB's
experimental blast furnace ‐ the learning curve. in The 5th European Coke and
Ironmaking Congress. (2005), Stockholm, Sweden.
21. Hallin, M., Hooey, L., Sterneland, J., and Thulin, D., LKAB's experimental blast
furnace and pellet development. La Revue de Metallurgie (2002), p. 311‐316.
22. Hallin, M. BF Inside story through quenching. in John Floyd International
Symposium on Sustainable Developments in Metals Processing. (2005),
Melbourne.
23. Zou, G. and Hirsch, A., The Trial of the Top Gas Recycling Blast Furnace at LKAB's
EBF and Scale‐up. La Revue de Metallurgie (2009), p. 387‐392.
24. Hallin, M. Evaluation of ferrous burden properties in an experimental blast furnace
after quenching and dissection. in International BF Lower Zone Symposium
(2002), Wollongong, Australia.
25. Fan, Z., Natsui, S., Ueda, S., Yang, T., Kano, J., Inoue, R., and Ariyama, T.,
Transient behavior of burden descending and influence cohesive zone shape on solid
flow and stress distribution in blast furnace by discrete element method. ISIJ
International 50 (2010), p. 946‐953.
26. Nouchi, T., Sato, T., Sato, M., Takeda, K., and Ariyama, T., Stress field and solid
flow analysis of coke packed bed in blast furnace based on DEM. ISIJ International
45 (2005), p. 1426‐1431.
27. Zhou, Z.Y., Zhu, H.P., Yu, A.B., Wright, B., and Zulli, P., Discrete particle
simulation of gas‐solid flow in a blast furnace. Computers and Chemical
Engineering 32 (2007), p. 1760‐1772.

114
Chapter 5
Conclusions and recommendations

5.1. Conclusions
The work done in this thesis has developed a model based on the Discrete Element Method
(DEM) in combination with Computational Fluid Dynamics (CFD), which can be used to give
insight into the process conditions inside the blast furnace and more specifically the
cohesive zone, having interactions with ascending gas flow simulated with CFD.

The use of DEM modelling has increased greatly in recent years due to the increase in
computing power. For the application of DEM on blast furnace ironmaking three types of
simulations can be found: the burden flow, the raceway and burden charging. At the
current level of development of DEM models it can be most readily applied to the last type,
the burden charging. For both raceway and burden flow the complexity is much greater
and requires additional sub‐models.
From literature a number of models can be found which have been used to
simulate the burden flow in the blast furnace. Some include gas flow and some also include
a liquid hearth level. The work from literature focusses on solid flow patterns and the
resulting layer distribution of the burden. Simulations show the development of a stagnant
deadman and a dendritic force network, highest in the deadman and hearth. All the models
require a number of simplifications to run; the full scale of an industrial blast furnace is too
large to simulate, if real physical size of burden materials (pellets, sinter, and coke) are
used.
Due to its complexity the cohesive zone, when studied by modelling in the
literature, is pre‐defined and simplified. This is where the present study tries to improve on
the current state of modelling. The end goal is to simulate the behaviour of the cohesive
zone based on the temperature profile in the furnace, the chemical composition and
reduction reactions of the burden, and the softening and melting behaviour of the ferrous
materials.

Because of the large computational effort required, the geometry needs to be simplified to
reach a workable number of particles. The first models have made it very clear that the
results are very strongly influenced by the geometry used. A full scale cylindrical blast
furnace is too large to be modelled, but the choices made to simplify it can have a large
effect on the observed solid flow. The DEM modelling technique is not yet efficient enough
or the computational power large enough to be able to run larger scale models which make
these size and shape simplifications unnecessary.

The DEM method uses perfectly spherical particles in the simulations, which in reality is
never the case. In the present study non‐spherical particles are created by combining
spherical ones and using these to approximate the blocky coke particles. This shows a large

115
influence on the particle flow and resulting layer structure. It is shown that the non‐
spherical particle layers have a larger resistance to deformation than the spherical particle
layers. The shape of a non‐spherical particle restricts the rotation of the surrounding non‐
spherical particles. Clusters of non‐spherical particles are also more resistant to shear due
to the particles restricting each other’s movement by interlocking. This results in the initial
layer structure as charged in the top to remain more constant while descending through
the furnace.

An investigation of the optimal simplified shape of the geometry concludes that the full
diameter slot model is the best trade‐off between the amount of particles and realism. The
pie‐slice model has the major disadvantage of becoming an edge at the centre line of the
furnace when using a small angle for the geometry. This edge and the decreasing thickness
of the geometry towards it create problems for the model: the walls have a very large
influence on the flow, especially when the thickness will become one particle diameter and
less. The pie‐slice model would be much more useful for a large angle geometry, such as a
90 slice. The second geometry shape which was investigated was related to the diameter
of the slot model. The amount of particles can be halved by simulating only half of the
furnace diameter. Results show that the centre wall along which the geometry is halved has
a large influence on the flow. Its presence restricts the flow and causes the formation of a
small deadman not being reflected in the 3D simulation or full width simulations which are
thought to be closer to the reality.

The present study indicates that the ascending gas flow through the furnace is strongly
influenced by the layer structure. By charging more coke into the centre of the furnace the
porosity is higher, causing more gas to flow through the centre line of the furnace and at
higher velocities. This breaks up the pellet layers and curves the layers upwards. Gas flow
causes upward forces on the particles and loosens the structure of the bed; the
compressive force network in the hearth is much smaller than in the area around the
cohesive zone. These effects are highly dependent on furnace or model scale, and gas
velocity. A full scale industrial blast furnace has a much larger and thus heavier burden
pushing down. The gas flowing through an industrial furnace has a much higher
temperature and pressure. If insight is to be gained from these models, these factors have
to be taken into the model.

During the development of the comprehensive DEM‐CFD coupled model, it has been
continuously applied to simulate solid flow and coupled solid and gas flow cases. The result
of each of these cases has been the basis for the following one; thus the simulations have
steadily grown in complexity. The conclusions of the cases with increasing complexity can
be listed as follows:

116
For the lab scale slot models:
 Three identified solid flow regions: plug flow in the shaft, funnel flow between
the deadman and the wall, and a low velocity deadman
 Good agreement for the general flow patterns with experimental results from
literature, but the simulated deadman size is comparatively smaller than I
experiments and there is less particle mixing
 Flow pattern strongly influenced by the boundary conditions and model
parameters
 Non‐spherical particles have a large influence on the solid flow

For the pie‐slice model of the scaled‐down blast furnace:

 The pie‐slice model generates similar layer structures to the slot model from the
previous simulations
 Increasingly non‐spherical particles generate a larger deadman
 Layers consisting of non‐spherical particles are more resistant to deformation
than spherical particle layers.
 Deadman size is strongly influenced by the front‐ and back wall boundary
conditions

For geometry study:

Compared to half‐diameter and pie‐slice geometries, the full‐diameter slot
geometry is preferred.

For modelling the experimental blast furnace:

 For the EBF a slot model, a geometry of the full width of the actual furnace was
used, with periodic front and back walls.
 Increasing gas velocity from raceway causes layers to increasingly curve upwards
at the centre line of the furnace and cause the pellet layers to split
 Gas flow through the centre is highest due to the lower porosity of the pellet
layers
 Increasing the gas velocity reduces the compressive force networks in the hearth
of the furnace
 Above the cohesive zone the layer structure is similar to the LKAB EBF, but at
the cohesive zone level they differ due to the simplified cohesive zone in the
simulations.

For charging and burden descent models:

 Gas flow loosens the particle packing causing smaller angles of repose. Particle
shape and density have a significant influence on the resulting layer structure
after charging

117
  Wall accretions cause the layer structure to curve upwards directly above them
and layers are mixed below the accretions. Wear of the furnace lining causes the
layer structure to skew downwards and mixing occurs locally

5.2. Recommendations
Though the results from these specific simulations are not directly applicable to industrial
practice, they show clearly important flow characteristics of the burden and gas in the blast
furnaces. To accurately simulate burden descent is very complex; various factors are not
yet taken into consideration in the present simulations. However, a general model of
burden flow can increase the understanding of the processes taking place in the blast
furnace. Although the blast furnace process has been used for a long time, it is often
unclear what exactly takes place inside the furnace. Blast furnace dig‐outs have contributed
to the available knowledge; however, these are undertaken on a quenched furnace and
thus do not reflect precisely the conditions during operation.

From the smaller Experimental Blast Furnace (EBF) of LKAB more data is available on
internal conditions thanks to frequent dig‐outs. Data on the location, shape and size of the
cohesive zone are used to pre‐define the cohesive zone in this thesis research. For a further
developed model which can simulate the formation of the cohesive zone, the data from the
EBF can be used for validation of the model. An example is a study on the influence of
particle size distribution of the pellets on the operating conditions of the EBF. The results
show a clear improvement of the process efficiency and burden descent stability when a
narrow distribution of small particles is used. For a model which can successfully simulate
the reaction kinetics, thermal equilibriums, softening and melting as well as the gas and
solid flow, these experiments are excellent for validation.

Realistic DEM simulations of materials charging would require more accurate geometries,
experimentally derived particles properties and validation. Above the cohesive zone, where
no softening or melting is present, DEM has a great potential for realistic simulations when
the previously mentioned requirements are met. For models including the cohesive zone
the simulations are currently severely limited by the computational demand of the
interconnected sub‐models required to calculate the thermodynamic equilibriums,
softening and melting characteristics, and reaction kinetics.
The coupled gas‐solid model could be developed in steps to include more detail
by using simplified sub‐models. This increases the simulation capability of the models at
reduced computational cost compared to the very demanding detailed sub‐models. The
first such step can be a model which includes a hot blast and a predefined energy sink.
There are many problems with such a model, but it shows the large influence it has on the
gas velocity. One problem is the much faster particle flow compared to reality in order to
decrease the simulation time. There is no speed up applied to the gas flow as it would make

118
the drag on the particles unrealistic. This complicates heat transfer between the gas and
solid flow.

This thesis is part of a larger project with the aim to simulate the blast furnace cohesive
zone. In parallel work sub‐models for thermodynamic equilibrium and reaction kinetics
have been developed. With these models the chemical reactions and the softening and
melting of particles can be simulated. Coupled with the CFD‐DEM model this results in a
model with the capacity of determining the shape and location of the cohesive zone.
However, the simulations are so computationally demanding that only a very short process
time can be simulated, insufficient for the simulation to reach a steady state.
Even though there are issues with the current model as it is, it has succeeded in
coupling the sub‐models with the DEM‐CFD model. This allows for accurate solid flow
modelling in combination with kinetic, thermodynamic and softening and melting models.
It is the first time all the physical and chemical phenomena are included in a comprehensive
DEM‐CFD blast furnace model. Extremely high computational demand of the simulation
severely limits the model at this stage, but the methodology development was a great
success.

At this point there are a number of possible paths to continue the research on DEM
modelling of the blast furnace. The focus can be shifted from the full burden including the
cohesive zone to the burden charging, including stock house, top bunkers and the upper
burden layers. This can be done using only the DEM model coupled with the CFD gas flow,
without the need for sub‐models. Because these simulations do require a large number of
particles, the scale remains a problem. Results can be more directly applicable to the real
industrial blast furnace because no process simplifications are needed, but this does require
more information about of the particle properties e.g. shape, size distribution and friction
factors.

Because the current DEM‐CFD model with the coupled sub‐models is too computationally
demanding, there are two options for further development. The first is to incrementally
develop the DEM‐CFD coupled model with simplified sub‐models which are less
computationally demanding than the currently developed sub‐models. This model in its
different stages of development can be applied to undertake meaningful parameter
studies, e.g. particle size distribution affecting layer porosity or the effect of unsymmetrical
cohesive zone on particle and gas flow.
The second option is further development of the full coupled model to solve the
high computational demand. A large part of this demand is due to the constant switching
and transfer of data between the different sub‐models. If this coupling can be optimized
the simulation will probably be considerably faster, but will remain very time consuming.
The ever increasing computational power and the improvement of DEM software will make
larger and more realistic simulations possible. It is a matter of time before the full coupled
model can be used for large simulations.

119
120
Appendix A
Parallel development of sub‐models
The results presented in this thesis on the combined solid and gas flow burden model. The
goal of the total project is to create a model with the ability to simulate the location and
shape of the cohesive zone. In order to accurately simulate the cohesive zone additional
models are required besides the solid and gas flow models:
 Chemical reactions taking place in the blast furnace require a kinetic model
 Softening and melting can be simulated using a thermodynamic model based on
the chemical composition and its changes of the burden

Figure A.1 shows the project framework discussed in section 1.2, in which the softening and
melt formation model forms the second part of the final model, developed by two post‐
doctoral researchers Yuko Enqvist and Vilas Tathavadkar. This work has focussed on two
parts: the softening and melting behaviour of pellets and sinter, and the reduction of the
ferrous burden by the reduction gas. An additional sub‐model is required to determine the
thermodynamic equilibrium of the burden composition which determines the liquidus or
melting temperature of the pellets and sinter.

The combination of these models creates a tool with which comprehensive simulations of
the blast furnace can be done. Though the DEM‐CFD coupled model can generate
interesting results and is a powerful tool in itself, the real benefit of this project lies in the
coupling with the sub‐models. Combined it can simulate much more of the process in the
blast furnace, most importantly the cohesive zone which has a very large influence. It
therefore requires less pre‐defined data and thus has a better descriptive or even predictive
ability.

Reaction kinetic model

The reaction kinetic model is based on the literature available on this subject. At this stage it
contains the indirect reduction of the pellets using the shrinking core model [2] and the
solution loss reaction (the formation of CO from CO2 and coke) using a model proposed by
Yagi [3]. The reactions in the raceway and the presence of water and hydrogen were not
included in the model.

Softening and melting model

In the softening and melting model the behaviour of the pellets in the cohesive zone is
described. This behaviour is very complex a dependent on a large number of variables.
Softening starts with melt formation inside the micro‐structure of the particles. This hinders
reduction of the ferrous oxides by covering them, so called reduction retardation. The melt
formation reduces the strength of the particle and deformation will take place. Reduction of
the oxides to iron from the outside creates an iron shell which determines the amount of

121
deformation of the particle. The shell is filled with solid ferrous material and an increasing
amount of liquid slag. At a certain point the shell cracks and the ferrous oxide melt will
exude from the particle and get reduced by the gas. The formed iron will remain solid due to
its higher melting point and fill up the pores between the particles. At approximately 1400
C the iron will also start melting.

Blast Furnace Cohesive Zone Model Post Doc.

DEM + CFD + Thermodynamic + Kinetic

Burden Flow Model Cohesive Zone Softening and Melt Formation Model
DEM + CFD Thermodynamic + Kinetic

PhD

Solid Flow Gas Flow Thermodynamic Kinetic Model Iron

Model DEM Model CFD Melt Formation Ore Reduction

Figure A.1. Project framework

This behaviour depends on reduction degree, pellet porosity, basicity, mineralogy,

morphology, phase chemistry and distribution as well as temperature and load. Initially a
model of the softening of pellets is developed. It is assumed that the pellets consists fully of
Fe2O3 and is reduced according to the shrinking core model. For the strength of the pellets a
study using Finite Element Modelling, FEM, has been done in which the deformation of an
iron shell is investigated at different temperatures. This results in a correlation between
deformation, temperature, reduction degree and load. This correlation does not yet agree
with experimental results very well and further work is needed. The model describing the
softening and melting is still in development and requires a lot of additional parameters to
describe the process accurately.

Experiments have been done to study the softening and melting behaviour and supply
additional data by Dr. Tathavadkar. The material was quenched to investigate the slag
composition and microstructure. The experimental equipment is shown in Figure A.2, where
a 0.09 m high sample can be heated while being reduced and under pressure. During the
heating and reduction of the sample the pressure drop and the reduction of the bed height
are measured, shown in Figure A.3 with a sample taken from the furnace.

122
 Figure A.2. Experimental equipment at Tata Steel, Jamshedpur

Figure A.3. Pressure drop and bed shrinkage as a function of temperature (by Dr.
Tathavadkar)

123
 Figure A.4. Liquid solid ratio as a function of
temperature and reduction degree

According to Tathavadkar bed shrinkage starts at approximately 800 C and is 50% at 1200
C, which is defined as the softening temperature Ts. At this temperature the pressure drop
over the sample increases indicating liquids filling the void between the particles. This gives
an idea of how these specific pellets behave in the blast furnace. The softening and melting
process of pellets and sinter are very different and further work has been done in which the
liquid to solid ratio of pellets and sinter as a function of temperature is determined at
different reduction degrees, shown in Figure A.4. At lower reduction degrees there are more
oxides present, which have a lower melting point and thus the liquid/solid ratio is higher.
Sinter contains more slag materials than pellets and therefore also generates more liquid
slag.
The thermodynamic equilibrium model is based on the thermodynamic software
FactSage (GTT Technologies) with which the equilibria are calculated at different reduction
degrees and temperatures by Dr. Enqvist. For the coupling with the CFD model ChemApp is
used which allows the FactSage data to be used for thermodynamic calculations in the CFD
cell.

124
 EDEM ‐ FLUENT model

Custom coupling module

Figure A.5. Model integration [4]

The coupled blast furnace cohesive zone model [4]

A very challenging part of the work consisted of combining all these models. The reaction
kinetic and thermodynamic equilibrium models were applied in the CFD model and the
softening and melting model in the DEM model. Figure A.5 shows the integration of the
models. In addition to the EDEM‐Fluent coupling module a custom coupling module is
required for additional data transfer between the CFD and DEM models not allowed by the
commercial coupling module. This custom coupling module was originally developed by
Cristoph Kloss from the Johannes Kepler University in Linz [5, 6]. The simulations are based
on the LKAB EBF on which the simulations in this thesis are also based. Periodic boundary
conditions are used for both front and back walls. Two sets of initial conditions for the
temperature profile were used for the simulations: one based on a reversed V shaped
cohesive zone as observed in the EBF, and one based on steady‐state calculation results
from MOGADOR. MOGADOR is a 2D CFD model which can simulate and predict the
internal state of the blast furnace [7]. In both simulations the iron ore particles are removed
at 1573 C, the melting point of iron.

125
 (a) (b) (c)
Figure A.6. CFD Simulation results[4]: (a) gas velocity, m/s;
(b) gas temperature, K; (c) CO2 mole fraction

Total simulated time was 15 s, and when the gas flow stabilised after 5 s the thermodynamic
equilibrium and softening and melting models were included. Figures A.6 and A.7 show the
coupled calculation results of the CFD and DEM simulation, with the MOGADOR initial
temperature profile. These show that the coupling is working successfully; the gas is cooling
due to the burden heat up and reactions, and the CO2 fraction increases in the ore layers
due to the reduction of iron ore and decreases in the coke layers where it is reverted to CO.
The CO2 fraction is similar to MOGADOR results. For the softening and melting model in the
DEM simulation the deformation degree is successfully integrated. It is calculated from the
softening viscosity and the load, the former is determined by the thermodynamic
equilibrium model.

126
 (a) (b)
Figure A.7. DEM Simulation results: (a) Reduction degree;
(b) deformation degree

The main issue with the model at this stage is the extremely heavy computational load
which allows for only a short time to be simulated. This is far too short for the simulation to
establish a steady‐state and the initial conditions have a large influence. The evolution of
the softening and melting in the cohesive zone needs more time to be simulated, but to
simulate 1 hour of physical time would require 10 years of calculation time. Besides the
computational load the stability of the solution was a problem; unpredictable errors and
numerical instabilities caused frequent crashes. For the current software to be successfully
applied to simulation of the cohesive zone it requires greater stability and a very significant
speed‐up.

However, this is the first time all these physical and chemical phenomena are included in a
comprehensive DEM‐CFD blast furnace model. This allows for accurate solid flow modelling
in combination with kinetic, thermodynamic and softening and melting models, only
limited by computing power and efficiency, and the refinement of the sub‐models

127
References
1. Hooey, L., Sköld, B.‐E., Sundqvist Ökvist, L., Seppänen, M., and Zuo, G. LKAB's
experimental blast furnace ‐ the learning curve. in The 5th European Coke and
Ironmaking Congress. (2005), Stockholm, Sweden.
2. Omori, Y., ed. Blast furnace phenomena and modelling. (1987), Elsevier Applied
Science.
3. Yagi, J., Sasaki, K., and Muchi, I., Theoretical investigation on the blast furnace
operations with the aid of mathematical model. Tetsu to Hagane 54 (1968), p. 1019‐
1031.
4. Enqvist, Y., Modelling of the blast furnace cohesive zone using a coupled DEM‐CFD
approach, in M2i Technical Report. (2012).
5. Kloss, C., Goniva, C., Aichinger, G., and Pirker, S., Comprehensive DEM‐DPM‐CFD
simulations ‐ model synthesis, experimental validation and scalability, in Seventh
international conference on CFD in the minerals and process industries. (2009):
Melbourne, Australia.
6. Kloss, C., Aichinger, G., and Pirker, S., Multi‐scale modelling of particle motion by
means of DEM and DPM, in Multi‐scale modeling symposium (2009): Melbourne,
Australia.
7. Danloy, G., Mignon, J., Munnix, R., Dauwels, G., and Bonte, L. A blast furnace
model to optimize the burden distribution. in ICSTI Ironmaking Conference. (2001).

128
Summary

The focus in the steel industry is on increasing production and efficiency of the existing blast
furnaces. The driver for this development is economic as well as environmental; reduction
of greenhouse gas emissions has become a critical issue. Within the European Union several
large projects have been running to study CO2 reduction in adapted blast furnaces, e.g. top
gas recycling, or alternative technologies, e.g. the HISARNA smelting reduction process.

In order to improve the process or to develop new technologies for the blast furnace, a good
understanding of the processes inside the blast furnace is essential. One of the most
important parameters determining the process is the cohesive zone where the softening
and melting of the ferrous materials take place. Due to its low permeability and cohesion it
has a very large influence on the gas flow resistance and thus also on the descending solids.

Details on the cohesive zone are very difficult to measure and predict. Measurements on the
state of the process inside the blast furnace can be done indirectly on the outgoing flows,
e.g. top gas or hot metal temperature and composition. Another way is by probing into the
furnace, which is very difficult to do on a regular basis. Prediction of the cohesive zone
behaviour is very difficult due to the complexity of its behaviour and many influencing
factors. Some knowledge on the cohesive zone has been gained by rare dig‐outs of
quenched blast furnaces. These give an indication of its shape, location and composition.
But they provide only a snapshot and do not give any information on how the cohesive zone
reacts to changing process parameters.

A tool that can improve this understanding is mathematical modelling. This thesis describes
the application of a Discrete Element Method (DEM) model coupled with a Computational
Fluid Dynamics (CFD) model to study solid and gas flows in the blast furnace. In parallel,
work was done to develop softening and melting models, and reaction kinetic models for
burden materials. The full project, Prediction of physical and chemical properties of the
burden materials in the cohesive zone, was conducted with the financial support of the
Materials innovation institute (M2i), in cooperation with Tata Steel Europe.

The use of DEM modelling for granular flows has increased greatly in recent years due to
the increase in computing power. For the application of DEM models in blast furnace
ironmaking, three types of simulations can be identified: the burden flow, the raceway and
burden charging. At the current level of development of DEM modelling, it can be most
readily applied to the third type, the burden charging. For both raceway and burden flow
the complexity is much greater and requires additional sub‐models.

The results of the present study show that the combined DEM‐CFD method can be
successfully applied to the simulation of the burden flow in the blast furnace. The burden
porosity distribution as created by the DEM model influences the gas flow, which in turn
influences the solid flow and the layer structure. The simplified geometry of the furnace

129
used in the simulation is shown to have a large influence on the solid flow and on the
resulting deadman shape and size. The presence of a deadman above the tuyeres was found
to be mostly caused by the geometry, rather than the solid flow. Simulations show that the
optimal simplified geometry closest to a cylindrical model is a full diameter slot model.
In blast furnace DEM simulations, spherical particles are generally used. It is
shown here that non‐spherical particles significantly influence the solid flow. Layers of non‐
spherical particles are more resistant to deformation due to interlocking and resistance to
rotation. The initial layer structure as charged remains more constant while descending
through the furnace.
The non‐spherical coke particles create layers with a higher porosity than the
spherical pellets. The centre line of the furnace is also filled with coke, as is common blast
furnace practice. This creates a higher gas flow through the centre, resulting in locally
higher particle drag on the pellets, causing the pellet layers to curve upwards towards the
centre with increasing gas velocity. Simulations illustrate the influence of the gas flow on
the force networks in the burden and how it interacts with the furnace wall. The
compressive force networks in the hearth become smaller with increasing gas velocity.
In simulations of burden charging both the gas flow and the particle shape and
density have a large influence on the resulting layer shape. Accretions on the blast furnace
wall cause the layers to curve upwards and cause particle mixing immediately below.

The parallel work which has been done in the full project, Prediction of physical and chemical
properties of the burden materials in the cohesive zone, links the gas and solid flows with
softening and melting, reaction kinetics and thermodynamic equilibrium within the burden.
What severely limits the model at this stage is the very high computational power it
requires. However, the larger model with the sub‐models includes all the data required to
successfully simulate the cohesive zone based on the particle composition, softening and
melting, chemical reactions, and heat transfer as well as gas and solid flow.

130
Samenvatting
De focus in de staalindustrie ligt op de verhoging van de productie en de efficiency van
bestaande hoogovens. Deze ontwikkelingen worden gedreven door zowel economische als
ecologische factoren; het terugdringen van de uitstoot van ozongassen is de laatste jaren
een belangrijk onderwerp geworden. Binnen de Europese Unie zijn enkele grote
onderzoeksprojecten gaande naar de reductie van CO2 uitstoot van hoogovens,
bijvoorbeeld door top gas recycling, of alternatieve technologieën zoals het HISARNA
proces.

Om hoogovenprocessen te verbeteren of nieuwe technologieën te ontwikkelen is een goed

begrip van de processen die binnenin de hoogoven plaatsvinden essentieel. Een van de
belangrijkste parameters die bepalend is voor het hoogovenproces is de cohesieve zone,
waar het verweken en smelten van de ijzerhoudende grondstoffen plaatsvindt. Door de lage
permeabiliteit en cohesie heeft de cohesieve zone een grote invloed op de weerstand van
de gasstroming en daarmee ook op het zakken van de belading.

Details van de cohesieve zone zijn erg moeilijk te meten en voorspellen. Metingen van de
procescondities binnenin de oven kunnen alleen indirect plaatsvinden op de uitgaande
stromen, bijvoorbeeld top gas of ruwijzer temperatuur en samenstelling. Een andere wijze
is sonderen, hoewel dit erg moelijk is om op regelmatige basis uit te voeren. Het
voorspellen van het gedrag van de cohesieve zone is erg moeilijk door de hoge complexiteit
van het gedrag en de grote hoeveelheid factoren die hierop van invloed zijn. Enige kennis
van de cohesieve zone is vergaard bij vrij zeldzame uitgravingen van afgekoelde hoogovens.
Deze geven een indicatie van de vorm, lokatie en samenstelling van de cohesieve zone,
maar geven geen informatie over het dynamische gedrag van de cohesieve zone tijdens het
hoogovenproces.

Wiskundig modelleren is een hulpmiddel dat tot een beter begrip van de cohesive zone kan
leiden. Dit proefschrift beschrijft de toepassing van de Discrete Element Method (DEM)
gekoppeld aan Computational Fluid Dynamics (CFD) om de vaste stof‐ en gasstromingen in
de hoogoven te bestuderen. Parallel hieraan zijn modellen ontwikkeld om het verweek‐ en
smeltgedrag en de reactiekinetiek te kunnen beschrijven. Het volledige project, Prediction
of physical and chemical properties of the burden materials in the cohesive zone, is uitgevoerd
met financiele ondersteuning van het Materials innovation insitute (M2i) in samenwerking
met Tata Steel Europe.

Het gebruik van DEM modellering voor deeltjesstroming is de laatste jaren enorm
toegenomen door de vooruitgang in rekenkracht van computers. Toepassingen bij
hoogovens concentreren zich op drie gebieden: de deeltjesstroming in de oven, de raceway
en het beladen van de oven. Met het huidige niveau van ontwikkeling van de DEM methode

131
is de meest praktische toepassing de laatste. Voor de raceway en de stroming in de oven is
de complexiteit veel groter en zijn extra sub‐modellen nodig.

Het resultaat van deze studie laat zien dat de gekoppelde DEM‐CFD methode succesvol kan
worden toegepast op de simulatie van deeltjesstroming in de hoogoven. De
porositeitsverdeling van de belading zoals die door het DEM model is berekend beïnvloedt
de gasstroming, die vervolgens weer de deeltjesstroming en de lagenstructuur beïnvloedt.
De versimpelde geometrie die gebruikt is voor de simulaties heeft een grote invloed op de
deeltjesstroming en de uiteindelijke vorm en grootte van de dodeman. De aanwezigheid
van een dodeman boven het tuyereniveau bleek voornamelijk te zijn veroorzaakt door de
geometrie en niet zozeer door de deeltjesstroming. Simulaties laten zien dat de optimale
geometrie, die qua resultaat het meest op een cylindrisch model lijkt, een platte doorsnede
is van de volledige diameter van de oven.
In hoogoven DEM‐simulaties worden over het algemeen perfect ronde deeltjes
gebruikt. Het wordt aangetoond dat het gebruik van niet ronde deeltjes de resultaten
significant beïnvloedt. Lagen van niet‐ronde deeltjes zijn beter bestand tegen deformatie
doordat deeltjes in elkaar grijpen en grotere weerstand tegen rotatie hebben. De
lagenstructuur zoals die wordt gestort blijft beter behouden tijdens het zakken in de oven.
De niet‐ronde kooksdeeltjes creëren lagen met een grotere porositeit dan de
ronde pellets. Het hart van de oven is van boven tot beneden gevuld met kooks, wat bij
hoogovens een gebruikelijke procesvoering is. Hierdoor is er een grotere gasstroming door
het hart, wat resulteert in lokaal hogere krachten op de pellets waardoor de pelletlagen
naar boven buigen richting het hart. Simulaties laten ook de invloed van de gasstroming op
de krachtennetwerk in de oven zien en wat de interacties hiervan zijn met de wand. De
krachten op de deeltjes in de haard worden kleiner naarmate de gassnelheid toeneemt.
In simulaties van het storten van de belading hebben zowel de gassnelheid als de
vorm en grootte van de deeltjes een grote invloed op de uiteindelijke lagenstructuur.
Aanzettingen aan de wand zorgen ervoor dat lagen naar boven buigen en dat de deeltjes
eronder zich vermengen.

Het werk wat gelijktijdig werd uitgevoerd in het kader van het bredere project, Prediction of
physical and chemical properties of the burden materials in the cohesive zone, verbindt de gas‐
en deeltjesstroming met verweken en smelten, reactie kinetiek en thermodynamische
evenwichten binnen de belading. Wat dit model sterk beperkt in zijn toepassing zijn de erg
hoge eisen die het stelt aan computervermogen. Desalniettemin is het gelukt om het grote
model te verbinden met de sub‐modellen, inclusief alle data die nodig is om de cohesieve
zone te kunnen simuleren.

132
List of Publications

Adema A.T., Yang Y., Boom R., Discrete Element Method ‐ Computational Fluid Dynamic
simulation of the materials flow in an iron‐making blast furnace, ISIJ International, 50, (2010)
p. 954‐961

Adema A.T., Yang Y., Boom R., Coupled DEM ‐ CFD Modelling of the Iron‐Making Blast
Furnace, Journal of Iron and Steel Research, International, 16 (Part 2, Suppl. 2), (2009), p.
1081‐1085.

Adema A.T., Yang Y., Boom R., Coupled DEM‐CFD modelling of the iron‐making blast
furnace, Proceedings of the 5th International Congress on the Science and Technology of
Ironmaking, ICTSI'09, (2009), Shanghai, China, The Chinese Society for Metals, pp. 1085‐
1089

Adema A.T., Yang Y., Boom R., Coupled DEM‐CFD modelling of the iron‐making blast
furnace, Proceedings of the 7th International Conference on CFD in the Minerals and
Process Industries, (2009), Melbourne, Australia, CSIRO

Adema A.T., Yang Y., Boom R., DEM Simulation of solid flow in an ironmaking blast furnace –
influence of non‐spherical particles, International symposium on Recent progress on
mathematical modeling in ironmaking 2008, (2008) October 16‐17, Tokyo, Japan, Iron and
Steel Institute of Japan, pp. 130‐135

Adema A.T., Yang Y., Boom R., DEM Simulation of solid flow in an ironmaking blast furnace –
influence of non‐spherical particles, Proceedings of the XXIV International Minerals
Processing Congress, IMPC 2008, (2008) 24‐28 September, Beijing, p. 2422‐2429

Adema A.T., Yang Y., Boom R., Prediction of the physical and chemical properties of the
burden materials in the cohesive zone of the ironmaking blast furnace, Proceedings of Discrete
Element Methods 07, DEM 07, (2007) August 27‐29, Brisbane, Australia

Adema A.T., Yang Y., Boom R., Prediction of the cohesive zone properties in the ironmaking
nd
blast furnace, Proceedings of the 2 International Conference of Simulation and Modelling
of Metallurgical Processes in Steelmaking, Steelsim 2007, (2007) September 12‐14,
Conference Centre Graz/Seggau, Austria

133
134
Acknowledgements

I would like to express my very great appreciation to Yongxiang Yang, who has been my
daily supervisor during my PhD research work. But before that has also supervised my MSc
thesis work and helped me arrange an internship at the Anhui University of Technology in
Ma’anshan, China. Thank you for your support and advice during all these years.

My promoter Rob Boom has been an inspiration, our discussions always gave me new
energy and ideas on the project at times when the going was tough, and helped me
overcome any difficulties.

The support and discussions with Jan van der Stel and Mark Hattink at the industrial partner
Tata Steel R&D have been invaluable for this project. I would also like to thank Luc Bol,
Teun Bakker and Hans Jak from the Blast Furnace Department for giving me the
opportunity to finish writing the thesis after joining Tata.

For their friendship I would like to thank my friends and colleagues at the TU Delft; Christa,
Abbas, Floris, Tanya, Ali, Demian, Zhan, Yuko, Liang, Andrey, Tungky, Liliang, Bogdan,
Zhiguang, Yulia, Lambert, Sander, Andrea, Jun, Qingshi, Yingxia. I have greatly enjoyed
your company during much needed coffee breaks and lunches during my PhD days in Delft,
and also outside the university.

Finally I would like to thank my family for the support and the faith they have always had in
me. Without them I would have never gotten to this point.

Last but not least I would like to thank Sylvia for the patience, encouragement and love she
gave me.

135
Curriculum Vitae
Personal details
Allert Tjipke Adema
Surhuizum, 27‐07‐1979

Work experience
2011‐present Process Technologist, Blast Furnace department,
Tata Steel Europe
2006, Feb‐July Internship at Rössing Uranium Ltd., Namibia, Graduate Metallurgist for
Rio Tinto Technical Services
2003, Feb‐April Research internship at Anhui University of Technology, China

Education
2006‐2010 PhD Researcher for the Materials innovation institute (M2i) at Delft
University of Technology
2003‐2004 European Mineral Engineering Course (EMEC),
Combined curriculum at Delft University of Technology, RWTH Aachen,
Helsinki University of Technology and Imperial College London
1998‐2006 Delft University of Technology, Applied Earth Sciences,
MSc in Resource Engineering

Honors and awards

2008 Young Achievers Award, Presented at the XXIV International Mineral
Processing Congress, September 27, Beijing, China, by the International
Mineral Processing Council (IMPC)

Extra‐curricular
2007‐2013 Secretary of the KIVI NIRIA (Royal Institute of Engineers), section Mining

136
I
SBN:
978-
94-
91909-
06-
1