Detailed modelling of packed-bed gas clogging due to thermal-softening of iron ore by Eulerian–Lagrangian approach

Chemical Engineering Journal xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Chemical Engineering Journal journal homepage: www.elsevier.com/locate/cej

Detailed modelling of packed-bed gas clogging due to thermal-softening of iron ore by Eulerian–Lagrangian approach Shungo Natsuia,b, , Shingo Ishiharaa, Tatsuya Konc,d, Ko-ichiro Ohnoe, Hiroshi Nogamia ⁎

a

Institute of Multidisciplinary Research for Advanced Materials, Tohoku University, Katahira 2-1-1, Aoba-ku, Sendai 980-8577, Japan Division of Materials Science and Engineering, Faculty of Engineering, Hokkaido University, Kita 13 Nishi 8, Kita-ku, Sapporo 060-8628, Japan c Magnetic Powder Metallurgy Research Center, National Institute of Advanced Industrial Science and Technology, 2266-98 Anagahora, Shimo-Shidami, Moriyama-ku, Nagoya, Aichi 463-8560, Japan d Research Center for Steel, Kyushu University, Motooka 744, Nishi-ku, Fukuoka 819-0395, Japan e Department of Materials Science and Engineering, Faculty of Engineering, Kyushu University, Motooka 744, Nishi-ku, Fukuoka 819-0395, Japan b

HIGHLIGHTS

GRAPHICAL ABSTRACT

3D dynamic model for a packed • New structure with coke and ore is proposed.

can predict the packed bed perme• Itability, even with softening-ore. permeability characteristics of • Gas ore-layer shapes in ironmaking process discussed.

local pressure drop was estimated • The from the geometrical data of softening packed-bed structure.

formation of ineffective voids in • The the gas-flow path increase the local pressure drop.

ARTICLE INFO

ABSTRACT

Keywords: Packed bed Iron-ore softening Gas permeability Discrete element method (DEM) Eulerian–Lagrangian coupling Computational fluid dynamics (CFD)

Eulerian–Lagrangian numerical scheme is applied for analysing packed-bed-structure constructions involving non-spherical solids, such as metallurgical cokes and ferrous ores, and the high-temperature softening characteristics of such beds. 3D scanning is applied for determining the coke and ore shapes, and a multi-sphere discrete element method is used as the functional scheme for non-spherical solid-particle motion tracking. The transient deformation behaviour of the softening ore is simulated using the advanced discrete element method, and the gas permeability characteristics exhibited by the ore shapes in the ironmaking process are discussed. Based on this model, cases with varied softening behaviour represented by the joint spring coefficient are investigated and the effect of the ore-softening behaviour on the gas permeability is evaluated. It is established that the pathway of the passing rivulet depends upon the softening-ore deformation behaviour.

1. Introduction In response to the recent global warming issue, steel industries across the world have been required to decrease carbon dioxide

⁎

emissions. The ironmaking process, which converts the iron oxide in the ore to metallic liquid iron, consumes approximately 70% of the total energy input to integrated steelworks [1]. The blast furnace, which is a large reactor with a height of 30 m and volume of 2000–6000 m3, and

Corresponding author at: Institute of Multidisciplinary Research for Advanced Materials, Tohoku University, Katahira 2-1-1, Aoba-ku, Sendai 980-8577, Japan. E-mail address: [email protected] (S. Natsui).

https://doi.org/10.1016/j.cej.2019.123643 Received 5 May 2019; Received in revised form 20 September 2019; Accepted 28 November 2019 1385-8947/ © 2019 Elsevier B.V. All rights reserved.

Please cite this article as: Shungo Natsui, et al., Chemical Engineering Journal, https://doi.org/10.1016/j.cej.2019.123643

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Nomenclature

v x

Symbols

A c d d D F g h i I k i ', j', k' l L m N nx ny nz P Q q q R r Re Sf

t T u v

velocity, m s 1 positional vector of an element in a rigid body, m

Greek letters

projected area, m2 cell type function, diameter, m relative displacement increment, m tube diameter, m force, kg m s 2 gravity, m s 2 height of the bed, m particle index, moment of inertia, kg m2 spring coefficient, kg s 2 cell indices in the x, y, and z direction, respectively tube length, m representative length, m mass, kg number of particles in the effective radius, dimensionless number of particles in a solid, maximum cell number in the x direction, maximum cell number in the y direction, maximum cell number in the z direction, pressure, kg m 1 s 2 quaternion normalised interparticle length by the influence radius, components of the quaternion rotation matrix position vector, m Reynolds number, dimensionless deformation degree due to softening under gravity, time, s torque, kg m2 s2 gas velocity, m s 1 void-type function,

µ

void fraction, dashpot coefficient, kg s 1 rotational angle, rad Darcy friction factor, viscosity, kg m 1 s 1 density, kg m 3 sphericity, dimensionless angler velocity, rad s 1

Subscripts 0 a c eff f h i init j g max min n p proj s v void wet

contains layered iron-ore and coke (carbon source) packed beds, is the dominant method for the production of primary iron, despite the increase in direct-reduced-iron production. The cohesive zone, in which the ore is softened and melted, is focussed upon; it drastically affects the gas-flow pattern in the furnace and has considerable effect on not only the productivity of the process but also the carbon dioxide mitigation efficiency [2]. The phase-changing ores can be deformed due to powder pressure and exclude the gaps between the ferrous burden pieces. In the cohesive zone, the softening and melting oxide phase (molten slag) and the metallic-iron phase coexist. The softening and melting of slag are affected by the existing oxide quantity, along with their distribution, morphology, [3], and degree of reduction [4]. In general, the first

inlet joint spring contacted effective frictional hydraulic particle index initial neighbouring particle index aroundi centre of gravity maximum minimum normal particle projection shear equivalent volume void wetted perimeter

oxide-melt is formed at the interface of the lowest melting point between FeO and another oxide. If the ore reaches the softening temperature with a higher degree of reduction, a porous solid Fe shell is formed, which confines the solid and liquid oxides, after which it becomes a semisolid material. Thus, after the liquid slag in the core has wet all the oxide particles, the mechanical strength of the core is considerably reduced, and the resistance to deformation is determined by the metallic Fe shell [5]. These phenomena are depicted in Fig. 1. With the increase in the liquid volume fraction, a yield point is reached, where the metallic Fe shell can no longer hold the liquid. Thus, ore softening is considered to be the moment, when the metallic burden can no longer resist the action of mechanical forces [6,7]. Due to the

Fig. 1. Schematic illustration of the softening mechanism of reduced iron ore. 2

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

complexity of the packed bed system, the flow in such coke beds has generally been modelled on a macroscopic scale by a semi-empirical model assuming homogenous packing [8]. The complexity of the flow field increases significantly, when the softening of the solid is coupled with the hydrodynamics. Although the average porosity determines the distribution of the passing gas flow, the void structure and connectivity, which are strongly modified by softening and melting, determine the local flow and pressure loss. This study aims to track the changes in the packed-bed structure due to the initial ore softening deformation, based on computational dynamics. Recently, numerical analyses using the Lagrangian approach, which can track moving calculation points, have been applied for multiphase flow problems in a high-temperature metallurgical process with a complicated interface [5,9-17]. Furthermore, the multi-sphere (MS) extended DEM has been adopted as an applicable scheme for the simulation of non-spherical solid-particle motion [18-23]. This method can simulate the behaviour of the solid-particle aggregate with nonspherical complicated shapes; the gas flow model for the packed bed, namely, the grid-based Eulerian continuous fluid model, is often used for analysing heat and mass transfer. Recently, the Euler–Lagrange (EL) coupling scheme, e.g., the CFD-DEM scheme, which enables the investigation of the interactions between the gaseous and solid phases, was proposed [26-36]. The advantage of this method is the explicit tracking of the time-varying interfaces of each phase, which is difficult to observe experimentally. The E-L coupling scheme has been extended to analyse the flow field at sub-pore length-scales in the packed bed for mineral processing [37], petroleum engineering [38], catalytic reactors [39] etc.; meanwhile, there are few studies in which phenomena causing large non-steady elastic deformation such as ore softening are applied to packed beds [25]. The advanced discrete element method (ADEM), which discretises a continuous viscoelastic phase by moving particles [24], is suitable for the analysis of the interfacial deformation of the ore phase. The schematic of the E-L coupling simulation scheme for a packed bed including softened deformation ores is shown in Fig. 2. Based on this approach, the deformation of the softened ores in the coke bed, under gas flow conditions, in accordance with the strength and properties of the individual ore particles can be directly determined. In this study, a new 3-dimensional dynamic model for a packed structure with two different types of particles (coke and ore) is proposed. This model can treat deforming ore particles and hard coke ones. With this feature, this model can reproduce the above-mentioned complicated packed-bed structure. Using such a detailed structure, the permeability of the bed can be accurately predicted. In particular, a pressure-drop

Fig. 3. Geometrical representation of a) quaternion rotation operation and b) relative position of the connected particles in a joint spring approximation.

disordering process is expected. Therefore, the purpose of this study is to investigate the influence of the ore thermal softening behaviour and 3-D gravitational deformation on the gas-flow characteristics in detail. 2. Mathematical modelling 2.1. Non-spherical solid motion calculation procedure The DEM [40] has already been expanded to handle the motion of

Fig. 2. Schematic illustration of the Eulerian–Lagrangian coupled model including the ore-particle deformation behaviour. 3

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

mg

Ig

Nc

d vg

=

dt

d

Fc, i + Fg .

(1)

i=1

Nc

g

=

dt

Tc , i.

(2)

i=1

Solid motion depends upon Newton’s second law, and Fc was calculated with the “soft sphere” assumption, based on the Hertzian contact theory [40]. DEM is popular scheme for predicting the movement of spherical particles by solving motion equations that consider the contact forces among the particles. These contact forces are described by a Voigt model that consists of a spring, dashpot, and slider. Fc was calculated using the following equations: (3)

Fc = Fc, n + Fc , s Fc , n =

kn dn

ks ds

Fc , s =

dn . t

n

s

µf Fc, n;

ds ; t

(4)

|Fc, s | < µf |Fc, n | otherwise

.

(5)

where k is the spring coefficient, is the dashpot coefficient, d is the relative displacement increment at discrete time, t , and µ f is the frictional coefficient. Further information on the calculation of k and can be found in a previous paper [42]. Next, we calculate the translational velocity vector of aggregate, vg , the angular velocity, g , and position vector, x g , with respect to the centre of gravity of each “solid”, as follows:

vg =

1 n

n

vi .

(6)

i=1 n

g

= Ig

1

{(xi i=1

xg =

Fig. 4. 3D-scanning STL data of the ore and coke shapes with diameters ranging from 21 to 27 mm.

1 n

x g) × vi}.

(7)

n

x i. i=1

(8)

The rotational angle, , was considered to describe the rotational motion as transformed by the quaternion, Qi . The quaternion is commonly used to minimize the numerical error of the rotational angle, and is composed of four components including a scalar and vector part. With quaternion representation, when the relative position vector from the centre of the mass of one solid to each sub-particle representing the complex-shaped solid is set to qi = xi x g , at time, t , quaternion, Qi t = (qinit , qi) = (qinit , qx , qy , qz ), is related to the conjugate quaternion, Q , as follows, considering the time steps:

arbitrary-shaped solids. It generally arrays small spherical particles and expresses complex-shaped solids to realise low calculation costs and intuitive mounting [41]. Here, we refer to the elements or calculation points constituting an aggregate as “particles.” The particles i within a solid conform to the equations governing the respective translational and rotational motions of the aggregate. When the interparticle contact force, Fc, and gravity are considered to be the forces acting on each particle, the governing equations are as follows:

Qi t = QQi t

tQ

(9) Fig. 5. Comparison of the ore-surface shapes derived from a photograph, STL data, and particle assembly, respectively. To obtain c), a level set function was employed, related to the distance from the central position of each particle arranged in a square lattice at the interval-particle diameter to the surface of the triangular polygon embedded in the volume.

4

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 6. Changes in the calculated packed behaviour of 100 numbers of coke (grey colour) and 50 ores (orange colour) of 15 different shapes, i.e., 10 pieces of 10 representative coke samples and 10 pieces of 5 ore samples were numerically generated. Time marching of one of the packed behaviours is displayed. When the velocities of all the particles were less than 0.01 m/s, the calculations were considered static and all the solids were fixed in space as the initial condition for the next simulation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

rotational axis. The transformation to the rotation matrix, R, using quaternion components is expressed as follows [39]:

Table 1 Shape properties of each sample. Solid #

O1 O2 O3 O4 O5 C1 C2 C3 C4 C5 C6 C7 C8 C9 C10

Calculation points per particle, n/-

Equivalent sphere diameter, dv/m

Average projected area, A/m2

Sphericity, ϕ/-

6703 6388 6039 6495 6547 9546 8587 10,062 8837 9239 9978 5219 8981 4922 5113

0.02339 0.02302 0.02259 0.02315 0.02321 0.02632 0.02541 0.02679 0.02565 0.02603 0.02671 0.02152 0.02579 0.02110 0.02137

0.0005440 0.0005416 0.0005029 0.0005646 0.0005302 0.0006758 0.0005998 0.0008082 0.0006252 0.0006383 0.0006827 0.0004049 0.0006900 0.0004066 0.0003922

0.7901 0.7685 0.7973 0.7455 0.7980 0.8051 0.8452 0.6971 0.8266 0.8339 0.8207 0.8984 0.7570 0.8604 0.9148

1

R = 2qx qy 2qx qz

g

| g|

sin 2 , Qi t is equal to Qi t

t

rotated by

2qz 2 2qx qy 2qinit qz 1

2qx 2

2qinit qy 2qy qz

2qinit qz 2qx qz 2qz 2 2qy qz 2qinit qx 1

2q x 2

2qinit qy 2qinit qx . 2qy 2

(10)

As shown above, at time, t , quaternion, Qi t , is expressed by its position relative to the centre of mass/gravity:

Qi t = (0, Rq it

(11)

t).

Therefore, the transformation from the quaternion to the particle position vector, x i t , comprising the solid can be expressed as follows:

xi t = x g + qi t = x g + Rqi t

t.

(12)

As the non-spherical solid takes an arbitrary shape and position, the inverse matrix I 1 of the solid’s inertia tensor at each step must be determined in order to solve Eq. (7). Here, for simplification, at time, t , the inertia tensor, Ig t , is obtained using the initial value, Ig , init , of the solid inertia tensor and the transposition matrix, R , of the rotation matrix.

Fig. 3a shows the schematic of the rotation operation with quaternion representation. Because the quaternion around the rotational = | g | t at axis, g | g |, is expressed as the rotational angle

qinit = cos 2 , qi =

2qy 2

It = RIg , init R .

around the

5

(13)

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 7. Calculation results of the dynamic equilibrium state of each softened-ore condition. Softened ores and cokes bed that became static at each joint spring constant, are displayed. The orange, and dark and light grey colours indicate the ores, coke, and gas phase, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

As Ig , init and Ig , init 1 are symmetrical tensors, the inverse matrix, of the inertia tensor can be expressed as

Ig

1t ,

Ig

1t

2.2. Ore-softening model In this study, to simulate the dynamic behaviour of the softened ore particles, ADEM simulation [46] was used, in which neighbouring particles are connected to one another by a joint spring, as shown in Fig. 2a. Theses joint springs represent the normal and shear interaction forces between the connected particles. The viscosity is represented by a dashpot inserted in parallel with the joint spring. The interaction force by the joint spring is similar to the contact force in ordinary DEM; however, there is a difference in definition of the relative displacement. The interaction force due to the joint spring, Fa , was calculated using the following equations:

(14)

= RIg , init 1R

Therefore, there is a need to find the value of Iinit determined by the initial stance of the solid. Here, it is practically impossible to determine Ig , init for an irregular-shaped solid by the analytical method. The method of integrating the inertia of the particles constituting an irregular-shaped rigid body is generally used to obtain the strength of the rotation tensor [44,45]. For the mass-point inertia tensor in relative vector, qi = (qx , qy , qz ), from the centre of the mass, we assume a cube with side length, dp, and find the initial inertia tensor, Ig , init , of a single solid as follows [43]:

Ig , init =

n i

mi d p 2 1 0 0 0 1 0 + mi 6 0 0 1

qy 2 + qz 2 qx qy qx qz

qx qy qx 2 + qz 2 qy qz

qx qz qy qz

(16)

Fa = Fa, n + Fa, s.

.

qx 2 + qy 2

Fa, n =

ka, n d a, n

Fa, s =

k a, s d a, s

a, n va, n . a, s va, s .

(17) (18)

Here, ka is the joint spring coefficient and a is the dashpot coefficient of the joint spring. Fig. 3b shows the relative displacement of the connected particles. In order to express the partial deformation behaviour of cluster particles, relative position changes between the connected primary particles are allowed. The motion of the entire cluster particle is caused by the propagation of the relative position change between the local primary particles. The relative displacement of the joint spring, da , was calculated as follows:

(15) If the resolution is insufficient, i.e., without sufficient small-particle diameter, the numerical error of the rotational movement of the solid increases. For the precision calculation of the solid-shape representation, dp should be as small as possible. In this study, dp = 1.00 mm was adopted as a constant value considering the computational cost also. 6

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

da = rij '

d a, n = 1

d a, s = d a

(19)

rij.

|rij | |rij '|

d a, n

r 'ij.

(20) (21)

Here, rij is the initial relative position vector of the joint spring and rij ' is the relative position vector after the movement of the particle. The connecting length of the joint spring is set as 1.5 times the particle diameter; i.e., the particles within this length from the particle of interest in the initial condition are connection by i springs. This value was determined considering the veracity and stability of the calculation. As the joint spring constant is correlated to the Young’s modulus, it is possible to represent various physical states, such as soft and hard, depending upon the magnitude of its value. In ADEM, the softening and deformation of a cluster particle is treated as the interaction force between the primary particles by the joint spring; on the other hand, collision detection and collision response in the contacts between different cluster particles are calculated based on the primary particles. The contact force between cluster particles is calculated based on the Hertzian contact theory, similar to conventional DEM. In Hertzian elastic contact theory, it is assumed that the size of the contact portion is “much smaller” compared to the size of the object. If the Young's modulus of the softening solid was given as lower value than before deformation in order to evaluate the contact force of the softened body, the amount of overlap between particles becomes large, which violates this assumption. Therefore, in order to avoid this paradox, the Young's

Fig. 8. Dependence of the ore-layer height characteristics on the equilibrium state joint spring constant.

Fig. 9. Calculated phase map of the horizontal cross-section for each height in the dynamic equilibrium state, for each softened-ore condition. The orange, and dark and light grey colours indicate the ores, coke, and the gas phase occupied region, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

7

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 10. Relationship between the void fraction and joint spring coefficient in each arbitrarily packed bed structure. The 3D views on the right-side correspond to the 21 i’ 80 and 21 j’ 80 regions with various height levels of each softened-ore condition. The orange and dark grey colours indicate the ore and coke, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

modulus was used as constant value, which is known before softening, for calculation of the contact force, and we evaluated deformation only using the joint springs. In this method, Hertzian elastic force correctly evaluates the load under gravity, and the deformation behavior is solved by relying only on the joint spring coefficient.

discretised using the 2nd order Adams–Bashforth method, 3rd order upwind difference (UTOPIA scheme) [48], 1st order forward difference, and centre difference scheme with 2nd order accuracy, respectively. For determining the pressure of each calculation point, the Poisson equation of the pressure correction at each time step derived from the continuity equation requires iterative calculation and was solved by applying the Bi-CGStab method [49]. It was assumed that the influence of the momentum exchange from the gas phase to solid phase is sufficiently small to not fluidise the solid. As shown in Fig. 2b, the discretised cell includes cubic structural lattices, and the volume of the solid phase particle is embedded in the cell by defining the cell with the centre of the particle as “filled”. The velocity on the boundary between cells occupied by the solid and gas phases was set zero as the boundary condition. Although this simple classical boundary condition has lower resolution compared to the body-fitted mesh, it has the advantage of faster computation and has been applied for the geometric analysis detailed in section 4. The momentum exchange from gas to solid was not considered by this “1way” scheme gas phase flow passing through the packed bed. Therefore, it was assumed that the influence of the gas phase flow on the softening solid deformation is sufficiently small.

2.3. Gas-flow simulation In order to evaluate the gas permeability of the packed bed structure obtained using the numerical scheme described in the previous section, gas-flow simulation reflecting the gas-solid boundary was used as per the Eulerian approach. The local void shape can be obtained from the solid locations. The space was discretized into square lattices and each cell was occupied by a gas or solid. In the extensively used DEM-CFD method, the void fraction and mean-particle diameter information are generally provided for each cell; however, in this study, a number of solid shapes must be reflected accurately to handle solid deformation. Here, the gas flow was obtained by solving the continuity equation and the Navier–Stokes equation assuming incompressible flow and constant properties: (22)

u = 0. u + t

u u=

1

P+

µ

2u

2.4. 3D capturing of the solid shape (23)

Considering the particle size distribution immediately above the blast furnace raceway, sieving was performed using mesh sizes of 15 mm and 35 mm to obtain the representative coke and ore samples (indexed from C1–C10 and O1–O5, respectively). Using a 3D laser scanner (Matter and Form), the surface points of each sample were obtained. Approximately 300,000 surface points for each coke sample were obtained with a minimum resolution of 0.43 mm. The obtained coordinates were converted into standard triangulated language (STL), and the surface shape was polygonate with a triangular mesh. The results of the representative 3D shape of each sample obtained by this

Here, u and P are the gas velocity and pressure, respectively. Considering the container width of the calculation domain, L, as the representative length scale and the inflow velocity, u 0, as the representative velocity scale, the Reynolds number, Re = u 0 L µ . Here, the dimensionless pressure was obtained using P0 = u 02 . The finite difference method, using a staggered grid, was applied with the Cartesian coordinate system, and Eq. (23) was discretised based on the unsteady simplified marker and cell (SMAC) method [47]. The nonsteady, advective, pressure gradient, and viscosity terms were 8

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 11. Calculated gas velocity distributions in the vertical cross sections for each softened-ore condition.

method are shown in Fig. 4. In order to apply this shape information for MSDEM and ADEM simulation, the spherical calculation particles were filled within the scanned surface shells of coke and ore to conform to the non-spherical shape of each sample. We employed a level set function that is related to the distance from the centre position of each spherical particle arranged in a cubic lattice at an interval, dp, to the surface of the triangular polygon embedded in the volume [50]. In this study, dp = 1.00 mm was adopted as a constant value because of the trade-off between the computational cost and analysis accuracy [13]. Fig. 5 displays ore sample #O3, the STL polygon image, and particle assembly image; the three surface shapes are in very good agreement with the actual sample shape.

and packing simulation into a rectangular container with a 0.10-m-side square cross-section was performed. The 50 pieces of ore were then charged similarly to form a layered structure. Coke samples with a density of 1,050 kg/m3, dynamic friction coefficient of 0.43, Young's modulus of 5.4 GPa, and Poisson’s ratio of 0.22, and ore samples with a density of 3,950 kg/m3, dynamic friction coefficient of 0.43, Young's modulus of 35.0 GPa, and Poisson’s ratio of 0.24 were employed [51]. Although the values of the Young's modulus and Poisson’s ratio are directly reflected in the spring constant in the contact force model, they are insufficient to specify the high-temperature elastic mechanical properties of the ore under various degrees of reduction. Therefore, in order to analyse the influence of the softening degree on the packed-bed structure, the joint spring constant was set to ka < 1.0 N/m, considering the magnitude of the softened ore 's elastic force with respect to the gravitational force, which contributes to deformation (see definition of Sf number in Appendix A). In this study, the connections between the connected primary particles do not break, even if the joint spring extends; hence, the deformation alone was considered. When softening progresses, it is conceivable that the solid will eventually melt, and forces such as the liquid crosslinking force will act. The phase change from solid to liquid and the breakage of the connection of primary

2.5. Analytical condition Using the particle assembly data obtained in the previous section, 150 pieces (100 cokes, 50 ores) of “solid” were generated, including 10 pieces of 15 representative samples. The total number of calculation particles including the wall particles was 1,351,068. To obtain the initial-packed bed structure, the position and rotation angle of each of the 100 coke pieces were first arranged using a pseudo random number, 9

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 12. Isobar planes in a packed bed with different joint spring constants. P was normalised by Pmax atka = 0.001 N/m.

particles will be addressed in future. Fig. 6 depicts the packing behaviour. The layered bed of coke and ore was successfully simulated. Similar packing tests were repeated five times, and the average projected area of each sample, Aproj , and the volume equivalent diameter, d v [12], were calculated. Carman's shape factor, (= d v 2 4Aproj ), which is one of the indices representing the sphericity [52], was obtained from the 3D surface data of all the coke samples, and summarised in Table 1. To evaluate the gas flow field, cubic meshes were generated in the packed bed space and a buffer space below the packed bed was installed for stabilizing the gas flow. We assumed a cell length, x = y = z = 1.000 mm, t = 10-7 s, maximum cell number, nx:ny:nz = 100:100:400, and inflow velocity, u 0 = 0.1 m/s.

penetration into the coke bed at k’ < 150 in Fig. 6. Consequently, we derive two the deformation mechanisms of ore bed as follows: (i) The exclusion of the void space among the ore particles by the deformed ore particles in the initial ore layer and (ii) penetration of the fused ore particles into the void space in the coke bed in a more softened condition. With these two void reduction effects in the ore layer in range (A) and the softened-ore penetration into the coke bed in range (C), the ore layer thickness, hmax hmin, is the least in the range, 0.01 < ka < 0.1 N/m. As per the vertical z-direction variation, the densest ore layer should be formed in this range. Fig. 9 shows the distribution of ore and coke in the horizontal cross-section x-y plane at various heights. At k’ = 180, in the range, 0.01 < ka < 0.1 N/m, small voids exist dispersedly in the ore layer, but at ka = 0.001 N/m, these voids disappear because the layer is completely occupied by softened ore, whereas the void area near the wall is enlarged. At k’ = 160, penetration of the softened ore into the coke bed is noticeable when ka < 0.1 N/m, and voids also exist in the vicinity of the coke surfaces. Fig. 10 shows the effect of ka on the void fraction in each arbitrary packed-bed structure; here, the changes in the void fraction of region ‘a’, where cokes and ores were present in the initial arrangement, and in the divided regions b–e are depicted. The differences were caused by the thickness of the coke layer up to the ore bed. The void fraction changes depicted in the figure are due to the presence of voids in the ore layer in the region where ka > 0.1N m , and the penetration effect of the fused ore into the coke bed in the region where ka < 0.1N m . As shown in Fig. 10(b), the bed comprising ore alone easily clogs with the progress of slight softening. Fig. 10(c) depicts approximately half the highest first coke layer, and it is obvious that the voids around coke remain, even in the ore softening condition. In the lower part of the coke bed, ore penetration occurs in the later stage of softening (Fig. 10(d) and (e)).

3. Results and discussion 3.1. Dynamic behaviour of the softening ore Fig. 7 shows the variation of the softening-ore-bed static morphology with the joint spring constant (herein referred to as ka ). The deformation of the ore increases with the gravitational force (along the z-axis direction) as ka decreases. The vertical cross-sections indicate that the void shape in the softened ore bed markedly decreases down to ka = 0.1 N/m, and penetration of the softened ore into the coke packed layer occurs with the further decrease in ka . The z-direction deformation profile is quantitatively depicted in Fig. 8. The maximum height, hmax , and the minimum height, hmin , of each ka were obtained, and the difference, hmax hmin , was considered as the representative thickness of the layer. The maximum height hmax decreases with decrease in ka in the range, 0.1 < ka < 1 N/m (range (A) in Fig. 8), and shows negligible variation with further decrease in ka . This can also be established by comparison with the height, k’ = 200, in Fig. 7. Here, the cell indices in the x, y, and z directions are represented by i', j' , k ' (1 i' nx , 1 j' ny , and 1 k ' nz ). On the other hand, hmin decreases in the range, 0.001 < ka < 0.01 N/m ((C) in Fig. 8), corresponding to the

3.2. Simulation of the gas flow in a packed bed containing softening ore In this section, we discuss the gas flow behaviour in a packed bed 10

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 13. Calculated gas velocity distributions in the horizontal cross sections for each softened-ore condition.

containing softening ore. Fig. 11 shows the gas velocity distribution in each vertical section. In this figure, the packed-bed structure at each joint spring constant, ka, corresponds to that of Fig. 6. Although the gas velocity field in the coke bed (region k’ (1 0 0) shows negligible change, under any ka value, the flow in the region occupied by the softening ore is strongly suppressed under small ka , while the velocity near the wall (j’ = 1) increases. In addition, as ka decreases, the local flow velocity increase in the ore layer within the packed bed becomes significant. Fig. 12 depicts the pressure distribution for each value of ka . Reflecting the shapes of the solids and the void space among them, the isobar surfaces are irregular. Generally, the pressure-drop increases as ka decreases, but in the case of ka = 0.01 N/m, the total pressure-drop is lower than that at ka = 0.05 N/m, which is considered to be the effect of the void distribution, including near the wall region, as discussed in the previous section. This suggests that there is an indeterminacy due to the shape of the void. To predict the pressure-drop in the packed bed in which the void is present, it is necessary to consider not only the void but also the influence of the shape of the void [53]. This indicates that in a packed bed with a dispersed phase (softened ore, in this case), all the voids are not present equivalently, for breathability. Fig. 13 shows the gas flow velocity distribution in each horizontal cross-section. This figure corresponds to the phase map in Fig. 9. In Fig. 13, in the whitecoloured region, u = 0 m/s. A major part of the white region is occupied by the solid. For ka = 0.05 N/m, small voids are dispersed over the cross- section at k’ = 160 in Fig. 9, but in Fig. 13, these voids are

present at u = 0 m/s. Hence, almost no gas circulates in the remaining voids in the ore bed, which are generated as softening occurs. The gas flow in the packed bed should be influenced by the upstream flow structure and a specific intensive flow path must be formed. However, the existence of voids “not” contributing to the gas flow renders the analytical correlation between the pressure-drop and ka complicated. Therefore, to correlate them in the next section, we attempt to geometrically analyse the relationship between the continuity of the local voids and the pressure loss. 4. Evaluation of the gas permeability 4.1. Pressure drop estimation method by multi-diameter circular tube approximation The pressure drop in the packed bed was estimated from the geometrical data of the calculated softening packed-bed structure. In this study, the void space in the packed bed was assumed to have a continuous structure of multi-diameter circular tubes (MDCT) with varying diameters in the flow direction, and the local pressure drop was calculated using the Darcy–Weisbach equation:

P=

l D

u2 . 2

where P is the pressure drop, 11

(24) is the Darcy friction factor, l is the tube

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 14. Schematic illustration of the effective-void decision.

length, D is the tube diameter, is the fluid density, and u is the gas flow velocity. The local P is integrated over the packed bed height to estimate the overall P . To calculate the pressure drop using the Darcy–Weisbach equation, the equivalent hydraulic diameter, Dh , was used as the tube diameter.

Dh =

4A . L

veff (i', j' , k ') =

c (i ' , j ' , k ' ) =

(25)

'

veff (i , j , 1) = 1.

j' + 1 veff (a, b, k ' b = j' 1

i'+ 1 a = i' 1 '

'

1) = 0

'

.

or [c (i , j , k ) = 0] '

1;

where the flow area, A, is the cross-sectional area of the gas flow and the wetted perimeter, L, is the length of the line on which the fluid and solid contact in the cross-section perpendicular to the main gas flow direction. In order to obtain D h from the flow area, the wetted perimeter in each cross-section was obtained by the following method. As shown in Fig. 13, within the deforming packed bed, voids in which no gas flows exist. Among the voids in the packed bed, those in the flow-path were defined as effective voids, whereas those that did not contribute to the flow were defined as ineffective voids. The void classification procedure is schematically shown in Fig. 14. Initially, all the computational cells were classified as solid cells and void cells. The void cells were further classified as effective and ineffective cells. This classification was made with reference to the state of the upstream cells, i.e., the continuity of the effective cell was considered. In this study, the main flow direction was set as the positive z-direction; thereby, classification was made using the following equations. '

c (i', j' , k ') = 1 and

0;

'

'

c (i , j , k ) = 1 and

i'+ 1 a = i' 1

0(solid) . 1(fluid)

j' + 1 veff (a, b, k ' b = j' 1

1)

0

(27)

(28)

where veff = 1 is an effective void and veff = 0 is an ineffective one. In the flow simulation, a buffer zone was created beneath the packed bed and above the inlet boundary. The cells in this zone were set as the effective voids. For the other cells, effective/ineffective-void determination was performed using Eqs. (27) and (28) shows whether each cell is occupied by a solid or fluid, and is defined using the calculated packed-bed structure. The flow area, A(k’), on the x-y plane at k’ was defined as the total projection area of the effective void cell, and expressed by the following Eq. (29): (29)

A (k ') = n void (k ') l c 2. n void (k') =

nx

ny

i '= 1

j '= 1

veff (i ', j', k ').

(30)

The wetted perimeter, L (k') , was set as the sum of the cell sidelength between the effective void and solid cells, l wet (i', j', k') , on the x-y plane.

(26) 12

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 15. Calculation result of the cross-sectional area-ratio change in each phase by ore softening.

L (k ) =

nx

ny

i'= 1

j' = 1

l wet (i', j' , k ').

(31)

D h (k ' ) =

l wet (i', j', k') was calculated using the number of effective voids contacting with the solid cell and the lattice constant, lc .

l wet (i', j', k') =

0; c (i', j', k') = 1 nc lc ; c (i ', j', k ') = 0

nc = veff (i '+ 1, j', k') + veff (i ' (i ', j '

1, k ').

4A (k ') L (k ' )

.

(34)

In this study, the Darcy friction factor, λ, was obtained as per Eq. (35):

=

(32)

1, j', k') + veff (i', j '+ 1, k') + veff

64 Re'

Re ' =

(33)

uD h ( k ' ) . µ

(35)

(36)

The flow velocity, u , was defined as the average velocity in the effective void, and expressed by the following equation:

Substituting A (k ') and L (k') in Eq. (25), the equivalent hydraulic diameter, D h (k'), in the packed bed was obtained. 13

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 16. Change in the ore-volume ratio distribution due to softening.

u = u0

nx ny n void (k')

.

Fig. 14, with the progress of softening, the ore distribution moves downwards, confirming that the height at which the maximum area ratio occurs hardly changes. The softening from the initial condition to ka = 0.05 N/m mainly exhibits a tendency in which the upper ore shifts to the lower part of the ore layer. On the other hand, when softening proceeds beyond ka = 0.05 N/m, the ores penetrate into the coke layer, and at ka = 0.05 N/m, the ore distribution is the narrowest. This may be the reason for the reduction of the ineffective void, when ka changes from 0.05 to 0.01 N/m. The comparison between the pressure-drop obtained by gas-flow simulation and MDCT approximation is depicted in Fig. 17. The pressure-drops shown in this figure are respectively normalized by their maximum values in each condition and method. This figure also shows the average cross-sectional flow velocity, u, and the equivalent hydraulic diameter, D h , obtained by MDCT approximation. For any joint spring coefficient, the location of the maximum pressure-drop obtained by MDCT approximation coincides with the one obtained by gas-flow simulation. Additionally, for the non-softening condition, MDCT approximation reproduces the periodical change of the pressure-drop and its variations indicating that MDCT approximation as well as gas-flow simulation can reflect the periodic change in the packed-bed structure. However, it should be noted that the pressure profiles obtained by MDCT approximation and simulation differ as ka becomes smaller because of the complex flow field in the softening ore layer. In future, this MDCT model may be extended to adapt the complex flow field for the accurate prediction of the effective-void judgment, which refers to the 3-dimensional velocity. Fig. 18 shows the relationship between the pressure-drop throughout the entire packed bed and the joint spring coefficient. For

(37)

Substituting Eqs. (34)–(37) in Eq. (24), the local pressure-drop at k' in the packed bed was obtained:

P (k ' ) =

32lc µ u 0

nx ny

D h (k ')2 n void (k ')

.

(38)

4.2. Pressure-drop distribution obtained using the MDCT approximation model and comparison with the simulation results Fig. 15 shows the vertical (z-direction) distribution of the material and void fractions calculated using the procedure explained in the previous section. In the initial condition, the area ratio of the ore and coke changes periodically with the solid diameter. This periodic variation in the ore-layer disappears, when ka decreases to 1.0 N/m. Moreover, the height range in which the ore exists becomes narrow, and the ore area ratio increases in this range. With the decrease in ka , the distribution of ore area moves downwards. The movement of the lower boundary of the ore layer is less with the decrease in ka down to 0.05 N/m; additionally, the maximum solid area fraction in the oreexisting area increases. On the contrary, the existing region of the ore decreases considerably, with the further decrease in ka . This reflects the penetration of the softened ore into the coke bed. The maximum solid ratio shows negligible change. The ore distribution is shown in Fig. 16. The ore volume ratios are plotted against the relative height of which the origin is set at the height the ore area ratio is the maximum in the initial condition. As shown in 14

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Fig. 17. Comparison of the pressure distribution in the packed bed between the numerical simulation value and analytical results of the MDCT approximation.

comparison, the pressure-drops obtained by gas-flow simulation and MDCT approximation are divided by the pressure-drop at the initial condition; i.e., the pressure- drop in each packed bed is normalised with reference to a packed bed without deformation. With the progress of softening, the pressure-drop initially increases in the range of the joint spring coefficient down to 0.05 N m−1. Further, the pressure-drop slightly decreases. When the joint spring coefficient is less than 0.01 N m−1, the pressure-drop increases again. Although MDCT approximation provides a higher increment of the pressure-drop than the gas-flow simulation, the variations in the pressure-drop agree well. This tendency can be explained as follows. In the initial stages of softening (range (A) in Fig. 18), the ore layer above the coke layer collapses. The void space in the ore bed becomes smaller and narrower as shown in Figs. 8, 15–17. Consequently, the pressure-drop increases with the progress of softening. In range (B), rearrangement of the softened ores occurs, and the ineffective void at 0.170 < z < 0.190m decreases from ka = 0.5 to 0.05 N/m. On the other hand, the total void space slightly increases from ka = 0.1 to 0.01 N/m due to the slight overall ore

descent, as shown in Fig. 15. This change also decreases the pressuredrop slightly. In range (C), the softened ore deeply penetrates into the coke bed as shown in Figs. 8, 15 and 16; this penetration eliminates the void space in the coke bed and increases the pressure drop. 5. Conclusion A new simulation model for predicting gas-clogging behaviour was developed, based on the thermal deformation modelling of reduced softened ores in the ironmaking process. The complex incompressible gas-flow in packed beds, affected by softened-ore gravitational deformation, was successfully reflected in the Eulerian–Lagrangian coupling numerical scheme. The numerical results and geometrical analysis of the flow path in the packed bed were compared. The following conclusions were obtained: I. The ores were deformed by gravity due to softening; consequently, the two-step mechanism of the ore bed rearrangement was derived 15

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

structure and formed a specific path. Even after ore softening, the main flow path did not change; however, when the flow in the vertical direction reached the critical-pressure value, path branching occurred. III. Using geometric analysis, the influence of the change in the packedstructure due to the deformation of the softened ore on the shape of the void and the pressure distribution were investigated. The pressure distribution in the packed-bed obtained using the multidiameter circular tube approximation model agreed well with that obtained by E-L simulation. These results suggest that the formation of ineffective voids in the gas-flow path increase the local pressuredrop as softening progresses. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments We would especially like to thank Prof. R. O. Suzuki, Dr. T. Kikuchi, and the eco-processing laboratory members in Hokkaido university, where S. Natsui had conducted this research. S. Natsui was partially supported by the Iron and Steel Institute of Japan (ISIJ) Research Promotion Grant, Amano Institute of Technology Research Promotion Grant, Tanikawa Fund Promotion of Thermal Technology Research Promotion Grant, and the Steel Foundation for Environmental Protection Technology (SEPT). S. Ishihara was partially supported by the Grants-in Aid for Scientific Research “KAKENHI” (17K14846) of Japan Society for the Promotion of Science, and Ishihara-Asada Research Fund through the Iron and Steel Institute of Japan. K. Ohno and S. Natsui were partially supported by the Grants-in-Aid for Scientific Research “KAKENHI” (15H04168) of Japan Society for the Promotion of Science.

Fig. 18. Relationship between the normalised pressure drop and joint spring coefficient, ka .

as follows: In the first step, the individual ore shape was lost by deformation in the initial ore layer, and in the second, the fused ore particles penetrated into the void space in the coke bed in a more softened condition. The softened ores occupied the voids in the orelayer. Isolated voids did not contribute to the gas flow. As softening proceeded, the ore penetrated the coke layer. II. The gas flow in the packed bed was influenced by the upstream bed

Appendix A. Basic dimensional analysis of the gravitational deformation of softening material Let us consider generalising the amount of deformation in an elastic body, which is acted upon only by gravity. If described from the deformation of the softening material to breakage, we must consider the correlation between the hydrostatic pressure and yield stress of the softening material. In order to avoid considering the yield stress of the material because of the difficulty in presenting it parametrically, only a single solid sphere is assumed here for simplification. It is simpler to obtain the deformation between individual joint springs as the equivalent deformation of the solid and consider only one dimension in the gravity direction. Supposing that x is the shrinkage from the natural length of the individual spring and k is n 1 the spring constant in the gravitational direction, the elastic energy, Us = 2 kx 2 , is balanced by an initial potential energy, Ug = m gh , and X = i = 1 x i is the amount of shrinkage of all the springs in the system; X can be expressed as

X=

2m gh . k

(A1)

We define the ratio between the gravitational force and elastic force acting on the softening material as the dimensionless degree of deformation:

Sf0 =

mg = kX

mg 2km gh

.

(A2)

Assuming h = dp and m =

Sf0 2 =

gdp 2 12k

dp3 6

, we get

.

(A3)

Thus, we define Sf again as the dimensionless single-sphere deformation degree due to softening:

Sf =

gdp2 k

.

Therefore, because Sf =

(A4) 12

Sf0 2 , this simplification can describe the following mode change of the deformation as

16

Chemical Engineering Journal xxx (xxxx) xxxx

S. Natsui, et al.

Sf <

12

: Elastic force domination, gravity can be ignored . Sf =

12

12

: Elastic force and gravity balance.

.

< Sf : Gravity domination, softening material is deformed. (A5)

k

Given the average ore diameter, dp = 0.023m, and its density, 5.4N m .

= 3950kg m3, the spring coefficient to balance the gravity is estimated as

[28] F. Municchi, S. Radl, Consistent closures for Euler-Lagrange models of bi-disperse gas-particle suspensions derived from particle-resolved direct numerical simulations, Int. J. Heat Mass Transf. 111 (2017) 171–190. [29] R. Kasper, J. Turnow, N. Kornev, Numerical modeling and simulation of particulate fouling of structured heat transfer surfaces using a multiphase Euler-Lagrange approach, Int. J. Heat Mass Transf. 115 (2017) 932–945. [30] H. Wu, N. Gui, X. Yang, J. Tu, S. Jiang, Numerical simulation of heat transfer in packed pebble beds: CFD-DEM coupled with particle thermal radiation, Int. J. Heat Mass Transf. 110 (2017) 393–405. [31] S. Bellan, K. Matsubara, H.S. Cho, N. Gokon, T. Kodama, A CFD-DEM study of hydrodynamics with heat transfer in a gas-solid fluidized bed reactor for solar thermal applications, Int. J. Heat Mass Transf. 116 (2018) 377–392. [32] H. Wu, N. Gui, X. Yang, J. Tu, S. Jiang, A smoothed void fraction method for CFDDEM simulation of packed pebble beds with particle thermal radiation, Int. J. Heat Mass Transf. 118 (2018) 275–288. [33] M. Trofa, G. D’Avino, L. Sicignano, G. Tomaiuolo, F. Greco, P.L. Maffettone, S. Guido, CFD-DEM simulations of particulate fouling in microchannels, Chem. Eng. J. 358 (2019) 91–100. [34] C. Moliner, F. Marchelli, N. Spanachi, A. Martinez-Felipe, B. Bosio, and E. Arato, CFD simulation of a spouted bed: Comparison between the Discrete Element Method (DEM) and the Two Fluid Model (TFM), Chemical Engineering Journal, in press. [35] M. Maestri, G. Salierno, S. Piovano, M. Cassanello, M. A. Cardona, D. Hojman, and H. Somacal, CFD-DEM modeling of solid motion in a water-calcium alginate fluidized column and its comparison with results from radioactive particle tracking, Chemical Engineering Journal, in press. [36] S. Bellan, T. Kodama, K. Matsubara, N. Gokon, H. S. Cho, and K. Inoue, Thermal performance of a 30 kW fluidized bed reactor for solar gasification: A CFD-DEM study, Chemical Engineering Journal, in press. [37] J.P. Latham, A. Munjiza, X. Garcia, J. Xiang, R. Guises, Three-dimensional particle shape acquisition and use of shape library for DEM and FEM/DEM simulation, Miner. Eng. 21 (11) (2008) 797–805. [38] H.S. Rabbani, V. Joekar-Niasar, T. Pak, N. Shokri, New insights on the complex dynamics of two-phase flow in porous media under intermediate-wet conditions, Sci. Rep. 7 (1) (2017) 4584. [39] S.P. Sullivan, F.M. Sani, M.L. Johns, L.F. Gladden, Simulation of packed bed reactors using lattice Boltzmann methods, Chem. Eng. Sci. 60 (12) (2005) 3405–3418. [40] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assembles, Geotechnique 29 (1979) 47–65. [41] W. Zhong, A. Yu, X. Liu, Z. Tong, H. Zhang, DEM/CFD-DEM modelling of nonspherical particulate systems: theoretical developments and applications, Powder Technol. 302 (2016) 108–152. [42] S. Natsui, S. Ueda, H. Nogami, J. Kano, R. Inoue, T. Ariyama, Gas–solid flow simulation of fines clogging a packed bed using DEM–CFD, Chem. Eng. Sci. 71 (2012) 274–282. [43] H. Gotoh, A. Okayasu, Y. Watanabe, Computational wave dynamics, World Scientific Publishing Co., Singapore, 2013, pp. 137–194. [44] J.F. Favier, M.H. Abbaspour-Fard, M. Kremmer, A.O. Raji, Shape representation of axi-symmetrical, non-spherical particles in discrete element simulation using multielement model particles, Eng. Computations 16 (4) (1999) 467–480. [45] H. Kruggel-Emden, S. Rickelt, S. Wirtz, V. Scherer, A study on the validity of the multi-sphere Discrete Element Method, Powder Technol. 188 (2008) 153–165. [46] S. Ishihara, J. Kano, ADEM simulation for analysis of particle breakage of irregular shaped particles, ISIJ Int. (2019) in press. [47] A.A. Amsden, F.H. Harlow, A simplified MAC technique for incompressible fluid flow calculations, J. Comput. Phys. 6 (1970) 322–325. [48] B.P. Leonard, A survey of finite differences with upwinding for numerical modeling of the incompressible convective diffusion equation. Computational Techniques in Transient and Turbulent Flow, vol.2, edited by Taylor, C. and Morgan, K., Pineridge Press, Swansea, U. K., (1981), 1-35. [49] H.A. Van der Vorst, Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 13 (2) (1992) 631–644. [50] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146–159. [51] S. Natsui, R. Shibasaki, T. Kon, S. Ueda, R. Inoue, T. Ariyama, Effect of high reactivity coke for mixed charge in ore layer on reaction behavior of each particle in blast furnace, ISIJ Int. 53 (10) (2013) 1770–1778. [52] P.C. Carman, Fluid flow through porous rock, Trans. Inst. Chem. Eng. London 15 (1937) 150–157. [53] S. Natsui, A. Sawada, K. Terui, Y. Kashihara, T. Kikuchi, R.O. Suzuki, Evaluation of coke degradation effect on flow characteristics in packed bed using 3D scanning for rotational cechanical strength test and solid-liquid-gas three-phase dynamic model analysis, Tetsu-to-Hagané 104 (7) (2018) 347–357.

References [1] T. Ariyama, M. Sato, T. Nouchi, K. Takahashi, Evolution of blast furnace process toward reductant flexibility and carbon dioxide mitigation in steel works, ISIJ Int. 56 (10) (2016) 1681–1696. [2] J. Ishii, R. Murai, I. Sumi, Y. Yongxiang, R. Boom, Gas permeability in cohesive zone in the ironmaking blast furnace, ISIJ Int. 57 (9) (2017) 1531–1536. [3] P.F. Nogueira, R.J. Fruehan, Blast furnace burden softening and melting phenomena: Part I. Pellet bulk interaction observation, Metall. Mater. Trans. B 35 (5) (2004) 829–838. [4] A. Kemppainen, K.I. Ohno, M. Iljana, O. Mattila, T. Paananen, E.P. Heikkinen, T. Maeda, K. Kunitomo, T. Fabritius, Softening behaviors of acid and olivine fluxed iron ore pellets in the cohesive zone of a blast furnace, ISIJ Int. 55 (10) (2015) 2039–2046. [5] M. Baniasadi, B. Peters, J.C. Pierret, B. Vanderheyden, O. Ansseau, Experimental and numerical investigation into the softening behavior of a packed bed of iron ore pellets, Powder Technol. 339 (2018) 863–871. [6] X.F. Dong, A.B. Yu, S.J. Chew, P. Zulli, Modeling of blast furnace with layered cohesive zone, Metall. Mater. Trans. B 41 (2) (2010) 330–349. [7] S. Ueda, T. Kon, H. Kurosawa, S. Natsui, T. Ariyama, H. Nogami, Influence of shape of cohesive zone on gas flow and permeability in the blast furnace analyzed by DEM-CFD model, ISIJ Int. 55 (6) (2015) 1232–1236. [8] P.R. Austin, H. Nogami, J. Yagi, A mathematical model for blast furnace reaction analysis based on the four fluid model, ISIJ Int. 37 (8) (1997) 748–755. [9] S. Ueda, S. Natsui, H. Nogami, J.I. Yagi, T. Ariyama, Recent progress and future perspective on mathematical modeling of blast furnace, ISIJ Int. 50 (7) (2010) 914–923. [10] T. Kon, S. Natsui, S. Ueda, H. Nogami, Analysis of effect of packed bed structure on liquid flow in packed bed using moving particle semi-implicit method, ISIJ Int. 55 (6) (2015) 1284–1290. [11] M. Tsuboi, K. Ito, Cold Model Experiment and Numerical Simulation of Flow Characteristics of Multi-phase Slag, ISIJ Int. 57 (7) (2017) 1191–1196. [12] S. Natsui, H. Takai, R. Nashimoto, K.I. Ohno, S. Sukenaga, T. Kikuchi, R.O. Suzuki, Capturing the non-spherical shape of granular media and its trickle flow characteristics using fully-Lagrangian method, AIChE J. 63 (6) (2017) 2257–2271. [13] S. Natsui, R. Nashimoto, T. Kumagai, T. Kikuchi, R.O. Suzuki, An SPH Study of Molten Matte-Slag Dispersion, Metall. Mater. Trans. B 48 (3) (2017) 1792–1806. [14] S. Natsui, A. Sawada, K. Terui, Y. Kashihara, T. Kikuchi, R.O. Suzuki, DEM-SPH study of molten slag trickle flow in coke bed, Chem. Eng. Sci. 175 (2018) 25–39. [15] S. Natsui, K.I. Ohno, S. Sukenaga, T. Kikuchi, R.O. Suzuki, Detailed modeling of melt dripping in coke bed by DEM–SPH, ISIJ Int. 58 (2) (2018) 282–291. [16] F. Bambauer, S. Wirtz, V. Scherer, H. Bartusch, Transient DEM-CFD simulation of solid and fluid flow in a three- dimensional blast furnace model, Powder Technol. 334 (2018) 53–64. [17] S. Kuang, Z. Li, A.B. Yu, Review on modeling and simulation of blast furnace, Steel Res. Int. 89 (1) (2018) 1700071. [18] A. Khazeni, Z. Mansourpour, Influence of non-spherical shape approximation on DEM simulation accuracy by multi-sphere method, Powder Technol. 332 (2018) 265–278. [19] R. Cabiscol, J.H. Finke, A. Kwade, Calibration and interpretation of DEM parameters for simulations of cylindrical tablets with multi-sphere approach, Powder Technol. 327 (2018) 232–245. [20] Y. He, T.J. Evans, Y.S. Shen, A.B. Yu, R.Y. Yang, Discrete modelling of the compaction of non-spherical particles using a multi-sphere approach, Miner. Eng. 117 (2018) 108–116. [21] Y. You, M. Liu, H. Ma, L. Xu, B. Liu, Y. Shao, Y. Tang, Y. Zhao, Investigation of the vibration sorting of non-spherical particles based on DEM simulation, Powder Technol. 325 (2018) 316–332. [22] S. Das, N.G. Deen, J.A.M. Kuipers, Multiscale modeling of fixed-bed reactors with porous (open-cell foam) non-spherical particles: Hydrodynamics, Chem. Eng. J. 334 (2018) 741–759. [23] B. Soltanbeigi, A. Podlozhnyuk, S.A. Papanicolopulos, C. Kloss, S. Pirker, J.Y. Ooi, DEM study of mechanical characteristics of multi-spherical and superquadric particles at micro and macro scales, Powder Technol. 329 (2018) 288–303. [24] S. Ishihara, Q. Zhang, J. Kano, ADEM simulation of particle breakage behavior, J. Soc. Powder Technol., Japan 51 (6) (2014) 407–414. [25] L. Kempton, D. Pinson, S. Chew, P. Zulli, A. Yu, Simulation of particle deformation within a compressed packed bed, Powder Technol. 320 (2017) 586–593. [26] S. Natsui, T. Kikuchi, R.O. Suzuki, Numerical analysis of carbon monoxide–hydrogen gas reduction of iron ore in a packed bed by an Euler-Lagrange approach, Metall. Mater. Trans. B 45 (6) (2014) 2395–2413. [27] S. Natsui, H. Takai, R. Nashimoto, T. Kikuchi, R.O. Suzuki, Model study of the effect of particles structure on the heat and mass transfer through the packed bed in ironmaking blast furnace, Int. J. Heat Mass Transf. 91 (2015) 1176–1186.

17