Discrete element simulation of the charge in the hopper of a blast furnace, calibrating the parameters through an optimization algorithm

The purpose of this study is to simulate the distribution of a coarse granular material discharged in a hopper via a conveyor belt. This simulation is intended to be a model calibration for an optimization that will be later performed to obtain a proper material distribution device. From the hopper, the material is discharged in a blast furnace. Hence, achieving an adequate distribution in the hopper is crucial, since that distribution is directly linked to how the material spreads in the blast furnace, and this heavily influences the efficiency of the whole steel-making process. The apparatus is modeled by online three dimensional Computer-Aided Design software Onshape. Rocky DEM, a Computer-Aided Engineering software based on Discrete Element Method, is used to simulate the charge. The parameters of the numerical model are calibrated through an optimization algorithm. This phase is realized thanks to the optimization platform modeFRONTIER, using an algorithm that exploits meta-models to reduce the computational time of the optimization. By comparing the simulated results with the visual data obtained from blast furnace plant, the goal is to validate the model and to better understand the behavior of the whole charging process.


Introduction
The blast furnace is today one of the key steps of the steelmaking process route. The aim of the blast furnace is to transform iron ore into pig iron via a complex set of multiphase (gas, liquid, solid) chemical reactions. Briefly, the iron ore (mainly, ferrous oxide) is fed into the top part of the reactor in solid phase (lumps), it is reduced inside the reactor while slowly travelling downward, finally it is collected in the bottom part as molten pig iron. The reducing gas is the carbon monoxide produced by the combustion of coke, also charged in the blast furnace in alternate layers with iron ore, and hot air, continuously injected in the bottom part of the reactor. The blast furnace can indeed be modeled as a multiphase, countercurrent, moving granular bed, partly batch and partly continuous chemical reactor.
The efficiency of the blast furnace, normally measured as total energy consumption per ton of product, is mainly due to the right contact between gas and solid phase. It means that in every part of the reactor (with volumes up to 6000 m 3 ), the gas must have the right pressure, temperature and composition to react efficiently with the iron ore. This is somehow controlled by means of the distribution of the feed material at the top of the furnace.
Up to more than 100 ton of material can be fed in a single batch at the top of the furnace and they spread across the cross section of the furnace, with an area that can be larger than 100 m 2 . To achieve a correct distribution of the material (iron ore and coke), different technological solutions have been developed during the last century.
One of the oldest and simplest is the single bell and hopper arrangement, applied in the blast furnace this study refers to (Acciaieria Arvedi Trieste BF#3 ). In this system, the feed material is charged via a conveyor belt into an hopper and then released into the top of the furnace rolling on a conical surface. This system is viable only for small blast furnaces, since it could distribute the material only in one single circumference and not in the whole cross section of the furnace. The material must be uniformly arranged on the entire base perimeter of the hopper, which is why a rotating type, to operate during the release of the material from the conveyor belt, was adopted in this blast furnace. The mechanism, however, was at the basis of premature wear of the sealing systems, with important consequences on emissions that were no longer acceptable according to the most recent normative dispositions. For this reason, it was decided to replace the rotating hopper with a fixed type for which the necessary uniformity of charge was to be achieved by means of a deflector.
This solution has the advantages to be simple, cheaper for construction and maintenance, without rotating parts, so more reliable, but not really effective in the uniform distribution of material.
Having acknowledged that the distribution of the material in the charging hopper could promote or demote the efficiency of the overall process of the blast furnace, next step was to model the charging process in order to understand if and how to optimize it.
A Discrete Element Method (DEM) solver, Rocky DEM, has been used to simulate the charging of coke and iron ore. The DEM approach describes the motion of each particle individually, with a special treatment of eventual collisions.
DEM simulations of the charging process of blast furnace can be found in the literature, such as [23,24,30], but focused on the bell-less type charging system, while in this work the atypical charging system, described above, is considered.
In this particular industrial application, a single phase granular flow occurs, contrary to the majority of industrial processes concerning granular materials, where a secondary-fluid phase, such as air, is present.
To calibrate the parameters related to the virtual particles' motion, so that they behave just like the real material, an optimizer has been used, which is quite original. As a matter of fact, literature is full of calibration examples based on trial and error heuristic method, [4,17,18,36] just to name a few. But, there are not so many examples of such a systematic approach to material calibration. Some works on DEM model parameters' optimization are [3,13,14]. It is also worth mentioning the use of design of experiments and response surfaces for calibration, such as Taguchi methods [19], Artificial Neural Networks [6], Latin Hypercube Sampling and Kriging [25,26], Random Forest algorithm [7,8]. However, those are mainly laboratory and not industrial applications, which are easier to design and control and with available experimental measurements for verification and validation.
In this study, as mentioned above, the application is not easy to control and the direct measurements of the parameters of interest are very difficult. So that, the numerical model becomes fundamental to understand the functioning of the apparatus, giving an explanation of the phenomena detected by the pressure measurements in the blast furnace and monitored by the operators. In this sense the whole process of parameters' calibration and charge simulation done in this work is very interesting.
The paper is organized as follows: in Sect. 2 the calibration process is described in detail and the results of parameters' optimization are given. In Sect. 3 there is the sketch of the charging equipment and charging phase, with the results of the DEM simulations of coke and iron charges, obtained with the parameters' values identified by optimization. Section 4 discusses future work, whose need has been highlighted by the outcomes.

Calibration process
In DEM software the parameters regulating the particle dynamics have to be set. These parameters have physical meaning, but in some cases they can be difficult to measure. Moreover, even assuming the physical micro parameters are known for the bulk material, in order to reduce the computational time of numerical simulation, the virtual material is frequently different from the real material, in terms of physical micro parameters, particle shapes and particle size distribution.
So, in general, a reverse calibration of numerical parameters is a key step to achieve a similar behavior between the virtual and real material. The numerical parameters have to ensure the macro properties of material, which are experimentally easy to measure. For granular material the angle of repose, which is the steepest angle of descend to which a material can be piled without slumping, is essential to describe the macro-behavior [2], so it is widely used for calibration purposes. The idea is to find the virtual material parameters which ensure a repose angle as close as possible to the value typical of the real material.
In order to simulate the process of hopper burden charging, the parameters of two materials have to be calibrated: coke and iron. The latter one is made up of 65% of total weight by sinter and 35% by pellet. Typical repose angles of the materials loaded in the hopper are given in Table 1.
Below there is the description of the workflow used for parameter calibration, in particular the approach used to calculate the repose angle, and the results of calibration process for materials loaded in the hopper. Figure 1 shows the workflow used for calibration of DEM solver parameters. The solver is Rocky DEM, which is a commercial software based on Discrete Element Method. Extensive description of the original DEM method can be found in [31,38]. For this study, the following models were selected: Hysteretic Linear Spring model for normal forces; Linear Spring Coulomb Limit model for tangential forces; Type 3 model for rolling resistance. In [15] detailed information related to these models is given.

Workflow of calibration
The calibration of DEM solver parameters is performed through an optimization algorithm, more efficient than a trial and error approach, allowing to explore the design of experiments with the focus on the search of the best parameters, which ensure a similar behavior between the virtual and real material. In the optimization problem under study, there are 4 variables and 1 objective.
The variables are a subset of the numerical parameters to be set in DEM solver. In [33] it is demonstrated that the friction coefficients mainly influence the behavior of particles and in [26] it was found there is no correlation between static and rolling friction coefficients, although this has sometimes been assumed in published research. Given these premises and after some numerical tests, it has been decided to optimize only a subset of the numerical parameters. The parameters to optimize are: The assumptions on the other parameters are described and explained in Sect. 2.1.2.
The objective of the optimization problem is to find the numerical parameters which minimize the difference between the repose angle of the virtual material ( repose_angle ) and the real material ( expected_repose_angle). Table 1 Angles of repose by [16] Material Repose angle

Optimizer
The algorithm used to optimize is pilOPT, which is a proprietary algorithm, designed to be multi-strategy and self-adapting, smartly combining local and global search. This particular application is single objective, so it does not fully exploit the potential of the method, capable to find the Pareto front for multi objective problems, but it has been chosen because it balances real and Response Surface Method (RSM)-based optimization. In this way it reduces the number of real experiments necessary to reach the optimum.

Assumptions
The real materials loaded in the hopper are not-cohesive and dry, as well as the virtual material used in DEM simulations. However, some assumptions are made about the particle geometries, size distributions and micro properties of the virtual materials.
In different papers, such as [5,10,32,34], the influence of particle geometries have been demonstrated on material macro properties and, in particular, [9,20] investigate the influence of particle shape and shape approximation on hopper discharge. But the particles of the materials loaded in the hopper have very irregular and different shapes, not even easy to identify and estimate. Even assuming a distribution of various geometries among the particles, the computational cost of the simulation would increase so much that it would be too expensive in an optimization process. So particles of materials under study are assumed spherical. The effects of real geometries will be obtained as much as possible thanks to the calibration of the virtual material.
The particle size distribution for real materials are known. They are provided in Fig. 2a by weight percentage.
Since the dimension of particles is one of the parameters that determine the time step used in the simulation, some assumptions have been done on this aspect in order to reduce the computational time. In particular, a scale-up of particle size has been introduced.
In literature, different approaches of particle scale-up have been developed and investigated [12,27,28,32]: • exact scaling, where particle number is constant while upscaling of all lengths; • coarse-graining, which consists in reduction of the number of particles while container size is constant; • scalping or cut off, which means the coarse-graining approach is applied just to the finer fractions of the real bulk material, obtaining a reduction of the number of particles while container size is constant.
In this application, both coarse-graining and cut off have been used. Analyzing the effect of particle size distribution of the coke, it turns out that considering particles with a diameter in the interval 10 ÷ 100 mm, instead of 5 ÷ 100 mm, it is possible to double the time step, reducing the computational time of the whole simulation. Moreover, the percentage of particles with a diameter in the interval 5 ÷ 10 mm is 1% , that is very low. So the cut off of coke particle size distribution has been applied for diameters lower than 10 mm and percentage of the particles with diameter within 5 ÷ 10 mm has been joined to that relative to the interval 10 ÷ 20 mm, obtaining the distribution in Fig. 2b, which has been used for the numerical simulations.
For iron, which is is made up 65% by sinter and 35% by pellet, both coarse-graining and cut off techniques have been used. Since the 88% of sinter particles have diameter lower than 40 mm, the remaining 12%, with diameters between 40 and 80 mm, has been included in the interval 30 ÷ 40 mm, so to reduce the difference between the minimum and the maximum diameter. However, the particle size for sinter and pellet is very low. Median value of particle diameter for sinter is about 14 mm and for pellet 8 mm, versus 50 mm of coke particles. In the simulation of iron hopper charge, the number of sinter and pellet particles would have been so high to make impossible to simulate the process in a reasonable computational time with the available hardware. For this reason, the operation of coarse-graining was performed and the real dimensions of particles were doubled. The assumed particle size distribution for sinter is given in Fig. 2b.
For pellet particle size distribution, 100% by weight of particles is considered to have diameter within the interval 10 ÷ 22 mm.
Not all parameters, to be set in DEM solver, have been calibrated. Below there are the assumed values for these parameters: • bulk density: equal to the bulk density of real material, which is 508 kg∕m 3 for coke, 1863 kg∕m 3 for sinter and 2088 kg∕m 3 for pellet; • Young modulus: a number of techniques exist for minimizing the computational cost of DEM simulations and one such method is the reduction of particle stiffness, which allows for bigger time steps and therefore fewer iterations in the simulation [21,22,35]. Being this parameter inversely related to the time step, it is set as low as possible: 10 7 Pa , as in [37]. The value set is much lower than the real one, both for coke, sinter, pellet and wall of the hopper. However, it is common practice to keep it low as long as the behavior is as expected, and here this is provided by the parameters' calibration; • Poisson ratio: set at 0.25 [1]; • particle-wall restitution coefficient: to simulate the different bouncing conditions that can occur when particles collide with the walls, the material is calibrated using a mix of two materials with particle-wall restitution coefficient set at 0.2 and 0.5 respectively.
Since only spherical elements are used, it is not possible to simulate in any other way the different possible bouncing conditions caused by colliding either an edge or a flat area of the particles; • damping coefficient for long-term contacts: set at 0.5 after conducting tests that highlighted the low influence on the resulting repose angle, while it was found useful in terms of avoiding infinitesimal fluctuations caused by long term contacts; • dynamic friction coefficients: both particle-particle and particle-wall dynamic friction values are set at 90% of the static value, based on the fact that typically these values are slightly higher than the dynamic values [7,15].

Angle of repose test
In [11] there is an interesting review of the calibration methods of DEM parameters proposed over the last 25 years. The angle of repose test has become a standard for the calibration and in literature a lot of different methods to measure it are discussed [10,18]. But, the most used method remains the conical pile test, both for non-cohesive and cohesive materials [29], which consist in the lift of a cylinder, previously filled with the bulk material, so that the material is free to dispose itself according to the angle of repose.
The phases to simulate to replicate the experiment are: In [27,28] it has been demonstrated that lifting velocity of cylinder must not exceed a limit value (10 mm/s) in order to not influence the final result.
In fact, some tests on our specific problem using this approach (Fig. 3), have confirmed that, even if the speed is kept low (around 30 mm/s), the result is influenced by the number of particles involved in the simulation and the repeatability of the experiment is not assured.
Since a lower velocity would increase the calculation time too much, in order to avoid the influence of the lifting Fig. 3 Different screenshots of the cylinder removal simulation in different simulated moments. Particles coloured per dimension velocity, it has been decided to make the cylinder disappear at once, which means, in practice, disable its geometry in the software. In this manner, after the cylinder is filled with bulk material and after the particles completely settled under gravity, the cylinder is instantly removed. This operation allows the particles to fall freely (Fig. 4). This choice also grants a computational time reduction because the rising cylinder phase is avoided. Another important decision, in order to favor this aspect, is that the base where the material settles, has limited diameter letting the exceeding particles to fall and exit the computational domain. Thanks to that, the solver does not consider them anymore and the computational cost reduces because of the gradually decreasing of the total number of particles for which to solve the motion.
Comparing the two approaches, that is the slow removal of cylinder and its instantaneous disappearance, the time interval to be simulated is about 31.5 s in the first case and 7.5 s in the second one. In fact, filling the cylinder requires 0.5 s and 1 s is the time needed for particles settlement under gravity for both approaches. The removal of cylinder with velocity equal to 30 mm/s, supposing to lift it up to 1000 mm to free all particles, takes 30 s, during which the angle of repose is reached. Whereas, in case of instantaneous disappearance, only 6 s are required for particles to settle and reach the repose angle. Let us underline we are referring to simulated time interval and not the duration of the simulation.
To ensure the repeatability of the experiment with the cylinder disappearance, various dimensions of cylinder and base were tested, resulting in different values for coke, pellet, and sinter. The chosen values for the diameter of the base and the diameter and the height of the cylinder are given in Table 2. The total time simulated in all test is equal to 8 s. Table 3 shows the computational time required by each simulation and the number of particles with which the cylinder is filled and the final number of particles which form the angle of repose with the underlying considered base.
Once the simulation of the conical pile experiment is over and the material is settled according to the angle of repose, the angle is measured through the following steps:   1. the pile of particles is divided into N circular sectors; 2. every circular sector is itself divided in the radial direction and, in each radial sector, the particle with the highest height is detected; 3. a regression line is determined using the heights and the radial positions of these particles in each circular sector so to measure the repose angle; 4. an average value of the repose angle among circular sectors is calculated from the results of the previous step.
In Figs. 5 and 6 an example of circular sector and an example of the distributions of the repose angle values in each circular sector are shown, respectively.
The number of sectors, in which the pile is divided, has been chosen after various tests performed in order to have enough particles in each sector. The number of sectors is different for the calibrated materials since the particle dimensions are different. In table 4 there is the used number of radial and circular sectors.

Results
In Table 5 there are the results obtained by the calibration of parameters, using the optimization workflow previously described. In the table the optimum sets of parameters are reported for coke, sinter and pellet. The objective of the optimization is to minimize the difference between the repose angle resulting from the numerical experiment and the expected one, that is the mean value of the range reported in Table 1. For coke calibration the target value of angle of repose is 36.5 • , for sinter 31 • and for pellet 25 • . The single objective optimization converges when the objective can not be further improved.
The final settlements of calibrated materials obtained by the conical pile test with the optimum parameters are shown in Fig. 7.
The optimum parameters have been reached by 26 simulations of the conical pile test for coke, 55 for sinter and 16 for pellet. In Table 6 there are the number of Design Of Experiments (DOE) simulations which have given feasible, unfeasible and error result, where feasible means the repose angle is within the experimental range, unfeasible means the repose angle is out of the experimental range and error means the simulation has crashed without result.    An analysis of the designs created during the optimization phase has been carried out, in order to understand the effects of the four parameters on the repose angle. The only information that could be gathered from the available data was the strong direct correlation between rolling resistance and repose angle, as expected, since this parameter also simulates the shape of the particles.
To verify the repeatability of the conical pile test, two further simulations are carried out with the optimum set of parameters. The resulting angles of repose are reported in Table 7 for each calibrated material.
Repeating the test, the obtained values belong to the ranges reported in Table 1, except for the pellet, for which one of the values are out of range. But, being the value very close to the expected range, it has been considered acceptable.
As already mentioned, the materials loaded in the hopper are coke and iron, which is a mix of sinter and pellet. Since no repose angle of the mix is available, but only the repose angles of the materials within the mix, to handle interactions between sinter and pellet particles, the parameters of the mix are obtained by averaging the values of the two materials, being aware this is a risky choice but it will be eventually justified by the results. The parameters for the mix, which constitutes the iron, are reported in Table 5.

Simulation of the charge
The simulation of the charge of the hopper has the purpose to better understand how the burden is currently distributed inside, both for coke and iron charges. The numerical results are compared with the information deduced by the pressure distribution measured in the furnace, the erosion found inside it (more evident where the piles are higher) and some photos, available only for coke charge. The voids of coke particles are highlighted in Fig. 12a, where are also pictured the position of the rod and the deflector, positioned just under the hatch where the conveyor belt releases the materials. Because of the presence of a two-skew deflector, two piles of material form. Based on what detected inside the furnace, piles currently form as pictured in Fig. 8, with coke having one of the piles much higher than the other, whereas for the iron they have similar heights.
The red circles in Fig. 8 are representative of peaks of the piles, while the black ones are the zones where voids were found, generated by the presence of the deflector and the rod that moves the bells, which creates a shadow cone. Note that in Fig. 8 the deflector is not drawn, because the plant design is prior to the addition of the device, which is placed inside the hopper exactly below the hatch through which the material passes.

Description of the apparatus
The mutual placement of the hopper and conveyor belt, as well as the velocity of the belt, are very important to correctly simulate the charge of material in the hopper. To position the belt the reference is the point given by the intersection between the symmetric axis of the rod at the center of the hopper and the section identified by the contact surface between the hopper and the lid. This point is highlighted in Fig. 9. The distance between this point and the mid-point of the axis of the drive drum define completely the position of the conveyor belt.
Next some features of the conveyor belt are summarized: • distance between the reference point and the center of the drive drum: 2050 mm (height distance), 1541 mm (distance orthogonal to drum drive axis), 927.5 mm (distance parallel to drum drive axis); • drive drum diameter: 750 mm; • tilt angle: 4 • ; • number of rollers: 3; • thickness: 12 mm; • speed: 1.79 m/s. Figure 10 shows the geometry of the entire system used for DEM simulation of the charge, designed with Onshape Computer-Aided Design (CAD) software and used in Rocky DEM solver.
The charges of two different materials are simulated. The quantities of material to be loaded are different  Table 8 shows the values for the coke and the iron charges.
The simulation has to take into account not only of the charge duration, but the time interval needed by the material dropped on the conveyor belt to reach the point where it leaves the surface of the belt and an additional time interval for the settlement of the material under gravity at the bottom of the hopper. So the simulations of the charges are extended by 10 s. The total simulated time for coke is 82 s and for iron 150 s.
Let us underline that coupled DEM-Computational Fluid Dynamics (CFD) analysis is not needed in this application because the interaction between particles and gases can be neglected, being the particle drag at least one order of magnitude lower than inertial forces. So, a single phase granular flow is investigated.

Results of the charge simulation
The charges of coke and iron within the apparatus described above have been simulated using the numerical parameters obtained by the calibration process.
The computational time needed to simulate the coke charge has been 5 h and half, where the number of particles involved is 280,000. For the iron charge the particles are 1,400,000 and the computational time has been 25 h.
These simulations were possible within a reasonable timeframe for industrial applications thanks to the exploitation of one Nvidia RTX-6000 GPU, able to reduce six times the computational time necessary for solving the same problems using an Intel(R) Xeon(R) E5-2630 v3 processor (8 physical cores).
The results for the coke charge inside the hopper are shown in Fig. 11. In Fig. 12 there is the comparison between the distributions of coke viewable in the available photos and obtained by DEM simulation, which confirms the reliability of the results.   Figure 13 shows the results for the iron charge inside the hopper.
For coke, the separation of the charge into two different piles is evident. The one which is the farthest from the conveyor belt is the highest. Two areas are visibly depressed with a much smaller amount of material. One is under the deflector and one is beyond the bell rod.
For iron, the charge separates into two different piles, too. But, unlike the result found for coke, the two piles have similar heights. The depression areas, on the other hand, are arranged in the same way as the coke.

Comparison between the distributions of coke and iron inside the hopper
The two charges, having different characteristics (flow rate, density, particle size, etc.) have different final distributions, mainly referring to the height of the piles. The reason is the difference of material distributions on the conveyor belt, as shown in Fig. 14. The amount of material and its maximum height is greater for coke. This implies that the coke particles furthest from the belt surface have a greater throw, impacting on the further face of deflector and contributing more to the creation of the pile further from the conveyor belt.
The comparison between the two simulations also showed selective behavior of the belt: the pile furthest from the conveyor belt has spheres with a larger diameter, while the nearest pile has smaller spheres, highlighted in Fig. 15.
This behavior is because the spheres with a smaller diameter form a layer on the surface of the belt, while the larger ones lay above the layer of small particles (see Fig. 14). The consequence of this distribution is that the large particles, being higher, have a greater throw, while the small ones have a smaller throw. So the large particles impact the further face of deflector and create the pile further from the conveyor belt, the small particles impact the closer face of deflector and create the pile closer to the conveyor belt.

Conclusions
The present paper investigated the charge distribution of the materials loaded in the hopper of an existing blast furnace. The purpose was to evaluate the effects of a deflector which has been inserted inside to make the distribution of both charges, coke and iron, more uniform.
Some assumptions have been made to simplify the process to be simulated, in order to both reduce the computational time and keep the simulation plausible. So the numerical parameters of both materials were calibrated. The numerical parameters had to ensure the angle of repose proper of material. An automated calibration framework based on an optimization algorithm, which exploits alternately the DEM numerical evaluations and a meta-model of the objective function, was used. This feature of the algorithm allowed to reduce the optimization run time.  11 Different views of the particle distribution obtained by the simulation of coke charge. The color of the particle represents its height The values of parameters, given by the calibration process, were used for the charge simulation of coke and iron in the hopper. So, the distributions of the materials were obtained.
The results highlighted the inability of the deflector to uniform the distribution, sufficiently to ensure that the blast furnace operates properly.
In particular, in addition to confirming the formation of two piles, it pointed out that the large and small particles are clearly separated during the charging phase. So that one pile consists mainly of small particles and the other of large particles. The shape of the deflector emphasizes this behavior. This last aspect was not so evident from the information available before the results of the simulations.
This study highlighted the need to replace the device with a more performant one and it has allowed to lay the foundations to proceed with the optimization of the current deflector geometry.
Moreover, this work allowed to validate a method for automatically calibrated materials requiring a reduced amount of time thanks to the optimized, in terms of the time, repose angle test, granting also the repeatibility of the output value. In fact, time is a very important variable in this kind of problems, since the simulations are quite heavy and it is difficult to obtain the results in the required time with, in many cases, limited computational resources.
Concluding, the main contribution of the paper is the methodological approach to a general problem, when a required material distribution must be obtained in a container, silo, hopper or similar, from complex material release conditions. This is the case of the system under study, that is the one bell-type charging system of an operational blast furnace at the time this research was Fig. 14 Comparison of the material distribution on the conveyor belt between coke (on left) and iron (on right). The color of the particle represents its height Fig. 15 Comparison of the distribution of particles by size between coke (on left) and iron (on right). The color of the particle represents its size conducted. In this application the achievement of a correct distribution of the material is fundamental, which is a requirement in many other contexts. So, the presented methodology can be extended to other applications and backgrounds, such as food, chemical, pharmaceutics industry.