ISIJ International
Online ISSN : 1347-5460
Print ISSN : 0915-1559
ISSN-L : 0915-1559
Special Issue on "Recent Approaches to Control of Cohesive Zone Phenomena and Improvement of Permeability in Blast Furnace"
Numerical Simulation of Coexisting Solid-liquid Slag Trickle Flow in a Coke Bed by the SPH Method with a Non-Newtonian Fluid Model
Shungo Natsui Akinori SawadaHiroshi NogamiTatsuya KikuchiRyosuke O. Suzuki
Author information
JOURNAL OPEN ACCESS FULL-TEXT HTML

2020 Volume 60 Issue 7 Pages 1445-1452

Details
Abstract

A simulation of the dripping behavior of solid-phase suspended molten slag in a coke bed was performed using the Bingham fluid model in the framework of smoothed particle hydrodynamics (SPH). SPH can track the motion of droplets which are suspended in dispersed solids. A case study was performed in which the yield coefficient varied with the viscosity coefficient of the molten SiO2–CaO–Al2O3 slag, with suspended solids approximated by the Bingham fluid model, and the holdup droplets were trapped at approximately the same site regardless of the yield value. When the yield value exceeded the threshold, the volume of each holdup droplet increased. This threshold is correlated with the decrease in the shear rate at the bottleneck. Increasing the solid phase ratio in the molten slag is also predicted to increase the yield value with viscosity; thus, the amount of holdup droplets increases at the specific holdup sites as a dispersed phase. This indicates that the increase in the yield value approaches the powder holdup mechanism that is closed one after another starting from a specific holdup site.

1. Introduction

In the blast furnace ironmaking process, there is an increasing demand to decrease the coke ratio because a large amount of CO2 is released from the ironmaking processes. When the volume of the coke, which acts as a spacer in the lower part of a blast furnace, is reduced, the resistance to the melt passage increases, which can cause instability. In other words, an operational design based on a phenomenological understanding of the slag drip behavior in the coke bed is important, and there have been numerous studies on such melt properties.1,2,3,4) Therefore, an accurate understanding of melt behavior characteristics is vital for the construction of a low-carbon blast furnace.

Typically, when the liquid flow through a packed bed is in a steady state, the volume fraction of the liquid in the packed bed is called “liquid hold up.” Liquid hold up, ht, is composed of static, hs, and dynamic (operating), hd holdups.5) A static holdup is the volume fraction of the liquid in the packed bed that remains in the bed when the supply of liquid is stopped, and it is much larger than the dynamic holdup, which is a dominant factor in determining ht. It is considered that hs should be determined by the balance between the gravitational force and the liquid interfacial force. Because molten slag has a higher density, viscosity, and surface tension than ordinary temperature liquid phase, it exhibits a different holdup behavior. Therefore, it would be difficult to predict the holdup behavior of high-temperature melts with results from an analysis that uses an experimental system limited to room temperature media.

In addition, when powder is present in the slag, the dripping zone will exhibit more complex behavior. When coke descends into the lower part of a blast furnace, the volume destruction and surface wear progress due to collisions between the furnace wall and the coke,6) and the generated coke powder comes in contact with the melt flowing down into the lower part of the blast furnace. Even if the coke degradation doesn’t occur, then oxides that have high melting points, such as 2CaO·SiO2, are thermodynamically liberated and become suspended in the molten slag.7,8) Consequently, the powder suspended in the slag phase transitions from being a Newtonian fluid, in which viscosity is proportional to shear stress and flow shear rate, and becomes a non-Newtonian fluid called a “Bingham plastic fluid” with a concomitant increase in the solid phase ratio.9,10)

Although the interaction between the melt and powder should not be ignored, few reports focus on the complicated dynamics of packed-bed behavior in detail.11) Numerical simulation analysis using the discrete element type model helps describe such dispersed-powder motion, and in recent years, such analyses have been applied to local void blockages where the liquid phase is mixed around the powder.12) It is also possible to analyze the liquid holdup in the packed bed by using smoothed particle hydrodynamics (SPH), a method that can analyze the dispersed-liquid phase seamlessly from the continuous phase.13,14) Application of the “coarse-graining” approach for SPH to non-Newtonian fluids have been reported, and it has become possible to analyze actual systems that cause engineering problems,15,16,17,18) but detailed reports on the trickle-flow behavior in the packed bed of a Bingham fluid with a free surface. In this study, a trickle-flow simulation in a coke-packed bed for a Bingham fluid was conducted by the SPH method to investigate the flow behavior and holdup characteristics.

2. Numerical Scheme

2.1. SPH Procedure for Fluid Motion

SPH is based on an interpolation scheme that allows the estimation of a vector or scalar function f at position r in terms of the values of the function at the discretization points.   

f( r ) f( r ) W( r,h ) dV , (1)
where W is the smoothing kernel function with a length of influence radius h. We chose Wendland’s kernel to avoid kernel artifacts such as particle clustering.13) By using the divergence theorem, the kernel function W(r, h)=0 on the boundary: therefore, the gradient form of Eq. (1) is represented as follows:   
f( r ) - f( r ) W( r,h ) dV . (2)
In this manner, the density of the particles could be expressed in terms of the sum of the kernel functions of N particles present within the radius of influence. Here, smoothing length for surface h/d=1.15 is assumed.   
ρ i = j=1 N m j W( r ij ,h ) (3)
However, this assumption does not yield a sufficiently accurate density estimation in the vicinity of the phase boundary. To avoid this, the multidimensional moving least-squares interplant with constraint condition (CLS) method is introduced.13)

2.2. Governing Equations for Incompressible Flow and Discretization

The governing equations for a weakly-compressible viscous flow are based on the relationship between the velocity of sound and the flow density under adiabatic conditions as well as the Navier–Stokes (N-S) equations:   

( Dp Dρ ) S = c 2 (4)
  
ρ Dv Dt =-p+η 2 v+ρg+ F s . (5)
On the right-hand side of Eq. (5), the first, second, third, and fourth terms denote the pressure gradient, viscous force, gravity, and interfacial force, respectively. A density function is introduced that considers the boundary particles to smoothen each term. The N-S equation can then be written by considering the simple relationship ρ=m/V as follows:   
D v i Dt =- 1 m i j=1 N ( p i V i 2 + p j V j 2 ) W ij + η i m i j=1 N ( V i 2 + V j 2 ) r ij | r ij | 2 v ij W ij +g - 2 δ ij σ i d p 2 m i ( j=1 N E( | r ij | ) ) -1 j=1 N ( | r ij |- d p ) ( | r ij |-2h ) . (6)
Using the correlation of Eq. (4), the purpose of producing the time derivative of pressure, Tait’s equation of state can generally be used.13)   
p i = c 2 ρ 0 γ { ( ρ i ρ 0 ) γ -1 } (7)

2.3. Non-Newtonian Fluid Model

Several non-Newtonian models have been applied to the SPH method to describe the rheological properties of fluids. The non-Newtonian flow can be modeled effectively by the Cross model.16,17) The advantages of the Cross model over the Bingham plastic fluid model is that the effective viscosity is a continuous variable, and numerical instability is avoided. Now, the viscosity of the Bingham plastic fluid, η, is described as follows:   

η 0 -η η- η = ( K γ ˙ ) m (8)
where η0 and η are viscosity at zero and infinite shear rates, respectively, and K and m are constant parameters. γ ˙ is the shear rate, and it is simplified to   
γ ˙ = 2{ ( v x x ) 2 + ( v y y ) 2 + ( v z z ) 2 }+  ( v z x + v x z ) 2 + ( v y x + v x y ) 2 + ( v y z + v z y ) 2   (9)
Generally, M=1 is taken in the Cross model, and then the following equation is obtained for expressing η:   
η= η 0 +K η γ ˙ 1+K γ ˙ . (10)
Although the Cross model requires these four rheological parameters, it must define these parameters in terms of the Bingham plastic parameters. Under the Bingham plastic fluid model, the fluid at rest is capable of resisting any shear stress below the yield stress. When the yield stress is exceeded, the flow structure changes and the material behaves like a Newtonian fluid driven by the excess between the shear stress and the yield stress. When the shear stress falls below the yield stress, the fluid structure changes again, and there is either plug flow or no flow at all. Accordingly, for numerical computation, the effective viscosity of a Bingham fluid in each particle can be represented as   
η i = η B + τ B γ ˙ , (11)
where τB and ηB are Bingham yield stress and viscosity after yield, respectively. In summary, the viscous behavior of the molten slag is evaluated by substituting Eq. (11) into Eq. (6). By comparing between Eqs. (3) and (4), when the condition of K γ ˙ ≫1, the following parameters become K=ηB/τB and η=ηB. For the remaining unknown parameter η0, Shao and Lo proposed η0=103η to maintain the accuracy of the numerical solution.17)

2.4. Calculation Condition

First, we derive the digital data of the packed coke structure simulating the upper part of the raceway. To obtain the packing structure in a box-type container with sides measuring 0.2 m, forty pieces of 15 representative coke samples were numerically generated; thus, 600 pieces were prepared, and the initial position and rotation angles of each coke sample were determined by using a pseudorandom number. Figure 1 shows one of the calculation results. The calculation procedure of coke packing is the same as reported previously.14,19) Next, we simulated the slag trickle flow in the coke-packed bed. The packed bed obtained was used as the boundary condition for a no-slip condition. It should be noted that the packed bed structure used was the same for all calculations and only the physical properties of the slag changed. Molten slag having a width of 0.20 m and height of 0.01 m was placed immediately above the packed bed. This arrangement was set to reduce computational costs while achieving a steady flow. The flow behavior of molten slag was analyzed by the SPH method incorporating the Cross model. We set the particle diameter, d=1.0×10−3 m, and the time step, Δt=1.0×10−6 s, for the calculation condition. The slag composition was assumed to be the typical 45SiO2–45CaO–10Al2O3 at 1873 K: density ρ=2600 kg·m3, surface tension σ=0.5 N·m−1, and viscosity coefficient ηB=0.1 Pa·s.14,19) The physical property dataset of a solid-liquid coexisting slag is limited, but in general, τB and ηB tend to increase as the solid phase ratio increases in the slag. In this study, we assumed that the yield value τB varied from 0 to 20 Pa, and the influence on the flow characteristics is evaluated.

Fig. 1.

Digitized packed structure of 600 cokes in 15 different shapes (i.e., 40 cokes of the same shape were packed).

3. Results and Discussion

Figure 2 shows a snapshot of the flow behavior of molten slag passing through the coke bed at each yield value condition. Similar to a Newtonian fluid (τB=0 Pa), a Bingham fluid (τB>0 Pa) flows down a specific pass in the packed bed as an icicle shape. The dripping velocities are almost constant, regardless of τB values. Furthermore, relatively large droplets are trapped in substantially the same area (hs site) in the packed bed, irrespective of the value of τB. We see from this that the packed structure greatly affects the hs site. However, in the case of τB>0 Pa, small droplets are captured at locations other than the hs site. Figure 3 shows the time change of molten slag volume distribution passing through the coke-packed bed in the vertical direction. From the beginning until t=0.4 s, the shape of the profile curve makes a similar trickle-flow pattern because it barely changes, regardless of the value of τB. On the other hand, when the yield value τB=20 Pa at t=2.0 s, the stagnant slag volume peaks at h=0.08 m. As will be described later, after the stretched icicle-like liquid disperses, the number of droplets increases. The hydrostatic pressure of individual droplets then decreases, making apparent the influence of the yield stress.

Fig. 2.

The liquid phase velocity distributions in a 3-dimensional view for each yield stress condition: (a) τB=0.0, (b) τB=0.2, (c) τB=2.0, and (d) τB=20.0 Pa. (Online version in color.)

Fig. 3.

Effect of yield stress on dripping profiles of molten slag: (a) τB=0, (b) τB=0.2, (c) τB=2.0, and (d) τB=20.0 Pa. (Online version in color.)

Figure 4 shows the distribution of holdup droplets under each τB value and the volume distribution in the height direction. From the increase of the peak of the profile curve with the increase of the τB value, it is clear that the stagnant droplet volume increases at the same holdup sites in the packed bed. When droplets pass through the narrow “bottleneck” region in the packed bed, the velocity decreases due to the viscous force.14) Consequently, the volume of stagnant droplets increases as the τB value increases because the viscosity force increases locally when a droplet passes through the inter-coke void (refer to Fig. 7 for details). Figure 5 shows the relationship between the total hs and the modified capillary number. The non-dimensional correlation was compared using the relationship between the static holdup of molten slag in the coke bed by Ohgusu et al.: hs=9.96 (C p m * ) -1.38 .20) This equation considers only the force balance between gravity and surface tension, and the modified capillary number C p m * = ρg (ϕε D p ) 2 | σcosθ | ( 1-ε ) 2 doesn’t depend on viscosity. In the case of Newtonian fluid, although the value of hs obtained in the simulation underestimates the non-dimensional equation, the two almost coincide. On the other hand, in the case of τB>0 Pa, hs shows a tendency to increase as τB increases, and in the case of τB=20 Pa, hs was 2.4 times larger than it is when τB=0 Pa. To clarify the hs site change between the Newtonian fluids and Bingham fluids, the following analysis was carried out. The coincidence of the Bingham fluid hold-up particle coordinates was tracked for the particles existing at the hs site for the Newtonian fluid (τB=0 Pa). Here, the control volumes with smoothing lengths around each SPH particle are used as a reference for the coincidence of particle coordinates. Figure 6 shows the volume fractions of each hold-up site normalized by the droplet volume at the hs site of the Newtonian fluid. In the case of τB=0.2 and 2.0 Pa, 74–79% of droplets stay at the hs site of Newtonian fluid, and approximately 40% stay at different sites. As a result, the volume exceeding the Newtonian fluid is held up. There is almost no change between τB=0.2 and 2.0 Pa. On the other hand, in the case of τB=20 Pa, 97% is the hs site and 176% is the non-hs site. Thus, Bingham fluid droplets coincide with the Newtonian fluid hs site with high probability and increase the hold-up occupying different spaces with increasing yield stress. The mechanism is shown as follows. Figure 7 shows the distribution of viscosity, ηi, in the case of τB=2.0 Pa. The change in yield stress acting on the droplet occurs owing to the magnitude of ηi (i.e. inversely proportional to γ ˙ ). When the velocity of the Bingham fluid becomes larger, ηi becomes asymptotic to ηB=0.1 Pa s; conversely when the velocity of the Bingham fluid is small, ηi increases according to Eq. (11). As shown in the enlarged view in Fig. 7, ηi changes oscillations at the bottom of the drop while dropping. Here, it is clear that in addition the gravity and surface tension applied to the bottleneck, the yield stress of the Bingham fluid prevents smooth flow.

Fig. 4.

Effect of yield stress on holdup profiles of molten slag: (a) τB=0, (b) τB=0.2, (c) τB=2.0, and (d) τB=20.0 Pa. (Online version in color.)

Fig. 5.

Comparison of the holdup derived from the estimated values from previous studies and the calculated results obtained in this study. (Online version in color.)

Fig. 6.

Volume fraction at each yield stress value based on Newtonian fluid holdup volume.

Fig. 7.

Liquid phase viscosity ηi distributions as functions of Bingham yield stress τB along with a shear rate for the τB=2.0 Pa condition. (Online version in color.)

Next, the holdup characteristics will be considered by focusing on the morphology of individual stagnant droplets and the local structure of the coke bed. Figure 8 shows the volume distribution of holdup droplets at each τB value. The proportion of small droplets with Vd<0.0005×10−6 m3 increases uniformly with τB>0 Pa. Droplets are captured in regions other than the regions of hs sites. Compared with the particles in the bulk liquid, the number of particles around the dispersed droplet surface decreases. Therefore, the shear rate decreases and the viscosity increases. In the SPH method, the interaction with surrounding particles is taken into consideration within the influence region, so the possibility that such a cut-off distance affects the calculation accuracy cannot be denied. Therefore, the difference in hs value between 0<τB<2 Pa in Fig. 5 is due to droplets with Vd<0.0005×10−6 m3, so it will not be discussed in this study. On the other hand, the volume per droplet Vd occupies most of the range 0.1×10−6<Vd<1.0×10−6 m3, but the number of large droplets increases only when τB=20 Pa. This enlargement of a single droplet will occur when τB is greater than a threshold value. To consider this mechanism, the following examination was conducted. In the case of a Newtonian fluid, there is a correlation between the orthogonal projection area of the coke in contact with the holdup droplets and the total droplet volume in contact with the coke.19) To investigate the droplet volume distributions and special characteristics, this correlation was examined for each τB value. Figure 9 shows the number of droplets and the volume sum distribution classified by the various coke projected areas Aproj. The number of droplets nd and the volume of droplets Vd are classified by the cokes in contact with all holdup droplets. At any τB value, although there is a deviation in the data, the number of droplets and the droplet volume tended to increase as Aproj increased. The proportional relationships between Aproj and both nd and Vd indicate that the droplets are uniformly dispersed and coke with a larger Aproj contributes to the formation of holdup sites. The broken red lines passing through the origin are approximate curves for the plots of Vd. When τB=20 Pa, the gradient of the straight line is at its maximum, and nd tends to increase. In this mechanism, one droplet stagnates at the bottleneck, which has a hydrodynamic radius that is less than the capillary length. Compared with the results in Fig. 6, when τB=20 Pa, large droplets with Vd>0.1×10−6 m3 are not trapped by only a specific size of coke but are dispersed uniformly. In this case, as is seen with powder holdup characteristics,21,22,23) droplets are trapped one after another starting from a few remaining droplets. Furthermore, since the actual Bingham fluid increases ηB as the solid phase ratio increases, the shear rate at bottleneck sites is decreased further, so the holdup amount increases due to the effect of ηB.

Fig. 8.

Histograms of droplet volume distributions in the state of static holdup. (Online version in color.)

Fig. 9.

Relationship between the projected area of coke and holdup droplet (a) τB=0, (b) τB=0.2, (c) τB=2.0, and (d) τB=20 Pa. (Online version in color.)

4. Conclusions

A simulation of the dripping of solid-phase suspended molten slag in a coke bed with poor wettability was performed using the SPH Bingham plastic fluid model to which the Cross model was applied, and the following knowledge was obtained.

When a case study was performed with the yield coefficient being varied with the viscosity coefficient of the slag as approximated by the Bingham fluid model, the holdup droplets were trapped at almost the same site, regardless of the yield value. In addition, as the yield coefficient increased, the enlargement of droplets trapped at the same sites was confirmed. Regardless of the yield value, a correlation was confirmed between the projection area of the coke and the holdup droplet. When the yield value exceeded the threshold, holdup droplet volume increased. This threshold is correlated with the decrease in shear rate at the bottleneck.

In summary, an increase in the solid phase ratio in the molten slag is predictive of an increase in the yield value with viscosity; thus, the amount of holdup droplets increases at specific holdup sites as a dispersed phase. This indicates that the increase in the yield value approaches the powder holdup mechanism that is closed one after another, starting from a specific holdup site.

Acknowledgments

Part of this research supported by the Iron and Steel Institute of Japan (ISIJ) Research Promotion Grant. A part of this research was supported by the Steel Foundation for Environmental Protection Technology (SEPT) of Japan.

Nomenclature

Symbols

A: area (m2),

Cp: capillary number (−),

c: sound speed (m·s−1),

D: coke diameter (m),

d: particle diameter (m),

E: interparticle potential (−),

f: scalar function (−),

g: gravity (m·s−2),

h: smoothing length of influence radius (m),

hs: static holdup (−),

m: mass (kg),

n: number (−),

N: number of particles in effective radius (−),

p: pressure (kg·m−1·s−2),

r: position vector (m),

t: time (s),

v: velocity (m·s−1),

V: volume, (m3),

W: kernel function, (m−3)

Greek letters

γ : adiabatic exponent (= 7.0) (−),

γ ˙ : shear rate (s−1),

δ: delta function (−),

ε: void fraction (−),

η: viscosity (kg·m−1·s−1),

θ: contact angle (rad),

ρ: density (kg·m−3),

σ: interfacial tension (kg·s−2),

τ: yield stress (kg·m−1·s−2),

φ: shericity (−)

Subscripts

0: reference (initial),

B: Bingham fluid,

d: droplet,

i: particle index,

j: neighboring particle index around i,

m: modified,

p: particle,

proj: projection,

S: isentropic

Appendix.

Observation on Free Surface Flow Characteristics of Solid Phase Suspension

In Bingham fluid simulations using the Cross model,16,17) there is little knowledge about the accuracy of free surface flow predictions. Therefore, a simple accuracy verification based on a comparison with room temperature model experiments was performed.

A glass conical container (funnel), having a height of 75 mm and an outflow hole with a diameter of 8 mm, was filled with a suspension and fixed with a rubber stopper. Simulating the dripping phenomenon at the packed bed bottleneck, the free surface shape of the suspension flowing out was recorded with a high-speed camera (FASTCAM SA3, Photron Co., Ltd., 500 fps). This experiment and the same-scale simulation results by the numerical model were compared. For the Bingham fluid suspension, silicone oil (Shin-Etsu Chemical Co., Ltd., density: 970 kg·m−3, viscosity: 1.0 Pa·s, surface tension: 2.1×10−2 N·m−1), and polyethylene beads (Sumitomo Seika Chemicals Co., Ltd., density: 920 kg·m−3, mean diameter: 162.5×10−6 m) were used. The volume fraction (solid-phase ratio) fs was changed from 0% to 60%, as reported by Sukenaga et al.24)

Figure A1 shows snapshots of the flow behavior of the suspension column when fs is changed. At low fs values, the suspension column quickly flows down. When fs>0.30, the liquid column flows down slowly, and it hardly flows for 4.0 s at fs=0.60. The suspension column diameter reaches a maximum value at the beginning of the flow and then adopts a flow behavior that repeatedly increases and decreases at fs>0.30. When fs increased, it is expected that the suspended particles exposed on the free surface increase and aggregate due to liquid bridging between the particles. Sukenaga et al. reported that the transition to Bingham flow characteristics occurs with yield stress at fs>0.30.24) Figure A2 shows a comparison between the experimental results and the simulation results. Figures A2(a)–A2(d) shows a comparison of the results for t=1.0 s when fs=0 and 0.45. A previously reported yield stress τB=2.3 Pa was used for fs=0.45.24) When fs=0, continuous flow was confirmed. On the other hand, when fs=0.45, the suspension column became discontinuous at times, and this tendency shows a good agreement with the experimental results. However, the suspension column became discontinuous at the upper part in the simulation results, but in the experimental results, the tendency was different.

Fig. A1.

Snapshots of dripping liquid columns in each solid fraction condition.

Fig. A2.

Snapshots of the free surface profiles of dripping liquid columns. (a)–(d) show comparisons of experimental and calculated results at t=1.0 s: (a) and (b) fs=0, (c) and (d) fs=0.45; (e)–(h) show comparisons of each surface tension condition: (e) σ=0, (f) σ=2.1×10−2, (g) σ=4.2×10−2, and (h) σ=8.4×10−2 N/m.

In the experiment, the force acting on the suspension surface changes because the beads are exposed on the surface of the liquid phase, and the surface area of the liquid phase decreases. Consequently, the local curvature changes on the free surface cannot be tracked at the minimum resolution of 1 mm. Note that the breaking position cannot be predicted when a calculation resolution sufficiently larger than the bead size (0.16 mm) is used for spatial averaging. Figures A2(e)–A2(h) shows the typical shape of the suspension column when the surface tension value is varied. As the surface tension increases, the position of the column discontinuity shifts in the direction of gravity. As shown in (h), when the value of the surface tension is increased approximately 4-fold, the shape becomes qualitatively close to (c). Thus, in the suspension phase, the apparent surface tension is different from that of the liquid phase.

References
 
© 2020 by The Iron and Steel Institute of Japan
feedback
Top