Analysis of the Feeding Behavior in a Bottom-Blown Lead-Smelting Furnace

Abstract

Computational fluid dynamics (CFD) software was used to simulate the feeding behavior in a bottom-blown lead-smelting furnace. The results show that when the particle size is less than 30 μm, 20% of the particles are suspended in the gas phase and do not enter the melt pool for smelting, thus resulting in material loss. When the particle size exceeds 75 μm, the particles settle in the metal layer. When the particle size is 40–60 μm, the particles are distributed in the slag and metal phases, and the material is uniformly distributed in the molten pool; additionally, the average velocity of the particles exceeds 1.4 m/s, the average temperature exceeds 960 K, and the particles exhibit better behavior within this range, thus rendering it the optimal range of particle sizes for feeding.

1. Introduction

With increasingly stringent environmental requirements, conventional lead-smelting processes are being gradually eliminated and replaced by more environmentally friendly and efficient lead-smelting processes [1,2]. These new processes include the QSL one-step lead-smelting process, KIVCET lead-smelting method, Kaldo smelting method, and oxygen-enriched top-blown lead-smelting methods (such as the Isasmelt and Ausmelt methods) [3,4,5,6,7]. Based on the introduction and assimilation of QSL processes, China has independently developed the oxygen bottom-blown blast furnace and smelter reduction lead-smelting process (SKS method) [8,9,10]. For this process, the energy consumption of a crude-lead production unit is 50% lower than that of conventional lead-smelting processes, the sulfur recovery rate is over 96%, and the sulfur capture rate exceeds 99%, thus effectively solving the problems of low-altitude sulfur dioxide pollution and the dispersion of lead-containing dust inherent in conventional processes [11].

Owing to its high efficiency and environmental friendliness, the bottom-blown lead-smelting furnace is the core piece of equipment in oxygen bottom-blown lead smelting and the most mainstream technology applied in modern lead smelting [12]. In the bottom-blown lead-smelting process, oxygen is injected into the molten pool from the oxygen lance at the bottom of the furnace, where it reacts chemically with lead sulfide in the molten pool to generate lead oxide. This process is accompanied by the release of significant amounts of heat, which increases the temperature in the furnace, thus resulting in self-heated smelting. The stirring of the oxygen lance not only promotes heat and mass transfer, but also causes the melt to form a complex flow pattern. This stirring action causes the substances in the melt to be mixed comprehensively, thus improving the homogeneity and efficiency of the reaction. Furthermore, it enables the newly added concentrate and solvent to disperse and react promptly in the melt pool [13,14].

Oxygen-rich bottom-blown lead smelting is a complex multiphase flow process that involves heat transfer, mass transfer, fluid flow, and chemical reactions. Owing to the harsh reaction environment inside bottom-blown lead-smelting furnaces, some important physicochemical phenomena occurring inside the furnace cannot be observed. This problem is solved using numerical-simulation technology. In particular, the reaction in the smelting furnace can be effectively investigated by performing simulations using CFD software (FLUENT 15). The work is summarized in Table 1.

Previous studies have demonstrated that VOF multiphase flow and standard κ-ε turbulence models can well describe the fluid-flow behavior in a melt pool, whereas the DPM is an extremely effective method for describing particle behavior. The effects of changes in process parameters on the physical field in a smelting furnace have been investigated extensively via numerical simulations, thus facilitating the optimization of the process conditions and structural parameters. However, most studies pertaining to bottom-blown furnaces involve the use of copper- and zinc-refining systems. The lead-refining system has a high degree of reaction complexity, and its smelting behavior and flow patterns in the furnace differ from those of other bottom-blown furnace smelting systems. In industrial production, the particle size of the feed particles is one of the most important factors affecting the efficiency of lead refining. Therefore, this paper discusses the effects of the size of feed particles on the oxygen bottom-blown lead-refining process based on the simulation results of the particle-flow behavior for different particle sizes under the same operating conditions.

2. Model and Assumption

2.1. Geometric Models

The geometric model used in this simulation included the exhaust gas outlet and the interior space of a bottom-blowing lead-smelting furnace. A schematic of the furnace structure is shown in Figure 1. The diameters of the furnace body were d₁ = 2.2 m and d₂ = 2.36 m. The total length of the furnace was 11.05 m. The initial height of the slag layer was 0.7 m and the initial height of the metal layer was 0.5 m. Hexahedral grids were used to mesh the model of the bottom blowing furnace. The model was discretized using block grids and local refinement to adapt to the distribution characteristics of the computational data. The total number of grids was 700,000.

The entire bottom-blowing smelting process can be described as follows: oxygen-enriched air is sprayed into the melt pool by an oxygen gun and undergoes intense heat, mass, and transmission processes with the high-temperature melt in the melt pool. At the same time, due to the stirring effect of the oxygen-enriched air on the melt pool, the concentrate and flux added to the upper part of the melt are sucked up, thus completing the metallurgical reaction. The generated flue gas is discharged through the flue gas outlet at the upper part of the smelting furnace and enters the subsequent waste heat recovery, dust collection, and acid production sections. The smoke and dust in the flue gas enter the smelting furnace through the return system for smelting. The multiphase flow process between the oxygen enriched air and high-temperature melt in the melt pool is the key part of this process.

In the physical modeling process, the simplification and assumptions used are: the influence of the discharge of coarse lead and high-lead slag on the flow field inside the flame smelting furnace is ignored, only three oxygen guns are considered as the inlet, the flue is used as the outlet, and the others are treated as adiabatic walls. Inlet 1 of the spray gun is set as a velocity inlet boundary condition, with an inlet air velocity of 70 m/s.

The furnace charge (i.e., discharge particles) is composed of lead sulfide concentrate, flux, smoke, coal powder, and other lead-containing materials. The furnace charge is granulated and added to the bottom-blowing furnace.

In the bottom-blowing furnace oxidation-smelting stage, lead concentrate is granulated and added to the bottom-blowing furnace. Oxygen is blown into the melt pool through the bottom oxygen gun, and lead sulfide concentrate reacts with oxygen to produce lead oxide. Lead oxide then interacts with lead sulfide in the furnace charge to obtain a coarse lead and high lead slag.

The initial conditions and all the simplifying assumptions of their model are given in Appendix A.

2.2. Mathematical Models

2.2.1. Continuous Phase Model

The simulation was conducted using Fluent 14.5. The control equations for the momentum, mass, and energy transfer were selected according to the characteristics of the smelting process. Their general forms are given by Equation (1):

(1)$\nabla ·\left(\rho v\phi \right)=-\nabla ·\left({\Gamma }_{\phi }{\nabla }_{\phi }\right)+S_{\phi }+S_{a\phi }$

where $\rho$ represents the density, v denotes the kinematic viscosity coefficient, $\phi$ denotes the generic dependent variable, ${\Gamma }_{\phi }$ denotes the transport coefficient, $S_{\phi }$ denotes the source term of the continuity equation, and $S_{a\phi }$ denotes the source term resulting from the reaction between the gas and the particles. The specific forms of the variables in solving different equations are presented in Table 2.

In Table 2, ${\mu }_0$ denotes the molecular viscosity, ${\mu }_T$ denotes the turbulent viscosity, $p_r$ denotes the static pressure, $-q_r$ denotes the heat from radiation or chemical reactions, and $-w_s$ denotes the rate of substance generation during combustion or chemical reactions.

It should be noted that the equations for the VOF model are given in Appendix B.

2.2.2. Turbulence Model

Multiphase flow problems can be solved using turbulence models. Because of its wide applicability and economically reasonable accuracy, the $\kappa -\epsilon$ model has a wide range of applications in industrial flow field and heat exchange simulations [22]. The standard $\kappa -\epsilon$ model is a semiempirical equation based primarily on turbulent kinetic energy and diffusivity. The $\kappa$ equation is a precise equation. The $\epsilon$ equation is an equation derived from an empirical formula. These two equations are as follows:

(2)${\displaystyle \frac{\partial }{\partial t}}\left(\rho \kappa \right)+\nabla ·\left(\rho \kappa \overrightarrow{\nu }\right)=\nabla ·\left[\left({\displaystyle \frac{{\mu }_{\mathrm{T},\mathrm{m}}}{{\sigma }_{\kappa }}}\right)\nabla \kappa \right]+G_{\kappa ,m}+G_{b,m}-\rho \epsilon$

(3)${\displaystyle \frac{\partial }{\partial t}}\left(\rho \epsilon \right)+\nabla ·\left(\rho \epsilon \overrightarrow{\nu }\right)=\nabla ·\left[\left({\displaystyle \frac{{\mu }_{\mathrm{T},\mathrm{m}}}{{\sigma }_{\epsilon }}}\right)\nabla \epsilon \right]+{\displaystyle \frac{\epsilon }{\kappa }}\left(C_1{G}_{\kappa ,m}+C_1{C}_3{G}_{b,m}-C_2\rho \epsilon \right)$

where $\rho$ denotes the density, $\overrightarrow{\nu }$ denotes the velocity, ${\mu }_{T,m}$ denotes the turbulent viscosity, $G_{\kappa ,m}$ denotes the turbulent kinetic energy generated by the laminar velocity gradient, and $G_{b,m}$ denotes the turbulent kinetic energy generated by the buoyancy force. ${\sigma }_{\kappa }$ and ${\sigma }_{\epsilon }$ are the turbulent Prandtl numbers of $\kappa$ and $\epsilon$, respectively. The constants in the model are $C_1$ = 1.44, $C_2$ = 1.92, $C_3$ = 0, ${\sigma }_{\kappa }$ = 1.0, and ${\sigma }_{\epsilon }$ = 1.3 [23,24].

2.2.3. Discrete Phase Model

The simulation of concentrate particles was performed using a discrete phase model (DPM). Detailed information is given in Appendix C. It should be noted that the interfacial tension for the gas, slag and metallic phase was 0.02, 0.40, and 0.28 N·m⁻¹, respectively. The Lagrangian method was used to solve the particle trajectories. Taking the X-direction of the Cartesian coordinates as an example, the control equations are shown in Equations (4)–(6):

(4)${\displaystyle \frac{d{u}_p}{dt}}=F_D\left(u-u_p\right)+{\displaystyle \frac{g_x\left({\rho }_p-\rho \right)}{{\rho }_p}}+F_x$

(5)$R\mathrm{e}\equiv {\displaystyle \frac{\rho d_p\left|u-u_p\right|}{\mu }}{F}_D={\displaystyle \frac{18\mu }{{\rho }_p{d_p}^2}}\cdot {\displaystyle \frac{C_{D}R\mathrm{e}}{24}}$

(6)$R\mathrm{e}\equiv {\displaystyle \frac{\rho d_p\left|u-u_p\right|}{\mu }}$

where $u_p$ denotes the particle velocity, $u$ denotes the velocity of the continuous phase, $F_D\left(u-u_p\right)$ denotes the force acting on the particles per unit mass, $F_x$ denotes the additional acceleration term of the particles, ${\rho }_p$ denotes the density of the particles, $\rho$ denotes the density of the continuous phase, $\mu$ denotes the viscosity of the continuous phase, and $R\mathrm{e}$ denotes the relative Reynolds number.

2.3. Physical Parameter Testing

In this study, to make the fluid parameters used in the simulation consistent with those of the actual fluid, the RTW-10 melt property tester was used to determine the viscosity, density, and surface tension of the slag obtained from a company’s bottom-blowing lead-smelting furnace for industrial production [25]. The temperature inside the bottom-blowing lead-smelting furnace is generally between 1200 and 1300 °C. Therefore, in this study, the viscosity, density, and surface tension of the melt were determined at 1250 °C. The results are shown in Figure 2a, b, and c, respectively.

In Figure 2a, the viscosity of the slag decreases with an increase in the temperature. The slag viscosity decreased from 17.88 to 1.45 Pa·s when the temperature increased from 1200 to 1300 °C. According to the Arrhenius equation, the measured viscosity results from the experiment were fitted, and the relationship between viscosity and temperature was obtained. The correlation coefficient R² was 0.981. The relationship is given by Equation (7). Figure 2b,c show that both the density and surface tension of the slag decreased with an increase in the temperature. Meanwhile, there was an approximately linear relationship between the density and temperature as well as between the surface tension and temperature. Therefore, the density and temperature as well as the surface tension and temperature were fitted linearly from 1200 to 1300 °C. Their correlation coefficients (R²) were both 0.998. The linear relationship is given by Equations (8) and (9).

(7)$\mu =1.3937\times {10}^{-11}\cdot \exp \left({\displaystyle \frac{277840}{RT}}\right)\left(\mathrm{P}\mathrm{a}\cdot \mathrm{s}\right)\left(1200-1300 °\mathrm{C}\right)$

(8)$\rho =3.9022-0.0015T\left(\mathrm{g}/\mathrm{c}{\mathrm{m}}^3\right)\left(1200-1300 °\mathrm{C}\right)$

(9)$\sigma =1422.48-0.9352T\left(\mathrm{m}\mathrm{N}\cdot {\mathrm{m}}^{-1}\right)\left(1200-1300 °C\right)$

Finally, combined with the field operating conditions, the physical property parameters used in this study are presented in Table 3.

2.4. Boundary Conditions and Solution Methods

The parameters were set: an inlet flow rate of 70 m/s, a turbulence intensity of 2%, and a hydraulic diameter of 0.05 m. The exhaust gas outlet was set as the pressure outlet, with a negative outlet pressure of −200 Pa, a recirculating turbulence intensity of 2%, and a hydraulic diameter of 0.78 m. Those parameters were selected based on past experience, as well as other researchers’ work [26,27,28,29]. The wall boundary condition was set as the fixed wall. It was considered that the velocity of the fluid at the wall was zero [30,31], and the standard wall function was used in the near-wall region.

Equations (1)–(6) were solved using the commercial software Fluent 14.5, with the SIMPLEC algorithm for pressure–velocity coupling, the PRESTO format for the pressure-based segregated solver, and the second-order upwind format for the momentum equations [32,33]. With the minimum timestep of Δt = 1 × 10⁻⁴ s for the solution of the VOF model and the standard $\kappa -\epsilon$ model, the total calculation time for the bottom-blowing lead-smelting furnace model was 20 s. This value was chosen because after multiple iterations of calculation, the final residual results of each control equation were all less than 1 × 10⁻³, and the calculation results no longer changed with the iterations. Therefore, it can be considered that the calculation results converged, that is, within 20 s.

3. Simulation Results

3.1. Particle Velocity Distribution

Figure 3 presents the velocity distribution of particles with different sizes in the furnace, which can be used to determine the distribution of particles. As shown in Figure 3, after entering the furnace from above the furnace body, the particles fell into the molten pool under gravity. Affected by the upward buoyancy generated by the gas injected by the oxygen lances at the furnace bottom and the force of the flow field within the molten pool, the particles moved around in the furnace. The gas inlet consisted of three nozzles located at the bottom of the furnace, with a total effective area of 284.456 cm². The distance between the oxygen guns was 1.25 m, and the distance between the two discharge ports was 1.25 m. The diameter of the oxygen lance was 20 mm. The physical properties of the slag were inputted as one of the boundary conditions for the calculation. The explicit mode was chosen for calculation because it can show the final shape of the multiphase flow after stabilization, as well as the instantaneous distribution of the multiphase flow during the flow process.

Numerical simulations of the three-phase flow process inside the bottom-blown lead-smelting furnace were first conducted, and then the feeding process of the bottom-blown furnace were studied. To accelerate the convergence process and reduce the calculation difficulty, the calculated flow field results were used as the initial flow field state in the discrete phase model during initialization, which was imported into the calculation process of the discrete phase model. This gave the discrete phase flow factors and coupled the discrete phase and flow field for calculations.

Figure 3a,b indicate that when the particle size was <30 μm, the materials dispersed as soon as they entered the furnace from the feed port. Particles were dispersed in the furnace, and some particles overflowed from the exhaust gas outlet. This is possibly due to the fact that the particle size was too small, resulting in the gas drag force and thermophoretic force significantly exceeding the gravity of the particle. Figure 3c–e may show that with an increase in the particle size, particles fell to a certain height and then diffused after entering the furnace. The position where particles began to diffuse changed from the gas layer to the slag layer. This may be explained by the fact that the gas force was stronger when the particles were closer to the oxygen lances at the furnace bottom. When the gravity of the particles was balanced with the gas force, the particles began to diffuse. This may also be the reason why a larger particle size corresponded to a lower diffusion position. When the particle size was 40–60 μm, the falling materials broke through the gas–liquid interface under gravity, entered the melting pool for smelting reactions, and diffused to all parts of the melting pool under the impact of the jet flow. Figure 3f–h show that when the particle size was >75 μm, the materials entered the furnace and sank to the bottom. When the particle size was too large, the gravity of the material should have exceeded the gas force; thus, the particles fell directly into the bottom of the molten pool and then diffused.

The maximum and average velocities of particles of different sizes are presented in Figure 4. As shown, the variation in the particle velocity was small. As the particle size increased from 20 to 130 μm, the average velocity of the particles decreased from 1.53 to 0.73 m/s. The maximum velocity of the particles also decreased. According to Equations (4)–(6), a smaller particle size corresponds to a stronger gas drag force and higher average and maximum particle velocities.

3.2. Particle Temperature Field Distribution

Figure 5 shows the temperature distribution of particles with different sizes. The particle temperature was higher when the particles were closer to the furnace bottom. This is mainly because the lances at the furnace bottom injected oxygen-enriched air. In areas closer to these lances, the gas content in the melt increased, and the reaction between the particles and gas became more intense. Therefore, more heat was generated. This is also the reason why the low-temperature zone of the particles was near the feed port, while the high-temperature zone of the particles was above the lances at the furnace bottom. As the particle size increased, the high-temperature area at the furnace bottom expanded. This is because when the particle size was too large, the particles directly contacted the air at the furnace bottom and reacted exothermically without being diffused to other places. This not only caused the temperature at the furnace bottom to be too high, which significantly reduced the service life of the refractory bricks, but also caused the particles to accumulate, resulting in insufficient reaction and reducing the smelting efficiency.

It should be noted that due to the complex reactions in the melt pool, we mainly study the flow inside the furnace without considering chemical reactions. During initialization, the temperature distribution throughout the furnace is uniform, and the influence of temperature on the gas phase is ignored. In industrial production, many oxidation–reduction reactions occur in the melt inside the furnace, releasing a large amount of heat. For the convenience of calculations, in the simulation, the heat release of chemical reactions is equivalent to the heat generation rate, which means that the melt can release heat at a certain rate. For thermal field calculations, FLUENT users can customize the addition of the source term shown in Appendix D; the thermodynamic theory and the selected reaction kinetic model are given in Appendix E; and the solution method is referred to in Appendix F. The initial temperature of metal bulk is 300 K.

The distributions of the average and maximum temperatures of the particles are shown in Figure 6. With an increase in the particle size, the maximum temperature of the particles decreased from 1039 to 1020 K. This is because as the particle size decreased, the reaction rate between the particles and gas increased, leading to an increase in the heat release rate of the particles and a faster temperature rise. Additionally, as the particle size increased, the particle volume increased. The contact area between the particles and gas decreased, and the reaction rate decreased, reducing the heat release rate and temperature. Analysis of the average temperature of the particles revealed different temperature variation rules. When the particle size was 20–50 μm, the average temperature of the particles increased with the particle size, reaching a maximum of 970 K at 50 μm. When the particle size was 50–130 μm, the average temperature of the particles decreased with increases in the particle size.

It should be noted that due to current testing conditions, it is not possible to monitor the fluid inside the furnace in real time. We verified the reliability of the mathematical model by analyzing and comparing the results obtained from water model experiments and numerical calculations; the numerical simulation results are in good agreement with the literature [34,35]. Additionally, the multiphase flow model was not selected, and the reasons for this are given in Appendix G.

When the particle size was small, the reaction rate was high; thus, the maximum temperature of the particles was high. However, because the particle size was small, the particle distribution was dispersed. Consequently, some particles did not participate in the smelting reaction or had an incomplete smelting reaction. Therefore, the average temperature of the particles was low. As the particle size increased, the reaction rate of the particles decreased, along with the maximum temperature of the particles. Additionally, the particle distribution became more concentrated, and the degree of reaction for the particles increased. Thus, the average temperature of the particles increased. However, when the particle size was too large, the volume and mass of the particles made the diffusion difficult. This led to particle accumulation, making the particle reaction incomplete. Therefore, the average temperature of the particles decreased.

When the particle size was 40–60 μm, the maximum temperature and average temperature of the particles were at a higher level. Therefore, the material particle size should be controlled within this range during the feeding process.

Under the influence of the flow field inside the furnace, the velocity and temperature distribution of the particles were consistent with the velocity and temperature distribution of the flow field inside the furnace. The values in Figure 4 and Figure 6 are both global values.

3.3. Particle Escape Rate and Distribution in Three Phases

Table 4 presents the distribution and escape rates of particles of different size. The escape rate is the ratio of the number of particles escaping from the furnace to the total number of particles; this can reflect the loss rate of the material. When the escape rate is >0, particles inevitably flow out of the furnace with the exhaust gas. As indicated by Table 4, the particle escape rate was >0 when the particle size was 20 and 30 μm. In particular, when the particle size was 20 μm, the escape rate was 9%. The particle escape rate was zero when the particle size was 40, 50, 60, 75, 100, or 130 μm. When the particle size was <30 μm, the particles escaped, and the material loss was severe. Therefore, during industrial bottom-blowing lead smelting, the particle size should be >30 μm to reduce the material loss and relieve the pressure on dust collection. In this way, the full reaction of the materials can be guaranteed, and the smelting efficiency can be increased.

For particles sized 20–30 μm, over 20% were in the gas phase, about 50% in the slag phase, and less than 30% in the metal phase. As the particle size increased, their proportion in the gas phase decreased, while it increased in the metal phase, with the slag phase remaining stable. Particles smaller than 30 μm tended to diffuse in the gas phase due to their small size and mass. This caused material loss and increased dust collection pressure. For particles sized 40–60 μm, none were in the gas phase, over 35% were in the slag phase, and less than 65% in the metal phase. As the size increased, particles shifted more towards the metal phase. Particles sized 75–130 μm were found entirely in the metal phase, indicating that larger particles fell directly into it.

When the particle size was too small, severe material loss occurred. When the particle size was too large, the material sank to the bottom, and the reaction was uneven. Therefore, only the appropriate particle size can ensure efficient production. When the particle size was 40–60 μm, the particles were distributed in the slag and matte. Materials can uniformly diffuse in the molten pool and undergo smelting reactions.

The validation of the model with the literature can be found in Appendix H. It should also be noted that particles were detected in the considered phases following these steps: The DPM model was first set up, and then particle arrival at different layers (such as the gas phase layer, slag phase layer, metal layer) was selected in the Results Reports Sample Trajectory interface; that is, the different surfaces that the particle diffusion can reach was selected. The calculation started from the pre-calculation point and ended at the calculation completion point. The file recorded the particle information that reached different surfaces during the entire calculation process. Then, in the post-processing stage, the data was imported into Excel for statistics to obtain the number of particles in different regions.

In theory, the simulation calculation of high-temperature heterogeneous systems should take as much time as possible; in our work, the minimum time step for solving the model was 1 × 10⁻⁴ s, and the total running time of the bottom-blown lead-smelting furnace model was 20 s. Within this time, the calculation results were able to achieve convergence conditions very well; this result is similar to some other researchers’ findings [36,37].

4. Conclusions

The following conclusions are drawn:

1. With an increase in the particle size, the average velocity of the particles decreased from 1.53 to 0.73 m/s. The maximum velocity of the particles decreased sharply in the interval of 20–40 μm, tended to level off in the interval of 40–70 μm, and decreased again in the interval of 70–130 μm. In addition, as the particle size increased, the diffusion position of the particles moved downward, from the gas layer to the metal layer.

2. When the particle size was 20–50 μm, the average temperature increased with size, peaking at 970 K at 50 μm. For particles sized 50–130 μm, the average temperature decreased as the size increased, with the maximum temperature dropping from 1039 to 1020 K.

3. The particle escape rate was above zero for sizes under 30 μm, with a high of 8.57% at 20 μm, causing significant material loss. Particles under 30 μm were found in the gas, slag, and matte phases, with over 20% in the gas phase. Particles over 75 μm were only in the matte phase and unevenly distributed. Particles sized 40–60 μm were in the slag and matte phases, which promoted particle reactions.

4. Based on particle velocity, temperature, escape rate, and three-phase distribution, maintaining the particle size between 40–60 μm during feeding in bottom-blowing furnaces is optimal for the smelting reaction.

Author Contributions

Conceptualization, K.S. and X.J.; methodology, Y.Z.; software, K.S.; validation, W.G. and D.O.N.; formal analysis, K.S. and Y.Z.; investigation, W.G. and D.O.N.; resources, X.J.; data curation, K.S. and K.E.W.; writing—original draft preparation, K.S. and H.M.; writing—review and editing, H.M., D.O.N. and K.E.W.; visualization, K.S.; supervision, H.M.; project administration, H.M. and K.S.; funding acquisition, H.M. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Conflicts of Interest

Authors Kena Sun, Xiaowu Jie, Yonglu Zhang, Wei Gao and Hao Ma were employed by the BGRIMM Technology Group. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Bottom-blown lead furnace structural model, 1—oxygen lances, 2—feeds, 3—exhaust gas outlet.

Figure 2. Measurement of the slag’s physical property parameters: (a) viscosity, (b) density, (c) surface tension.

Figure 3. Velocity distribution of particles with different sizes.

Figure 4. Distributions of the average and maximum velocity of the particles.

Figure 5. Temperature distribution of particles with different sizes.

Figure 6. Distributions of the average and maximum temperatures of the particles.

Table 1

CFD simulations performed by various researchers.

Author	Focus Area	Ref.
Wang et al.	Used the VOF model to simulate the fluid flow in a copper bottom-blown smelting furnace, and improved the stirring effect in the molten bath by optimizing the oxygen injector arrangement and blowing parameters	[15]
Shao et al.	Developed a mathematical model to describe gas–liquid flow and mixing behavior in a bottom-blown oxygen–copper smelting furnace, and the model validation was carried out through a water model experiment	[16]
Wang et al.	Established a 1/5 scaled physical slice model of a bottom-blown oxygen–copper smelting furnace to investigate the flow field characteristics in the pool and in the freeboard above the pool	[17]
Wang et al.	Used The VOF multiphase flow model coupled with the standard κ-ε turbulence model to study the effects of the effect of radially and axially inclined mushroom heads on the flow distribution and the splash rate a bottom-blown copper smelting furnace	[18]
Cheng et al.	Used the discrete phase model (DPM) to analyze the effect of swirling gas inlet design on particle motion and decomposition in magnesia flash calciner (MFC)	[19]
Schmidt et al.	Predicted the size-dependent particle sedimentation and the risk areas for flue dust accretions by establishing a three-dimensional CFD model	[20]
Rajabi et al.	Studied the effects of sulfide concentrate particle size on pollutant emissions from a flash smelting furnace through numerical simulation	[21]

Table 2

General Equations for Turbulent Transfer Processes.

Equation	$\mathit{\phi }$	${\mathit{\Gamma }}_{\mathit{\phi }}$	${\mathit{S}}_{\mathit{\phi }}$
Continuity equation	1	0	0
Momentum equation	v	${\mu }_0+{\mu }_T$	$-\nabla p+\nabla \cdot \left(\mu \nabla v\right)$
Energy equation	h	${\displaystyle \frac{{\mu }_0}{p_r}}+{\displaystyle \frac{{\mu }_T}{{\sigma }_h}}$	$-q_r$
Species transport equation	Y1	${\displaystyle \frac{{\mu }_0}{p_r}}+{\displaystyle \frac{{\mu }_T}{{\sigma }_h}}$	$-w_s$

Table 3

Physical Parameters.

Type	Numerical Value
Metal density/g·cm−3	11.3406
Metal viscosity/cP	189.1
Slag density/g·cm−3	1.9934
Slag viscosity/cP	1464
Gas density/g·cm−3	1.02 × 10−3
Gas viscosity/cP	5.35 × 10−2

Table 4

Distributions and escape rates of particles of different sizes.

ParticleSize/μm	NumberTracked	Escaped	Remained	GasPhase	SlagPhase	MetalPhase	ParticleEscapeRate/%
20	175	15	160	42	84	35	8.57
30	175	2	173	35	82	56	1.14
40	175	0	175	0	91	84	0
50	175	0	175	0	70	105	0
60	175	0	175	0	56	119	0
75	175	0	175	0	0	175	0
100	175	0	175	0	0	175	0
130	175	0	175	0	0	175	0

Appendix A

To accurately simulate the physical and chemical processes inside the bottom-blown lead-smelting furnace, a reliable mathematical model was selected. Among the models considered, the VOF model was selected for the multiphase flow model, the standard model was selected for the turbulence model, and the discrete phase model was chosen for the influence of material particle behavior. Before analyzing the cutting behavior, this article first conducted a three-dimensional numerical simulation of the three-phase flow inside the bottom-blown lead-smelting furnace to obtain the flow field and thermal field distribution inside the bottom-blown furnace.

Appendix A1. Assumptions and Simplification

(1) Slag and metal phases are homogeneous, stable, and incompressible liquid phases. In order to simplify calculations, in engineering, the compressibility of the gas phase flow rate can be ignored in situations where it is not high, so it is treated as an incompressible fluid. Therefore, during the calculation process, it is assumed that the physical properties such as the density and viscosity of the gas phase and melt remain unchanged;

(2) Due to the complex reactions in the melt pool, this article mainly studies the flow process inside the furnace, without considering chemical reactions. During initialization, the temperature distribution throughout the furnace is uniform, and the influence of temperature on the gas phase is ignored.

(3) Due to the distance between the furnace top flue gas discharge outlet and the liquid level, and the fact that it only maintains a very weak negative pressure, it is believed that the main purpose of the flue gas discharge outlet is to recover the flue gas and prevent leakage, and it has almost no effect on the movement of the melt. Therefore, the influence of the flue gas discharge outlet is ignored, and it is considered that the upper boundary of the three-phase flow field is the liquid level.

Appendix A2. Solution Method

Appendix A2.1. Boundary Conditions

After grid generation, first import ANSYS CFX^® to set the inlet and outlet ports, and then read them into FLUENT. (1) Inlet conditions: The nozzle inlet is set as a velocity inlet boundary condition, with an inlet air velocity of 70 m/s, turbulence intensity of 2%, and a hydraulic diameter of 0.05 m. (2) Export conditions: The boundary conditions at the flue outlet are set as a pressure outlet, with a negative pressure of −200 Pa, a reflux turbulence intensity of 2%, and a hydraulic diameter of 0.78 m. (3) Wall boundary conditions: The wall boundary conditions adopt a fixed wall, assuming that the fluid velocity at the wall is zero, and a standard wall function is used in the near wall region.

Appendix A2.2. Solution Method

Non-steady state calculations were applied to the flow model inside the bottom-blowing furnace. In order to improve the accuracy of the calculation results, the pressure-based separation solver used a SIMPLEC format, and the momentum equation used a second-order upwind format. To ensure the easy convergence of the calculation process, both the turbulent kinetic energy and turbulent dispersion rate used a first-order upwind format. The minimum time step for model solving was Δt = 1 × 10⁻⁴ s, and the total running time of the bottom-blowing lead-smelting furnace model was calculated to be 20 s. The setting of FLUENT in the model solving process is shown in Table A1.

Table A1

FLUENT Settings.

Solver Type	Setting Value
Multiphase Flow Model	VOF Multiphase Flow Model
Turbulence Model	Standard κ-ε Turbulence Model, Standard Wall Function
Discrete Format	First-order Upwind Format
Solving Method	Standard SIMPLEC Algorithm
Momentum Equation	Second-order Upwind Format
Turbulent Kinetic Energy	First-order Upwind Format
Turbulent Dissipation Rate	First-order Upwind Format

Appendix B

In the VOF model, a set of momentum equations is applied between different fluid phases, and the volume fraction occupied by each fluid phase in each computational unit contained in the flow field is recorded during the calculation.

The VOF formula relies on the absence of interpenetration between two or more fluids (or phases). In the unit, if the volume fraction of the q-th phase fluid is expressed as ${\phi }_q$, the following three situations will occur:

When ${\phi }_q$ = 0, the q-th phase fluid does not exist within the computing unit.

When ${\phi }_q$ = 1, the q-th phase fluid is filled in the computing unit.

When 0 < ${\phi }_q$ < 1, there exists an interface between the q-th phase fluid and other phase fluids within the computing unit.

(1) Volume ratio equation:

For the q-th phase, the continuous equation for the volume ratio is as follows:(A1)$${\displaystyle \frac{1}{{\rho }_q}}\left[{\displaystyle \frac{\partial }{\partial t}}\left({\phi }_q{\rho }_q\right)+\nabla \cdot \left({\phi }_q{\rho }_q{\overrightarrow{\nu }}_q\right)=S_{{\phi }_q}+{\displaystyle \sum _{p=1}^n\left({\dot{m}}_{pq}-{\dot{m}}_{qp}\right)}\right]$$

In the formula: ${\rho }_q$ represents the density of the q-th phase; ${\phi }_q$ represents the volume fraction of the q-th phase; $v_q$ represents the velocity of the q-th phase; $S_{{\phi }_q}$ indicates the source phase, which is zero by default; ${\dot{m}}_{pq}$ represents the mass transport from fluid phase p to fluid phase q, and ${\dot{m}}_{qp}$ represents the mass transport from fluid phase q to fluid phase p. The volume fraction equation is not solved for the main phase, and the calculation of the main phase volume fraction is based on the following constraints:(A2)$${\displaystyle \sum _{q=1}^n{\phi }_q}=1$$

(2) Momentum equation

In the VOF model, the velocity field obtained by solving a single momentum equation within the region is shared among all phases. The momentum equation is as follows:(A3)$${\displaystyle \frac{\partial }{\partial \mathrm{t}}}\left(\rho \overrightarrow{\nu }\right)+\nabla ·\left(\rho \overrightarrow{\nu }\overrightarrow{\nu }\right)=-\nabla p+\nabla \cdot \left[\mu \left(\nabla \overrightarrow{\nu }+\nabla {\overrightarrow{\nu }}^{{\rm T}}\right)\right]+\rho \overrightarrow{\mathrm{g}}+\overrightarrow{\mathrm{F}}$$

In the formula: $\rho$ represents density, $\overrightarrow{\nu }$ represents velocity, $\mu$ represents viscosity, $p$ represents pressure, and $\overrightarrow{\mathrm{F}}$ represents the source term.

(3) Energy equation

The energy equation is shared among the phases, and its representation is as follows: (A4)$${\displaystyle \frac{\partial }{\partial \mathrm{t}}}\left(\rho E\right)+\nabla ·\left(\overrightarrow{\nu }\left(\rho E+p\right)\right)=\nabla \cdot \left(k_{eff}\nabla T\right)+S_h$$

In the formula, $\rho$ represents density, $\overrightarrow{\nu }$ represents velocity, $p$ represents pressure, $k_{eff}$ represents thermal conductivity, and $S_h$ includes chemical reaction heat and other user-defined volume heat source terms. Among these, $E$ represented as:(A5)$$E={\displaystyle \frac{{\displaystyle \sum _{q=1}^n{\alpha }_q{\rho }_q{E}_q}}{{\displaystyle \sum _{q=1}^n{\alpha }_q{\rho }_q}}}$$

In the formula, $E_q$ is influenced by the specific heat and temperature of the q-phase.

(4) Surface tension

The VOF model can incorporate the influence of adjacent interfacial surface tension. The additional surface tension in the calculation results is the source term of the momentum equation. Its formula is as follows: (A6)$$\mathrm{F}={\displaystyle \sum _{\mathrm{pairs} \mathrm{pq},\mathrm{p}＜\mathrm{q}}{\sigma }_{\mathrm{pq}}}{\displaystyle \frac{{\alpha }_{\mathrm{p}}{\rho }_{\mathrm{p}}{\kappa }_{\mathrm{q}}\nabla {\alpha }_{\mathrm{q}}+{\alpha }_{\mathrm{q}}{\rho }_{\mathrm{q}}{\kappa }_{\mathrm{p}}\nabla {\alpha }_{\mathrm{p}}}{\frac{1}{2}\left({\rho }_{\mathrm{p}}+{\rho }_{\mathrm{q}}\right)}}$$

In the formula: ${\sigma }_{\mathrm{pq}}$ represents the surface tension coefficient, ${\alpha }_{\mathrm{p}}$(${\alpha }_{\mathrm{q}}$) represents the gradient of the volume fraction of the p-th (q-th) phase, ${\kappa }_{\mathrm{p}}$ (${\kappa }_{\mathrm{q}}$) represents the surface curvature of the p-th (q-th) phase, and ${\rho }_{\mathrm{p}}$(${\rho }_{\mathrm{q}}$) represents the density of the p-th (q-th) phase.

Appendix C

For the study of the particle distribution behavior in the bottom-blown lead-smelting furnace, the following mathematical models were selected for simulation calculations based on relevant theories, as shown in Table A2. FLUENT settings are shown in Table A3.

Table A2

Mathematical Model Settings.

Type	Setting Value
Multiphase flow model	Not selected
Energy model	Selected
Turbulence model	Standard turbulence $\kappa -\epsilon$ model, standard wall function
Component equation	Finite reaction rate model
Discrete phase model	Selected

Table A3

FLUENT Settings.

Type	Setting Value
Solver type	Based on pressure
Speed expression	Absolute speed
Gravity (m/s2)	X = 0, Y = −9.81, Z = 0
Operating pressure	Atmospheric pressure

The attribute settings of the discrete phase model are shown in Table A4.

Table A4

Discrete phase model attribute settings.

Type	Setting Value
Exchange source items	Selected
Particle calculation steps	5000
Integral scale	5
Resistance parameter	Sphere
Low-order format	Explicit
Whether to couple energy and matter	Not selected
Force model	Particle drag force, thermophoretic force, Stokes force, Brownian force
Jet type	Surface
Particle type	Burning
Particle size distribution	Unified diameter
Particle size (m)	20 μm, 30 μm, 40 μm, 50 μm, 60 μm, 75 μm, 100 μm, 130 μm
Temperature (K)	300
Number of particles	175
Random orbit	Random discrete distribution trajectory model
Surface chemical reaction	customized

The tracking method for particles is steady-state tracking, with the Wall boundary set to Reflect and the Outlet boundary set to Escape. This means that particles will bounce back when encountering the Wall and flow out of the flow field when encountering the Outlet. The wall boundary condition adopts a fixed wall, assuming that the fluid velocity at the wall is zero, and a standard wall function is used in the near wall region.

Appendix D

For thermal field calculations, FLUENT users can customize the addition of the source term DEFINE_SOURCE(). The added source items are as follows:

Figure 6. Distributions of the average and maximum temperatures of the particles.

Appendix E

The equilibrium equations for particles are:(A7)$${\displaystyle \frac{d{u}_p}{dt}}=F_D\left(u-u_p\right)+{\displaystyle \frac{g_x\left({\rho }_p-\rho \right)}{{\rho }_p}}+F_x$$ (A8)$${\displaystyle \frac{d{u}_p}{dt}}=F_D\left(u-u_p\right)+{\displaystyle \frac{g_x\left({\rho }_p-\rho \right)}{{\rho }_p}}+F_x$$ (A9)$${\displaystyle \frac{d{w}_p}{dt}}=F_D\left(w-w_p\right)+{\displaystyle \frac{g_z\left({\rho }_p-\rho \right)}{{\rho }_p}}+F_z$$ (A10)$$F_D={\displaystyle \frac{18\mu }{{\rho }_p{d}_p^2}}\cdot {\displaystyle \frac{C_{D}Re}{24}}$$

In the formula: $u$($v$, $w$) represents the velocity in the x (y, z) direction of the fluid phase; $u_p$$u$($v_p$, $w_p$) represents the velocity in the x (y, z) direction of the particles; $F_x$, $F_y$, $F_z$ represent additional forces, including the Brownian force, thermal surge force, and Saffman lift; $\mu$ represents the fluid dynamic viscosity; $\rho$ represents the fluid density; ${\rho }_p$ represents the particle density; $d_p$ represents the particle diameter; and $Re$ represents the particle Reynolds number.

During the diffusion process, the material particles will undergo surface reactions with gas-phase components. The thermal equilibrium equations and mass change rate for particle surface reactions are:(A11)$$m_p{c}_p{\displaystyle \frac{d{T}_p}{dt}}=h{A}_p\left(T_{\infty }-T_p\right)-f_h{\displaystyle \frac{d{m}_p}{dt}}{H}_{reac}+A_p{\epsilon }_p\sigma \left({\theta }_R^4-T_p^4\right)$$ (A12)$${\displaystyle \frac{d{m}_p}{dt}}=-4\pi d_p{D}_{i,m}{\displaystyle \frac{Y_{ox}{T}_{\infty }{\rho }_g}{S_b\left(T_p+T_{\infty }\right)}}$$

In the formula: $c_p$ represents the particle specific heat, $T_p$ represents the particle temperature, $h$ represents the convective heat transfer coefficient, $T_{\infty }$ represents the continuous phase temperature, $H_{reac}$ represents the surface reaction heat release rate, ${\epsilon }_p$ represents the particle blackness (dimensionless), $\sigma$ represents the Stephen Boltzmann constant, and ${\theta }_R$ represents the radiation temperature. When calculating particle trajectories, there will be mutual influence between the discrete and continuous phases, and heat, mass, and momentum exchange will occur between the two phases. The calculation formula for the heat exchange is:(A13)$$Q=\left[{\displaystyle \frac{\overline{m_p}}{m_{p,0}}}{c}_p\Delta T_p+{\displaystyle \frac{\Delta m_p}{m_{p,0}}}\left(-h_{fg}+h_{pyrol}+{\displaystyle \underset{T_{ref}}{\overset{T_p}{\int }}{c}_{p,i}dT}\right)\right]\stackrel{•}{m_{p,0}}$$

In the formula: $m_p$ represents the average mass of the particles in the control body, $m_{p,0}$ represents the initial mass of the particles, $c_p$ represents the specific heat capacity of the particles, $\Delta T_p$ represents the temperature change of the particles in the control body, $\Delta m_p$ represents the mass change of the particles in the control body, $h_{fg}$ represents the latent heat of the volatilization analysis, $h_{pyrol}$ represents the heat required for pyrolysis during the volatilization analysis, $c_{p,i}$ represents the specific heat of the volatile matter precipitation, $T_p$ represents the temperature at which the particles leave the control body, $T_{ref}$ represents the reference temperature corresponding to enthalpy, and $\stackrel{•}{m_{p,0}}$ indicates the initial mass flow rate of the tracked particles.

The particle mass change value is:(A14)$$M={\displaystyle \frac{\Delta m_p}{m_{p,0}}}\stackrel{•}{m_{p,0}}$$

The change in particle momentum is:(A15)$$F={\displaystyle \sum \left(\frac{18\beta \mu C_{D}Re}{{\rho }_p{d}_p^{2}24}\left(u_p-u\right)+F_{other}\right)}\stackrel{•}{m_p}\Delta t$$

In the formula: $\mu$ represents the fluid viscosity, ${\rho }_p$ represents the particle density, $d_p$ represents the particle diameter, $Re$ represents the particle Reynolds number, $u_p$ represents the particle velocity, $u$ represents the fluid velocity, $C_D$ represents the drag coefficient, $\stackrel{•}{m_p}$ represents the particle mass flow rate, $\Delta t$ represents the time step, and $F_{other}$ represents other interfacial forces.

Appendix F

The pressure-based separation solver uses a SIMPLEC format, and the momentum equation uses a second-order upwind format. To ensure the easy convergence of the calculation process, both the turbulent kinetic energy and turbulent dispersion rate adopted a first-order upwind scheme. The minimum time step for model solving was Δt = 1 × 10⁻⁴ s, and the total running time of the bottom-blown lead-smelting furnace model was 20 s.

Appendix G

Before discussing the process, we first conducted a three-dimensional numerical simulation of the three-phase flow inside the bottom-blown lead-smelting furnace to obtain the flow field and thermal field distribution inside the furnace. Then, the influence of the cutting process was analyzed. In order to accelerate the convergence process and reduce the calculation difficulty, the calculated flow field results were used as the initial flow field state in the discrete-phase model during initialization, which gave the discrete phase flow factor and coupled the discrete phase and flow field calculations. Therefore, the multiphase flow model was not selected. At the same time, in order to facilitate the convergence of the model, corresponding adjustments were made to the iterative algorithm and the sub relaxation factor of the calculation. The iterative algorithm for the calculation is shown in Table A5, and the sub relaxation factor settings are shown in Table A6.

Table A5

Iterative Algorithm.

Solution Method	Set Value
Pressure–Velocity Coupling Method	SIMPLE
Gradient	Least Squares Element Method
Other Discretization	First Order Upwind Format

Table A6

Sub-relaxation factor settings.

Factors	Set Value
Pressure	0.3
Density	0.3
Body Forces	0.3
Momentum	0.3
Turbulent Kinetic Energy	0.3
Turbulent Dissipation Rate	0.3
Turbulent Viscosity	0.3
Species	0.3
Energy	0.3
Discrete Phase Sources	0.1

Appendix H

When comparing the experimental results of the water model in the literature with the numerical simulation results [34], as shown in Figure A1, it is found that the trends in the shape of the air masses shown in Figure A1a,b are basically the same, and that the slight differences are caused by the differences in the inlet flow rate and physical parameters between the melt and water.

Figure A1. Bubble morphology diagram of numerical calculation and water model experiment, (a) Numerical calculation of bubble morphology at the outlet of the spray gun; (b) Bubble morphology at the outlet of the spray gun in water model experiment, Adapted from Ref. [34].

Ref [35] proposed a quantitative analysis method for the ratio of bubble size to spray gun diameter, and the water model established in Ref [34] is similar to this study. Therefore, the results of the ratio of the bubble size to the spray gun diameter in its water model were compared with the simulation results of this study. The results are listed in Table A7.

Table A7

Comparison of simulation results with literature data.

Result	d1/d2
Simulation	6.5
Literature	6.7 [34]
Error/%	2.99

d1: Bubble diameter; d2: Lance diameter.

From the above quantitative analysis, the numerical simulation results are in good agreement with the literature. Therefore, the mathematical model of the bottom-blown lead-smelting furnace established by the author of this article has high reliability and can simulate the flow behavior inside the bottom-blown lead-smelting furnace.

References

1. Sohn, H.Y.; Olivas-Martinez, M. Lead and zinc production. Treatise on Process Metallurgy; Elsevier: Amsterdam, The Netherlands, 2024; pp. 605-624.

2. Gao, G.; Wang, C.Y.; Yin, F.; Cheng, Y.X.; Yang, W.J. Situation and technology progress of lead smelting in China. Adv. Mater. Res.; 2012; 581–582, pp. 904-911. [DOI: https://dx.doi.org/10.4028/www.scientific.net/AMR.581-582.904]

3. Li, Y.G.; Jing, X.D. A decade’s production practice of JCC Kivcet technique. J. Phys. Conf. Ser.; 2024; 2738, 012009.

4. Li, J.D.; Zhou, P.; Liao, Z.; Chai, L.Y.; Zhou, C.Q.; Zhang, L. CFD modelling and optimization of oxygen supply mode in KIVCET smelting process. Trans. Nonferrous Met. Soc. China; 2019; 29, pp. 1560-1568. [DOI: https://dx.doi.org/10.1016/S1003-6326(19)65063-6]

5. Zhao, H.L.; Liu, F.Q.; Xiao, Y.D.; Sohn, H.Y. Computational fluid dynamics simulation of gas–matte–slag three-phase flow in an ISASMELT furnace. Metall. Mater. Trans. B; 2021; 52, pp. 3767-3776. [DOI: https://dx.doi.org/10.1007/s11663-021-02290-8]

6. Zhao, H.L.; Lu, T.T.; Yin, P.; Mu, L.Z.; Liu, F.Q. An experimental and simulated study on gas-liquid flow and mixing behavior in an ISASMELT furnace. Metals; 2019; 9, 565. [DOI: https://dx.doi.org/10.3390/met9050565]

7. Rezende, J.; van Schalkwyk, R.F.; Reuter, M.A.; To Baben, M. A dynamic thermochemistry-based process model for lead smelting in the TSL process. J. Sustain. Metall.; 2021; 7, pp. 964-977. [DOI: https://dx.doi.org/10.1007/s40831-021-00387-7]

8. Queneau, P.E. The QSL reactor for lead and its prospects for Ni, Cu and Fe. JOM; 1989; 41, pp. 30-35. [DOI: https://dx.doi.org/10.1007/BF03220844]

9. Song, K.Z.; Jokilaakso, A. CFD modeling of multiphase flow in an SKS furnace with new tuyere arrangements. Metall. Mater. Trans. B; 2021; 53, pp. 253-272. [DOI: https://dx.doi.org/10.1007/s11663-021-02362-9]

10. Song, K.Z.; Jokilaakso, A. CFD modeling of multiphase flow in an SKS furnace: The effect of tuyere arrangements. Metall. Mater. Trans. B; 2021; 52, pp. 1772-1788. [DOI: https://dx.doi.org/10.1007/s11663-021-02145-2]

11. Li, D.B.; Zhang, Z.X. The application of new lead smelting technology of oxidizing in bottom blowing furnace and reduction smelting in blast furnace. Nonferrous Met. Extr. Metall.; 2003; 5, pp. 12-13.

12. Bai, L.; Xie, M.H. Pollution prevention and control measures for the bottom blowing furnace of a lead-smelting process, based on a mathematical model and simulation. J. Clean. Prod.; 2017; 159, pp. 432-445. [DOI: https://dx.doi.org/10.1016/j.jclepro.2017.05.057]

13. Jiang, B.C.; Guo, X.Y.; Wang, Q.M. Analysis of melt flow characteristics in large bottom-blowing furnace strengthened by oxygen lance jet at different positions. J. Sustain. Metall.; 2023; 9, pp. 1704-1715. [DOI: https://dx.doi.org/10.1007/s40831-023-00759-1]

14. Chen, L.; Hao, Z.D.; Yang, T.Z.; Liu, W.F.; Zhang, D.C.; Zhang, L.; Bin, S.; Bin, W.D. A comparison study of the oxygen-rich side blow furnace and the oxygen-rich bottom blow furnace for liquid high lead slag reduction. JOM; 2015; 67, pp. 1123-1129. [DOI: https://dx.doi.org/10.1007/s11837-015-1375-y]

15. Wang, W.; Mu, L.Z.; Zhao, H.L.; Cai, X.Y.; Liu, F.Q.; Sohn, H.Y. CFD Study on improvement of non-uniform stirring in a large bottom-blown copper smelting furnace. Min. Metall. Explor.; 2024; 41, pp. 1421-1435. [DOI: https://dx.doi.org/10.1007/s42461-024-00968-6]

16. Shao, P.; Jiang, L.P. Flow and mixing behavior in a new bottom blown copper smelting furnace. Int. J. Mol. Sci.; 2019; 20, 5757. [DOI: https://dx.doi.org/10.3390/ijms20225757] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31744044]

17. Wang, D.X.; Liu, Y.; Zhang, Z.M.; Zhang, T.A.; Li, X.L. PIV measurements on physical models of bottom blown oxygen copper smelting furnace. Can. Metall. Q.; 2017; 56, pp. 221-231. [DOI: https://dx.doi.org/10.1080/00084433.2017.1310362]

18. Wang, W.; Cai, X.Y.; Mu, L.Z.; Lu, T.T.; Lv, C.; Zhao, H.L.; Sohn, H.Y. CFD simulation of the effects of mushroom heads in a bottom-blown copper smelting furnace. Metall. Mater. Trans. B; 2023; 55, pp. 694-708. [DOI: https://dx.doi.org/10.1007/s11663-023-02984-1]

19. Cheng, D.; Xu, H.; Zhao, L.; Dong, H.; Zhang, Z. Effect of swirling gas inlet design on particle motion and decomposition in magnesite flash calciner. Chem. Eng. Res. Des.; 2024; 206, pp. 386-396. [DOI: https://dx.doi.org/10.1016/j.cherd.2024.04.057]

20. Schmidt, A.; Montenegro, V.; Reuter, M.; Charitos, A.; Stelter, M.; Richter, A. CFD study on the physical behavior of flue dust in an industrial-scale copper waste heat boiler. Metall. Mater. Trans. B; 2022; 53, pp. 537-547. [DOI: https://dx.doi.org/10.1007/s11663-021-02389-y]

21. Rajabi, N.; Ghodrat, M.; Moghiman, M. Numerical simulation of the effect of sulfide concentrate particle size on pollutant emission from flash smelting furnace. Int. J. Environ. Sci. Technol.; 2021; 18, pp. 2925-2936. [DOI: https://dx.doi.org/10.1007/s13762-020-03055-7]

22. Duraisamy, K.; Iaccarino, G.; Xiao, H. Turbulence modeling in the age of data. Annu. Rev. Fluid Mech.; 2019; 51, pp. 357-377. [DOI: https://dx.doi.org/10.1146/annurev-fluid-010518-040547]

23. Patel, V.C.; Rodi, W.; Scheuerer, G. Turbulence models for near-wall and low Reynolds number flows-a review. AIAA J.; 1985; 23, pp. 1308-1319. [DOI: https://dx.doi.org/10.2514/3.9086]

24. Nagano, Y.; Hishida, M. Improved form of the κ-ε model for wall turbulent shear flows. ASME J. Fluids Eng.; 1987; 109, pp. 156-160. [DOI: https://dx.doi.org/10.1115/1.3242636]

25. OuYang, K.; Dou, Z.; Zhang, T.; Liu, Y. Effect of ZnO/PbO and Fe/SiO₂ ratio on viscosity of lead smelting slags. J. Min. Metall. Sect. B Metall.; 2020; 56, pp. 27-33. [DOI: https://dx.doi.org/10.2298/JMMB190622052O]

26. Jensen, B.B.B. Numerical study of influence of inlet turbulence parameters on turbulence intensity in the flow domain: Incompressible flow in pipe system. Proc. Inst. Mech. Eng. Part E J. Process Mech. Eng.; 2007; 221, pp. 177-186. [DOI: https://dx.doi.org/10.1243/09544089JPME124]

27. Zamora, B.; Kaiser, A.S.; Viedma, A. On the effects of Rayleigh number and inlet turbulence intensity upon the buoyancy-induced mass flow rate in sloping and convergent channels. Int. J. Heat Mass Transf.; 2008; 51, pp. 4985-5000. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2008.02.040]

28. Minkowycz, W.J.; Abraham, J.P.; Sparrow, E.M. Numerical simulation of laminar breakdown and subsequent intermittent and turbulent flow in parallel-plate channels: Effects of inlet velocity profile and turbulence intensity. Int. J. Heat Mass Transf.; 2009; 52, pp. 4040-4046. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2009.03.041]

29. Chattopadhyay, H.; Murmu, S.C. Effect of inlet turbulence intensity on transport phenomena over bluff bodies. Int. J. Fluid Mech. Res.; 2020; 47, pp. 485-499. [DOI: https://dx.doi.org/10.1615/InterJFluidMechRes.2020034333]

30. Su, Y.; Davidson, J.H.; Kulacki, F.A. Numerical study on mixed convection from a constant wall temperature circular cylinder in zero-mean velocity oscillating cooling flows. Int. J. Heat Fluid Flow; 2013; 44, pp. 95-107. [DOI: https://dx.doi.org/10.1016/j.ijheatfluidflow.2013.04.015]

31. Son, J.S.; Hanratty, T.J. Velocity gradients at the wall for flow around a cylinder at Reynolds numbers from 5 × 10³ to 10⁵. J. Fluid Mech.; 1969; 35, pp. 353-368. [DOI: https://dx.doi.org/10.1017/S0022112069001157]

32. Kaya, F.; Karagoz, I.J.C.E. Performance analysis of numerical schemes in highly swirling turbulent flows in cyclones. Curr. Sci.; 2008; 94, pp. 1273-1278.

33. Bal, M.; Kayansalçik, G.; Ertunç, Ö.; Böke, Y.E. Benchmark study of 2D and 3D VOF simulations of a simplex nozzle using a hybrid RANS-LES approach. Fuel; 2022; 319, 123695. [DOI: https://dx.doi.org/10.1016/j.fuel.2022.123695]

34. Yan, H.J.; Xia, T.; Liu, L.; He, Q.A.; He, Z.J.; Li, Q.Q. Numerical simulation and structural optimization of gas-liquid two-phase flow in reduction furnace of lead-rich slag. Chin. J. Nonferrous Met.; 2014; 24, pp. 2642-2651.

35. Yan, H.J.; Liu, F.K.; Zhang, Z.Y.; Gao, Q.; Liu, L.; Cui, Z.X.; Shen, D.B. Influence of lance arrangement on bottom-blowing bath smelting process. Chin. J. Nonferrous Met.; 2012; 22, pp. 2393-2400.

36. Jiang, B.C.; Guo, X.Y.; Wang, Q.M. Two-dimensional analysis of melt pneumatic stirring in large bottom-blowing furnace based on CFD. Nonferrous Met.; 2023; 6, pp. 39-50. (In Chinese)

37. Hu, H.; Yang, L.; Guo, Y.; Chen, F.; Wang, S.; Zheng, F.; Li, B. Numerical simulation of bottom-blowing stirring in different smelting stages of electric arc furnace steelmaking. Metals; 2021; 11, 799. [DOI: https://dx.doi.org/10.3390/met11050799]