Calculation of Coke Layers Situation in the Cohesive Zone of Blast Furnace

Coke is the only batch component that does not soften in blast furnace thermal conditions. It is especially important at the temperatures of the cohesive zone forming because coke layers are the only gas-permeable charge. The aim of the work described in this article is the identification of individual coke layers situation in the cohesive zone. Numerical calculations of the cohesive zone situation are based on the horizontal below burden probe measures, however, coke layers are calculated using analytical geometry. The results can be presented as a bitmap; the individual and total area of the coke layers passing gases through the cohesive zone is also calculated. This form of results allows for subjective but quick assessment of the blast furnace operation by its crew.


Introduction
Hot metal (HM) production is the most expensive step in the steel production cycle. Although the efficiency of blast furnace (BF) mostly depends on the charge materials quality [1,2], mathematical modeling of partial phenomena or the whole process can help to understand and highlight current shortcomings in BF technology [3,4].
The cohesive zone (CZ) is one of the crucial issues related to the modeling of the BF process. It is a bottleneck for the flow of reducing gases in the BF shaft, and its location and shape determine the correctness and continuity of the hot metal production process.
According to the definition, the cohesive zone is an area closed between the isotherms of the softening and melting of the charge oxide materials (sinter, pellets, flux, lump ore), which, being in a soft state, stick to each other and thus are impermeable to gases [5]. The only places permitting gas flow are the coke slits between cohesive ore charge layers ( Figure 1). Thus, the cohesive zone plays the role of a gas distributor in the charge materials solid-state area [6].
All the mentioned models have a common goal, namely, to describe phenomena that cannot be measured or observed directly due to the working conditions of a sealed object, which a blast furnace is. However, cold models contribute to an understanding of macro-scale particle motion phenomena, while on the other hand, DEM-CFD models help to visualize the interaction between the different phases on a micro-scale. Numerous experimental and numerical studies devoted to charging particle behavior or identifying the shape and location of the cohesion zone can be found elsewhere. Among them should be mentioned cold models [7][8][9] or mathematical models [10][11][12][13][14][15]. Also, more advanced modeling techniques like the discrete element method (DEM) [9,[16][17][18], computational fluid dynamics (CFD) models [19][20][21][22][23] and a combination of the two [24][25][26][27][28] are also reported.
All the mentioned models have a common goal, namely, to describe phenomena that cannot be measured or observed directly due to the working conditions of a sealed object, which a blast furnace is. However, cold models contribute to an understanding of macro-scale particle motion phenomena, while on the other hand, DEM-CFD models help to visualize the interaction between the different phases on a micro-scale.
The miscellaneous types of input data are used as an initial condition. Most of them are measured outside the burden by utilities like thermocouples in the wall, above burden probes, radar profilometer, or throat IR camera [13,15,29,30]; cooling water temperature may also be taken into account. However, there are measuring methods known to be more intrusive into the BF working space such as descending probe (DP) and horizontal below burden gas probe (BBP).
DP is a disposable pipe that dropped on the burden surface and allows for continuous measuring of gas composition and temperature 3-5 points from the wall to the axis of BF. Although the method is accurate, it takes a long time to reach the cohesive zone level (two-three hours), so any possible response to irregularities could be delayed. On the other hand, DP helped to learn and understand the phenomena in the BF shaft.
However, BBPs became more popular and have been widely used since the 1980s. Usually, BBP is situated 3-6 m below the burden surface and measures the gas composition and its temperature at 6-10 points (depending on the BF size). A measure lasts 3-4 min per point.
It can be seen that the mentioned models  differ from each other not only in the method used but mostly in the type and place of input data collection. The less in- The miscellaneous types of input data are used as an initial condition. Most of them are measured outside the burden by utilities like thermocouples in the wall, above burden probes, radar profilometer, or throat IR camera [13,15,29,30]; cooling water temperature may also be taken into account. However, there are measuring methods known to be more intrusive into the BF working space such as descending probe (DP) and horizontal below burden gas probe (BBP).
DP is a disposable pipe that dropped on the burden surface and allows for continuous measuring of gas composition and temperature 3-5 points from the wall to the axis of BF. Although the method is accurate, it takes a long time to reach the cohesive zone level (two-three hours), so any possible response to irregularities could be delayed. On the other hand, DP helped to learn and understand the phenomena in the BF shaft.
However, BBPs became more popular and have been widely used since the 1980s. Usually, BBP is situated 3-6 m below the burden surface and measures the gas composition and its temperature at 6-10 points (depending on the BF size). A measure lasts 3-4 min per point.
It can be seen that the mentioned models  differ from each other not only in the method used but mostly in the type and place of input data collection. The less intrusive the measure, the more sophisticated the modeling methods used. However, less complex models should be sufficient for quick assessment of BF operation, but they must be based on measurements taken as close as possible to the place where a specific phenomenon occurs.
To sum up, the selection of the modeling method for a specific phenomenon is based on the type and availability of input data, but also should depend on the speed of calculations, permissible error, and ease of implementation.
The current paper presents a mathematical model based on BBP measures. The Rist principles, namely the mass and heat transfer consideration in three separate zones, preparation, reserve, and producing, are taken as fundamentals. The work data of the blast furnace No.5 in the Krakow steel plant was used for the calculations. The current model has been gradually developed for almost four decades. Invented by Itaya et al. [10], adopted to the conditions of Polish blast furnaces by Benesch et al. [31], upgraded by Stachura et al. [32], was finally improved by us and implemented in 2018 on the blast furnace No. 5 in Krakow. Here, the latest form is presented.

Object of Research
The referenced blast furnace has an overall volume of 2000 m 3 and production ability is 1.2-1.4 million tons of hot metal annually. Hearth diameter is 9.75 m, along the perimeter of which is placed 24 tuyeres. The usual hot blast (HB) volume is about 190,000 m 3 per hour.
BF No.5 is equipped with Dango & Dienenthal BBP situated 4.9 m below the throat and measures gas composition and temperature in 8 points, the distance between which is 0.54 m. In accordance with BBP geometry, a working space of BF is divided into 8 concentric segments ( Figure 2). Therefore, we decided to define the levels of softening and melting isotherms occurrence for each segment and then to draw isotherms for the entire furnace. Next, it will be possible to calculate the coke slits situation. To achieve this goal, the work may be divided into the following main stages:

1.
Determination of the charge materials mass and blast for each of the eight segments and their balance with the parameters of the entire furnace; 2.
Determination of the preparation zone operating parameters and the levels of the 1000 • C isotherm as a thermal and chemical reserve zone; 3.
Determination of the producing zone parameters using the 1000 • C isotherm level as input conditions; 4.
Drawing of cohesion zone, coke layers, and calculating of their area.
upgraded by Stachura et al. [32], was finally improved by us and imp the blast furnace No. 5 in Krakow. Here, the latest form is presented.

Object of Research
The referenced blast furnace has an overall volume of 2000 m 3 a ity is 1.2-1.4 million tons of hot metal annually. Hearth diameter is 9 rimeter of which is placed 24 tuyeres. The usual hot blast (HB) volu m 3 per hour.
BF No.5 is equipped with Dango & Dienenthal BBP situated 4.9 and measures gas composition and temperature in 8 points, the dist is 0.54 m. In accordance with BBP geometry, a working space of B concentric segments ( Figure 2). Therefore, we decided to define the and melting isotherms occurrence for each segment and then to dra entire furnace. Next, it will be possible to calculate the coke slits situa goal, the work may be divided into the following main stages: 1. Determination of the charge materials mass and blast f segments and their balance with the parameters of the e 2. Determination of the preparation zone operating param of the 1000 °C isotherm as a thermal and chemical reserv 3. Determination of the producing zone parameters using t level as input conditions; 4. Drawing of cohesion zone, coke layers, and calculating

Radial Material Balance
Radial distribution of materials charged from the top can be found using gas analyses of BBP measures. For this purpose, a system of linear equations describing the balance of elements or compounds such as nitrogen (1), hydrogen (2), carbon (3) and oxygen (4) was compiled: where: (k)-index of one of eight concentric segments; CO 2 , CO, N 2 , H 2 -components of the gas measured on BBP level, in %; G N2tuy , G Ctuy , G H2tuy , G Otuy -masses of basic elements or compounds at the tuyeres level, in kg/m 3 of the hot blast(HB); V gas -a volume of gas at the BBP level, in m 3 /m 3 HB; C dr -a mass of carbon used in solution loss reaction (direct reduction of wustite), in kg/m 3 HB; P HM -hot metal production, in kg/m 3 HB; W-a mass of H 2 O produced in reduction by hydrogen, in kg/m 3 HB; G CHM -a mass of carbon solved in hot metal, in kg/kgHM; G CSiMnP -a mass of carbon used for the reduction of Mn, Si, and P, kg/kgHM; N 2coke , C coke , H 2coke -contents of nitrogen, carbon, and hydrogen in coke, in mass%; G OFe -a mass of oxygen coming from Fe reduction, in kg/kgHM; G OSiMnP -mass of oxygen coming from Si, Mn and P reduction, in kg/kgHM; Amd Ocoefficient of oxygen amendment resulting from occurring reduction processes over BBP level (will be discussed further). Variables in bold are unknowns. The system of Equations (1)-(4) is resolved for every eight concentric segments (k). Then, it is possible to calculate the share of charge materials masses, in kg/m 3 HB: where: i-means any material charged from the top except coke (i.e., sinter, pellets, etc.); G CCtuy -a mass of coke carbon burnt at the tuyeres level, in kg/m 3 HB; G ichg -a mass of i material in a single charge, in kg/chg; P HMchg -hot metal production from a single charge, in kg/chg. It can be seen that most variables are expressed per cubic meter of the hot blast, so the most important thing now is the calculation of the blast volume for particular segments. The distribution of the blast volume will depend on the geometrical dimension and porosity of each segment. In turn, the porosity index depends on the share of coke and other charge materials.
V HB (k) = V HB tuy · S(k) where: n-max number of input material (except coke); PI-segmental porosity index, in m 3 /chg; PI chg -porosity index of a single charge, in m 3 /chg; V HB -segmental hot blast volume, in m 3 /h; V HBtuy -overall hot blast volume at the tuyeres level, in m 3 /h; S-an area of segment cross-section at the BBP level, in m 2 .
After calculating the blast distribution, the segmental material distribution in kg/chg is calculated and then their sum is compared with the real mass of the single charge. In the first iteration, the difference is high, due to the fact that the calculations were carried out to the BBP level as if reduction processes take place only below the probe. On the other hand, it is known that in such a thick charge layer above the BBP, and at such a high temperature, the mass exchange must take place. For this reason, the carbon and oxygen balance Equations (3) and (4) are amended for oxygen.
where: CO thr ,CO 2 thr and t thr -gas composition and temperature of the gas at the throat level, in vol.%; 1020-the assumed gas temperature at the reserve zone level where charge temperature is 1000 • C (on this level CO 2 does not exist yet and CO volume is the same as the sum of CO thr and CO 2 thr ), in • C; Amd itr -iteration amendment, which during balancing of segmental masses with a real mass of single charge is changed numerically from 1 down to 0 by 0.01 for each iteration and calculations carried out in loop from Equation (1) until misbalance reaches max 2%.

Top Gas Temperature Radial Distribution
The most inconvenient for determination of the top heat exchange zone parameters is the lack of above burden probe (ABP). So, the gas temperature radial distribution at the throat level must be calculated numerically. It may be realized by transferring of the linearized temperature distribution at the BBP level to the throat level ( Figure 3). volume, in m 3 /h;VHBtuy -overall hot blast volume at the tuyeres level, of segment cross-section at the BBP level, in m 2 .
After calculating the blast distribution, the segmental material d is calculated and then their sum is compared with the real mass of t the first iteration, the difference is high, due to the fact that the calcu out to the BBP level as if reduction processes take place only below other hand, it is known that in such a thick charge layer above the BBP temperature, the mass exchange must take place. For this reason, the balance Equations (3) and (4) are amended for oxygen.
where: COthr,CO2 thr and tthr -gas composition and temperature of t level, in vol.%; 1020 -the assumed gas temperature at the reserv charge temperature is 1000 0 C (on this level CO2 does not exist yet an same as the sum of COthr and CO2 thr), in 0 C; Amditr -iteration amend balancing of segmental masses with a real mass of single charge is c from 1 down to 0 by 0.01 for each iteration and calculations carrie Equation (1) until misbalance reaches max 2%.

Top Gas Temperature Radial Distribution
The most inconvenient for determination of the top heat exchan is the lack of above burden probe (ABP). So, the gas temperature radi throat level must be calculated numerically. It may be realized by linearized temperature distribution at the BBP level to the throat leve  Firstly, an average gas temperature at the BBP level (10), linear dependence of the temperature distribution on the radius (11) and reference radius by which the transfer will take place are determined (12).
where: t BBP -the average gas temperature at the BBP level, in • C; R BBP -average weighted by segments squares radius at the BBP level, in m; R-radius of the cylindrical segment corresponding with BBP measuring points (so R(1)=0), in m; A,B-linear regression coefficients; R trsf -reference radius of temperature distribution transferring, in m. Then, the gas temperature distribution (13) and average gas temperature (14) at the throat level are calculated: where, t GAS -the gas temperature distribution at the throat level, in • C; t THR -calculated average gas temperature at the throat level, in • C. If calculated t THR is significantly different from t thr , the A parameter (tangent of the linear distributed temperature angle) is changed numerically.

Reserve Zone Situation
Preparation zone is characterized by following assumptions: 1.
here is no direct reduction reaction of wustite; 2.
Wall heat loss is neglected; 3.
At the lower boundary, equated to the reserve zone, the batch temperature is 1000 • C and gas is 1020 • C.
Thus, heat exchange could be written as a differential equations system: where:-t Z and T Z mean gas and batch temperature at the level Z respectively, in C; ω t and ω ψ T -heat capacities fluxes of gas and batch, in kJ/ • C·h; Z T -level of occurrence of T isotherm, in m; K 1 -volumetric heat transfer coefficient for preparation zone, in kJ/m 3 • C·h. Assuming the initial integration conditions t 0 = t GAS and T 0 = 10 • C (temperature of loading batch), the solutions of the above-mentioned equation system are: where: α-quotient ψ T /ω t which is assumed to be constancy in the whole preparatory zone for each segment k, wherein 0 < α < 1. From Equations (17) and (18) and using known BBP level, the K 1 and α values can be found by inverse method. For this purpose, it is assumed: where: C CO , C CO2 , C H2 , C N2 , C H2O -molar heat capacities of gas components calculated at the temperature of BBP level t gas , in kJ/ • C mol.
Next, with K 1 and α values, the level of 1000 • C is calculated according to Equation (17). Another value of ω t must be taken since at this level the gas temperature is 1020 • C and its composition is also different.

Producing Zone Operating Parameters
Producing zone is characterized by following assumptions (Figure 4):

1.
At the top boundary (corresponding with the reserve zone level) the batch temperature is 1000 • C and gas is 1020 • C; 2.
At the lower boundary the batch temperature is 1450 • C and gas is the same as raceway adiabatic flame temperature (RAFT); 3.
Heat loss intensity is changed with the gas temperature, it achieves max value at the tuyeres level and min at the reserve zone level; 4.
Intensity of direct reduction reaction is changed with the batch temperature according to the normal distribution; it achieves max value at 1225 • C.
where: CCO, CCO2, CH2, CN2, CH2O -molar heat capacities of gas components calculated at the temperature of BBP level tgas, in kJ/°C mol. Next, with K1 and α values, the level of 1000 °C is calculated according to Equation (17). Another value of ωt must be taken since at this level the gas temperature is 1020 °C and its composition is also different.

Producing Zone Operating Parameters
Producing zone is characterized by following assumptions (Figure 4): 1. At the top boundary (corresponding with the reserve zone level) the batch temperature is 1000 °C and gas is 1020 °C; 2. At the lower boundary the batch temperature is 1450 ° C and gas is the same as raceway adiabatic flame temperature (RAFT); 3. Heat loss intensity is changed with the gas temperature, it achieves max value at the tuyeres level and min at the reserve zone level; 4. Intensity of direct reduction reaction is changed with the batch temperature according to the normal distribution; it achieves max value at 1225 °C. In this zone, as a result of the direct FeO reduction reaction, the weight of the charge and gas is significantly decreased. Under the conditions of countercurrent motion, this phenomenon can be written as two differential equations: In this zone, as a result of the direct FeO reduction reaction, the weight of the charge and gas is significantly decreased. Under the conditions of countercurrent motion, this phenomenon can be written as two differential equations: dm(k) dT Z (k) = − 28 12 ·R DR (k) (24) where: m and M are mass fluxes of gas and batch, respectively, in kg/h; R DR is the intensity of carbon consumption by FeO direct reduction reaction in the producing zone at the temperature T Z , in kg/ • C·h. The minus sign means that gas and batch masses decrease downward.
The batch extracts heat from the gas for its heating, melting, and direct reduction, and the gas also covers heat losses. Thus, the heat transfer balance from gas to the charge can be presented in the form of the following differential equations: where: Q HL -heat loss intensity at the t Z temperature, in kJ/ • C·h; H DR -enthalpy of direct reduction reaction per 1 kg of carbon, in kJ/kg; K 2 -volumetric heat transfer coefficient for producing zone, in kJ/m 3 · • C·h. Value of Q HL depends on the gas temperature and max heat loss intensity for each segment (see Figure 4): RAFT−1020 dt z (k) (28) where Q PZ -heat loss for the whole producing zone, in kJ/h. However, the value of R DR depends on the batch temperature and max direct reduction intensity (see Figure 4): Producing zone calculations are conducted numerically starting from the reserve zone level to tuyeres level (25.4 m) with the step of 0.01 m. Firstly, must be found coefficient K 2 for each segment. Next, the softening and melting isotherms levels are calculated ( Figure 5).
Hot metal and slag chemical composition and temperature; 4.
Top gas composition and temperature.
In following tables (Tables 1-5) are shown the main measured data and modeling results representing the work of the blast furnace after BBP entering, which took place on 24 March 2018.   Table 1 shows the radial distribution of gas parameters obtained by BBP input and average gas parameters measured at the throat level for the whole furnace. Table 2 shows that the materials segmental distribution sum is balanced with the entire furnace at a satisfactory level. For all materials calculation error is below 1%, despite the fact that mass exchange which takes place between the BBP and throat levels is amended according to Equation (9). However, it would be better if the material distribution was calculated based on ABP data.
In turn, from Table 3 it can be seen that quotient α in the axis region is lower than at the wall. It means that gas flux is much higher than batch flux at the axis. This is as expected since the increased axial gas flow is a prerequisite for obtaining the cohesive zone desired shape. So, the numerical calculation of α can be accepted as satisfactory. However, the K 1 coefficient is also the result of numerical calculations, but its value is the result of closing the heat exchange balance between the throat and BBP levels. In the current model, the gas temperature radial distribution at the throat level is approximated according to Equation (13), so K 1 values may be affected by a calculation error. On the other hand, from the temperature values shown in Table 1, it can be seen that heat exchange between the BBP and throat levels could not be omitted. The ideal solution for the direct determination of K 1 is to use radius measurements on both levels with BBP and ABP. However, in lack of ABP the applied inverse method seems to be successful in the reserve zone level calculation. This is evidenced by the calculated gas temperature at the level of the burden isotherm 1000 • C; it is close to the assumed 1020 • C. Also from Figure 7, it can be seen that in the axis region the isotherm 1000 • C (red color line) is very close to BBP level and the measured gas temperature was 1024 • C. Calculations of other cases (see Supplementary Materials section) also show a good convergence of the mentioned parameters. Table 4 shows the K 2 coefficient values radial distribution, which is calculated as a heat exchange balance between the tuyeres and the isotherm 1000 • C levels. Because input data at these levels are known, the K 2 values seem to be reliable. Calculated on the K 2 basis levels of softening and melting isotherms are also in good convergence. However, in Figure 7, it can be seen that the shape of the cohesive zone (green color area) is generally the same as the shape of 1000 • C isotherm, which is assumed as the reserve zone level. It can therefore be concluded that in the present model the identification of the reserve zone level is of key importance. Figure 7 also shows the coke layer situation only in the active (permeable) part of the cohesive zone.
Summarizing the discussion, it should be noted that the current model uses BBP data as primary. However, when calculating the radial distribution of the charge and the heat transfer coefficient for the preparation zone, ABP was clearly missing and some simplifications were required. In the future, if possible, the model should be upgraded in accordance with the mentioned devices used, of which there should be at least two at different levels.

Conclusions
The paper describes the construction of a model that calculates the coke layers situation in the cohesive zone of a blast furnace; in particular, the fundamentals of mass and heat exchanges in blast furnace zones.
On the basis of measurements and chemical analyses obtained online, the model draws the situation of coke layers in the blast furnace vertical crosscut. Using only below burden probe radial data required several simplifications or assumptions in modeling. It Summarizing the discussion, it should be noted that the current model uses BBP data as primary. However, when calculating the radial distribution of the charge and the heat transfer coefficient for the preparation zone, ABP was clearly missing and some simplifications were required. In the future, if possible, the model should be upgraded in accordance with the mentioned devices used, of which there should be at least two at different levels.

Conclusions
The paper describes the construction of a model that calculates the coke layers situation in the cohesive zone of a blast furnace; in particular, the fundamentals of mass and heat exchanges in blast furnace zones.
On the basis of measurements and chemical analyses obtained online, the model draws the situation of coke layers in the blast furnace vertical crosscut. Using only below burden probe radial data required several simplifications or assumptions in modeling. It would be much more reliable to also use the above burden probe and the profilometer. However, the presented model is satisfactory to visualize the shape of the cohesive zone and to determine the nature of the gas flow through the furnace. The model was tested in real conditions and then implemented. It can be adjusted to any blast furnace equipped with BBP only.