CFD Modeling of Multiphase Flow in an SKS Furnace

The emerging bottom blown copper smelting (SKS) technology has attracted growing interest since it came into production. To further reveal the agitation behavior inside the bath and optimize the variable parameters, CFD simulation was conducted on a scaled down SKS furnace model with diﬀerent tuyere arrangements. The Multi-Fluid VOF model was used for the ﬁrst time in SKS furnace simulation and the simulated results show good agreement with an experimental water model reported in the literature, in terms of plume shape and surface wave. It was found that a low velocity region would appear on the opposite side of the bubble plume and persisted for a long time. To enhance the agitation in the low velocity region and reduce the dead zone area, an arrangement with tuyeres installed at each side of the furnace was recommended. Results suggested that a smaller tuyere angle diﬀerence would help to strengthen the agitation in the system. However, further investigation indicated that the diﬀerence in tuyere angle between two rows of tuyeres should be limited within a certain range to balance the requirements of higher agitation eﬃciency and longer lining refractory lifespan.


I. INTRODUCTION
MODERN copper making technologies have undergone long-term development as the demand for copper materials keeps expanding. [1]Some copper smelting furnaces that are widely used in modern industry are continuously being optimized. [2]Among these technologies, bottom blown copper smelting technology, referred to as SKS (ShuiKouShan) technology, has attracted growing interest since the second decade of this century, for its good performance on production efficiency, energy saving, and environmental protection. [3,4]s a new and emerging technology, the approaches to optimize the furnace structure for SKS technology are mainly derived from industrial practice, and relevant research involving reaction mechanisms and flow phenomena is comparatively scarce.To further reveal and optimize the operation of SKS furnaces, additional modeling research is needed.7][8][9] The water model studies on the SKS furnace have mainly focused on the parameters related to mixing efficiency.[7] In order to approach the real situation, CFD simulations have also been carried out and the flow phenomena revealed, both digitally and visually.The VOF (Volume of Fraction) model, capable of tracking the interface of a multiphase flow, was applied in the simulation of SKS furnaces by Zhang et al. [8] In their research, a single tuyere model of an SKS furnace was built and the relation between mixing efficiency and bubble characteristics was investigated. [8]Another CFD study based on a scaled down SKS furnace model was conducted by Shao et al., [9] in which the Eulerian model was used and the mixing efficiency under different conditions was reported.The DPM (Discrete Phase Model) model has currently been used to describe the motion behavior of bubbles in gas stirred ladles that show similar flow patterns to SKS furnaces. [10]Detailed bubble behavior might be further revealed when DPM combined with the Eulerian or VOF model is applied to the SKS system in the future.The published CFD simulation work will help to optimize the numerical method and contribute to advancing the SKS furnace structure to some extent.However, the results may only be applicable in some special cases, as the simulated tuyere diameters and injection speed are not widely used, and the adopted numerical setup such as VOF and Eulerian models do not always perform well in the simulation of bubble behavior. [11]urrently, the reported CFD modeling work of the SKS furnace is not sufficient to meet the requirements of the copper industry, as the furnace structure parameters vary slightly in different copper production companies, [12] resulting in very different furnace parameters being used in various research studies.Therefore, for this newly adapted technology, continuous modeling work is still necessary for different types of furnace structure and tuyere arrangements.Due to the lack of a universal SKS furnace simulation that could be adapted in most cases, CFD modeling focusing on the general situation of flow phenomena in the SKS furnace has been conducted in this work, adopting the FangYuan furnace size reported in 2010. [13]The tuyere diameter of the FangYuan furnace is in the common reported industrial range of 48 to 75 mm [12] ; the simulation results of which would be referable for other types of SKS furnaces with similar geometrical parameters.
Based on the current situation, an optimization of the numerical model for SKS furnace simulation is required.To better simulate the jetting flow, the Multi-Fluid VOF model is used for the first time in this work.The Multi-Fluid VOF [14] model is provided for coupling the Eulerian and the VOF models, and has been proved to be successful in the simulation of some complex multiphase flow fields in recent years. [15,16]Similarly to the reported cases, the benefits of the Multi-Fluid VOF model are also noticed in this simulation.It is confirmed that both the continuity of bubble plume and the interface characteristics, which are well described by the Eulerian model and the VOF model, respectively, are balanced in the Multi-Fluid VOF model.For the results in terms of furnace operation, the velocity distribution and wall shear stress in a scaled down SKS CFD model are calculated for different tuyere arrangements, seeking a balance between better bath agitation and less impact on the lining refractory.Such results would help to understand the agitation efficiency inside the SKS furnace thoroughly and in some cases might be a reference for the further improvement of the furnace structure.
The primary goal of this CFD work is to verify the CFD model qualitatively with the help of published water model results of bath and gas plume flow patterns.The purpose of this study is also to analyze the macroscopic flow field inside the SKS furnace, and to provide a numerical method with the chosen turbulence and multifluid models, by comparing the CFD results with published water model flow images.Furthermore, the water model geometry is simulated using matte as the liquid phase.The model performance and calculation results pave the way towards the next step where a new water model with a different nozzle arrangement will be suggested and modeled.This will subsequently be simulated with a CFD model, which will then be used for predicting agitation behavior in an industrial scale SKS furnace.

II. MATHEMATICAL MODEL
The Multi-Fluid VOF model provides a framework to couple the VOF and Eulerian models.In general, for the Multi-Fluid VOF model, the flow equations are treated separately for each phase, and the symmetric and anisotropic drag laws which are not included in the single VOF model are also introduced.Additionally, along with the drag law, the Multi-Fluid VOF model also provides a series of interface sharpening schemes, of which the Geo-Reconstruct scheme is applied to this work. [14]The main mathematical expressions of the flow equations are given below.

A. Governing Equations
The mass and momentum conservation equations for phase q (including the gas phase and liquid phase) are given as Equations [1] and [2], respectively.
In Equations [1] through [3], a is the volume fraction,q is the density, v ! is the velocity, p is the pressure shared by both two phases, a q q q g ! is the gravity term, and f ! is an external body force which is specifically defined as the drag force in the current system.The subscripts g and l represent the gas and liquid phases, respectively.

B. Drag Force
Some of the interphase forces, such as the turbulence dispersion force, may influence the plume shape to some extent when the injection speed is in a comparatively low range, according to the report of Lou et al. [17] However, in the present simulation, a huge amount of drag force is induced due to a very high gas injection speed, which is far beyond the reported range.Further, as the detailed bubble behavior is not the main focus of this work, the interphase forces are ignored to improve the calculation efficiency, with the exception of the drag force.The drag force is introduced based on the symmetric model provided by the Multi-Fluid VOF model. [14]For the symmetric model, the density and viscosity are calculated from volume averaged properties: and the diameter of bubbles or droplets is defined as: METALLURGICAL AND MATERIALS TRANSACTIONS B VOLUME 52B, JUNE 2021-1773 In turn, the drag function could be given as: where the relative Reynold number Re is: and the drag coefficient C D is: The drag force is only affected by the bubble or droplet diameter, since the other variables are all derived from the industrial parameters.In this simulation, compared to the common situation in other copper making furnaces [18][19][20] and gas stirred ladles, [10,21] the gas injection speed is very high, transforming the bubbly flow into an almost jetting flow below the liquid level.The bubbles aggregate into a bubble group or bigger bubbles of irregular shape, based on experimental photographs from Shui et al. [6] This makes the bubble diameter vary a lot and difficult to evaluate.Therefore, as the bubbles tend to aggregate and coalesce, and the behavior of individual bubbles are not discussed in this work, the continuity of the bubble plume is preferentially considered.More specifically, the bubble plume is regarded as a bubble ribbon.This enables the gas phase to be set as a continuous phase and hence the droplet diameter is defined instead.Such an assumption should not affect the drag force, because in equation [6], d g ¼ d l , there is only one dispersed phase. [14]The diameter value is determined based on whether the plume shape (or ribbon shape) agrees with the experimental images, and finally the default value of 0.00001 m was adopted.

C. Turbulence Models
The standard k-e turbulence model with standard wall functions is adopted for the present simulation.The kinetic energy k and its rate of dissipation e yields are as follows:

½11
where G k represents the generation of turbulence energy due to the mean velocity gradients, and G b is the turbulence energy due to buoyancy.The mixture density q m and mixture velocity v !m are defined as follows: The empirical constants are C 1e ¼ 1:44, C 2e ¼ 1:92, r k ¼ 1, r k ¼ 1:3.The turbulence viscosity l t is given as Equation [14].

A. Physical Model
The prototype SKS furnace is cylindrical, with 9 tuyeres installed at the bottom.Enriched oxygen is injected via the tuyeres and impacts the copper matte inside, giving rise to violent agitation in the system.Generally, the key parameters affecting the mixing efficiency mainly include the tuyere angle, tuyere diameter, gas flow rate, and pool depth.The tuyere arrangement is the variable in this simulation and the others are fixed.The tuyere diameter is determined as 60 mm, which is a normal size in the reported range of 48 to 75 mm. [13]A small tuyere size as used by Shao et al. was not adopted in this research because it has been reported that the lifespan of small tuyeres with corresponding high gas injecting flow velocity would be significantly shortened. [22]A comparatively higher value of gas flow rate and pool depth was chosen to make the results closer to the expected situation in industry, as these parameters tend to be improved to achieve a higher production capacity.
In this work, a simulated SKS furnace model is created for the prototype in the ratio of 1:12 and the relative parameters are shown in Table I.The water model is the model from the experimental research of Shui et al. [6] The simulated SKS furnace is based on the water model of Shui et al., with same model size and similar nozzle arrangements.The volumetric gas flow rate was calculated based on the modified Froude number Fr¢, which is given as Equation [15]: where the subscript m represents the simulated SKS furnace model, subscript p represents the prototype, u is the volumetric gas flow rate, and L is the characteristic length.
1774-VOLUME 52B, JUNE 2021 METALLURGICAL AND MATERIALS TRANSACTIONS B In the simulated SKS furnace model, nine tuyeres are arranged in two rows with four at a larger angle and the other five at a smaller angle.Different tuyere angle combinations of 0 and 14 deg, 7 and 22 deg, À 7 and 14 deg, À 7 and 22 deg, were used in the simulation and analyzed in this work.The definition of the tuyere angle is shown in Figure 1.The tuyere with a positive angle value, for instance as 7 or 22 deg, is installed on the right side of the furnace, whereas the negative angle value is installed on the left side, as seen from the feed end where the tuyeres are located.
To improve the calculation efficiency, the following simplifications of the physical model were implemented: (1) The furnace structure was simplified as a simple cylinder.(2) There is no heat and mass transfer in the physical model.This would not have a considerable influence on the simulation accuracy as the key parameters of high-temperature matte and enriched oxygen remained.
(3) There is no slag phase present in the water model as a small amount of slag is not able to significantly affect the bubble plume.Also, the results in the current work could be regarded as a simulation towards an initial operating stage in an SKS furnace.

B. CFD Modeling
The mesh of the simulated SKS furnace is shown in Figure 2. Multizone meshing was adopted with hexahedral mesh as the majority.As the flow fields are only analyzed macroscopically in this simulation, which does not require detailed bubble behavior, the gas inlets of very small size were simplified and made square to reduce the cell number and improve the mesh quality, according to previously reported CFD studies on SKS furnaces and Peirce-Smith converters. [9,23]A preliminary mesh test was conducted to ensure that the macroscopic flow field was independent of the mesh size and inlet mesh simplification (Appendix).The number of cells with the current meshing setup was eventually set within the range of 230,000-280,000.For each mesh of cases with different tuyere arrangements, more than 99.72 pct of the elements are of a skewness below 0.5, which translates into good or excellent cell quality.The others, less than 0.28 pct, are of a skewness below 0.581, and are evaluated as fair cell quality. [24]he gas inlets are located at the bottom and set as the velocity inlet.The gas outlet is located at the top surface at a distance from the agitation zone and is set as the pressure outlet.To balance the performance of calculation speed and result accuracy, the Courant number is 0.25 as default, the time-step is set as 5 9 10 À4 seconds, and the convergence is marked when the number of residuals is below the level of 1 9 10 À3 .All control equations were calculated by the commercial software Ansys Fluent 2019 R3.
Table I.The Parameters of the Prototype, Water Model, Simulated Water Model, and Simulated SKS Furnace Prototype [5,12] Water Model [6] Verification

IV. RESULTS AND DISCUSSION
A. Numerical Model Performance

Simulated air-water system
The numerical model in this work has good modeling performance on geometrical similarity with the published water model results.The simulated bubble plume and surface wave, which would significantly affect the velocity distribution and wall shear stress analyzed in this research, show good geometric agreement between the literature water model (photographs) and the CFD model with an air-water system (blue images in Figure 3).Unlike the Eulerian model, which mostly describes a bubble plume with a very symmetric geometry, [9,17] the Multi-Fluid VOF model gives a more distorted plume shape, while maintaining the plume continuity to some extent.As seen in the blue images in Figure 3(a), the shape of the bubble plume is well reconstructed in the verification model of the air-water system.Both the plume continuity and bubble characteristics are to some extent balanced.In addition, the position of the bubble plume and the surface wave of the simulated model are also in good agreement with the experimental water model, as shown in the blue images in Figures 3(b) and (c).
It should be noticed that the simulated strong asymmetric and symmetric standing waves in Figures 3(b) and (c) do not appear only occasionally, but frequently, and different types of surface waves could only happen at a corresponding range of gas flow speed, as reported by Shui et al. [6] The simulated wave frequency and amplitude were investigated by further calculation of the air-water system shown in Figure 3(b), and the results are shown in Figure 4.In Figure 4, the time when the waves reached the highest amplitude were recorded and are shown using rhombic points.The interval between adjacent points could be regarded as the whole period of the wave, and a dimensionless frequency is calculated by where f is the natural frequency (Hz), R is the radius of the model (m), and g is the acceleration of gravity (9.8 m s À2 ).In the current test, the dimensionless frequency was 1.566, which is very close to the reported frequency of 1.495.The difference may come from errors in both the CFD and experimental tests.Additionally, the reported experimental frequency was measured after 20 minutes when the bath motion tended to be extremely stable.The wave amplitude increases slightly, indicating that the simulated bath motion was still in an unstable state, as shown in Figure 4.A calculation for 20 minutes is not acceptable for simulation; however, the intervals between each point are almost the same in Figure 4, indicating a strong regularity in the wave period.under the following conditions: pool depth 100 mm and 130 mm, gas flow rate 1.08 m 3 h À1 and 1.26 m 3 h À1 .
The wave amplitude was also measured and is shown in Figure 4.At 30 to 50 seconds, the wave amplitudes agree well with the estimated range according to the work from Shui et al. [6] The amplitude might further increase beyond the reported value after 50 seconds, based on the current variation trend.This could be a result of the interaction force between water and the backwall, because in the reported experiments, the tested wave was generated in the scaled down model but not in the slice model, as that used in the current simulation.
In general, the simulated wave frequency and amplitude strongly resemble the reported water model.The well reconstructed standing waves in the simulated air-water system indicate that a very similar macroscopic flow field has been established.This is further evidence to support the fact that the setup adopted in this simulation would contribute to making the calculation results closer to the real situation.

Simulated air-matte system
As shown in the orange images of Figure 3, in the air-matte system the flow regime of the bubble plume closely resembles that of the water model, with a very similar plume shape, plume position, and surface wave.Therefore, as proved by the good modeling performance in the water model simulation in Section IV-A-1, the numerical model in this work is also applicable for the simulation of the macroscopic SKS furnace system.Compared with the water model, the only difference for the air-matte interface is that it shows a thicker bubble plume.An intermediate cause is that the gas flow rate in the air-matte system is around twice as large as that of the experimental water model, based on the Froude similarity law.Concerning the physical property difference, the growth of plume size here could also be attributed to a different surface tension coefficient and liquid density.According to bubble rising research reported by Tripathi et al., [25] when bubbles rise in a quiescent liquid, the shape change of bubbles can be estimated by two important indicators including the Eo¨tvo¨s number (Eo q l gR 2 =r) and the Galilei number (Ga q l ffiffiffiffiffiffi gR p R=l l ), where q l is the liquid density, g is the gravitational acceleration, R is the initial radius of the spherical bubble, r is the surface tension coefficient, and l l is the liquid dynamic viscosity.In the current simulation, the Eo and Ga numbers of the water model and SKS furnace system are close, indicating that the bubbles in these two systems have very similar features.From this point of view, the physical property difference would not be the reason for a stronger bubble plume.However, on the other hand, such close Eo and Ga number values further prove that a numerical model showing good agreement with the corresponding water model also tends to perform well in the simulation of the SKS furnace system.B. Scaled Down SKS Furnace Simulation 1. Velocity distribution First, the relation between velocity distribution and bubble motion should be focused on for better understanding.In the simulated SKS furnace, the matte phase is dragged by the injected high-speed bubble plume, moving upwards and even out of the matte level, then dispersing around.The droplets of the dispersing flow at the upper and back areas move downwards, forming a circulation flow where the matte velocity is slightly higher than that in the surrounding area.The velocity field change in this simulation would be easier to analyze combined with the related bubble motion.In this work, the velocity distribution is displayed at selected sections of the simulated furnace model, which are shown in Figure 5.These sections are taken from the middle of the agitation zone, to involve the interactions between adjacent bubble plumes.Velocity distributions in and between the bubble plumes are analyzed separately in this simulation.First, the general situation of the matte velocity field is investigated.Figure 6 shows the matte velocity distribution in Section I (Figure 5) at 2 to 10 seconds, with a tuyere angle of 7 deg.The high velocity area where the velocity value is over 0.7 m s À1 overlaps to some extent with the location of the injected bubble plume.As shown in Figure 6, the high velocity area was not fixed at a certain position, but changed slightly at different times, which would be the result of a swinging bubble plume.More specifically, the bubble plume would not always move in a certain direction, but sometimes, for instance at 3 and 5 seconds, it moved in a different direction compared with the situations at other points of time.The change in the bubble plume path makes the velocity distribution vary from time to time.However, a common phenomenon is that mostly the low velocity region (beneath the liquid level) is located on the opposite side of the bubble plume, possibly at the center of the circulation and the near wall region that the circulation flow does not reach.In the low velocity region and dead zone (the absolute value of velocity is close to 0), the mixing process would be dramatically weakened, leading to an unexpected 'buckets effect' in the system.Therefore, the low velocity area is of great interest and will be discussed further.
The location of the low velocity region will be strongly associated with the tuyere angle.The velocity distribution in Sections I and II (Figure 5) with different tuyere arrangements including the angle combinations of 0 and 14 deg, ) is generated at the far side of the bubble plumes.The largest low velocity region is found when the combination of tuyere angle is 7 and 22 deg.This suggests that such a tuyere combination, which has been used in industry since 2010, [13] is not an efficient arrangement for mixing.
The occurrence of the large low velocity area could be attributed to the insufficient impact from the bubble plume when the tuyere is installed relatively far on the side, such as the situation of tuyere angle 22 deg.To lower the appearance frequency of such a low velocity region, arrangements with tuyeres installed on each side of the furnace are worth further consideration.A clearer low velocity distribution could be revealed by a shortened range of 0 to 0.1 m s À1 , which is shown on the right-hand side of Figure 7(b).It can be noted that when the tuyeres are installed at -7 deg, the low velocity area found in the sections of tuyere angles of 14 and 22 deg is comparatively reduced, with a smaller low velocity area, including the dead zone, remaining in the near wall region.Apparently, in these conditions, the bubble plume appearing at a negative angle strongly agitate the matte and was able to influence the velocity field at its adjacent bubble plume.This indicates that tuyeres arranged on each side of the furnace contribute to eliminating or at least reducing the area of low velocity.The tuyere combination arrangement on either side of the furnace of À 7 and 22 deg should be noted in Figure 7(b), as it is the only combination whose performance seems better than the case of 7 and 22 deg.This suggests that the arrangement on either side of the furnace is not always beneficial.To further reduce the low velocity region, the selection of tuyere angle must be considered.The low velocity area is significantly reduced with tuyere angle combinations of 0 and 14 deg, and À 7 and 14 deg, as shown in Figure 7(b).This might be the result of the agitation overlap of adjacent plumes that are relatively closer to each other when the tuyere angle difference is decreased.To prove this assumption, further evidence is discussed in the following sections.
The velocity field changes all the time as the bubble plume is not in a fixed position.Therefore, for further illustration, it is necessary to monitor the velocity variation, especially in the low velocity region and the dead zone.As shown in Figure 7(a), points 1 to 8 located in or near the dead zone were selected, for which velocity variations from 2 to 10 seconds were investigated separately.The monitoring results are shown in Figure 8, where it can be seen that the absolute value of the velocity at each of points 1 to 8 fluctuated and almost all were near or below 0.15 m s À1 .This indicates that the low velocity region in the current conditions would not be strongly agitated by a swinging bubble plume and would persist for a long time.However, even though it seems that the low velocity area is not easily eliminated, it is still possible to improve its performance by comparing the velocity fluctuation amplitude with different tuyere arrangements.
As seen in Figure 8, the velocity value at 10 seconds for all the points is below 0.05 m s À1 , but the fluctuation is different in each chart.It is obvious that the fluctuation amplitude for points 1, 2, 7, and 8 in charts (a) and (d) is larger than that for points 3 to 6 in charts (b) and (c).In other words, it could be considered that bubble plumes with tuyere angles of 0 and 14 deg and À 7 and 14 deg have a stronger impact on the nearby low velocity region.The comparatively quiet low velocity region appearing at tuyere angles of 7 and 22 deg, À 7 and 22 deg could be a result of an overlarge tuyere angle or tuyere angle difference.When tuyeres are installed at 7 and 22 deg, the high-speed flow is only active in the side region and is not able to take the far side liquid into the circulation.As regards the tuyere combination of À 7 and 22 deg, the major difference compared with the others is that it has the largest difference in tuyere angles.This means that the distance between the two rows of tuyeres is the largest and the interactions between the adjacent bubble plumes are minimized.Therefore, to strengthen the agitation in the (a) ( low velocity area, it is suggested to install tuyeres with a relatively small angle difference at each side of the furnace. The interactions between adjacent bubble plumes should not be ignored.The velocity distributions in Sections III through V (in Figure 5), in the middle of the tuyeres, are shown in Figure 9. Compared to the velocity distribution in the bubble plumes, the high velocity region in the middle sections shrinks to a small area, and tends to be much smaller when a larger angle difference is adopted.For instance, as shown in Figure 8(a), the combination of À 7 and 22 deg forms a larger area in dark blue, which means a very quiet region, and the comparatively higher velocity area is the smallest of the different tuyere angle arrangements.A clearer low velocity distribution is shown in Figure 9(b), in which the velocity range is further shortened.A larger low velocity region (black) is found with the angle combinations of 7 and 22 deg, and À 7 and 22 deg, which means that, with such tuyere arrangements, the interaction of adjacent plumes is weakened.This is in good agreement with the results given by the velocity variation, shown in Figure 8.Therefore, it can be concluded that such large tuyere angle differences could weaken the interactions of bubble plumes, not only impairing the agitation in the bubble plume sections, but also resulting in a large low velocity region in the middle of the adjacent bubble plumes.Comparing all the velocity distributions in Figure 9, the tuyere arrangement of 0 and 14 deg, À 7 and 14 deg show better performance in minimizing the low velocity region beneath the liquid level.This suggests that tuyeres installed with a relatively small angle difference would contribute to further improving bath agitation efficiency.

Wall shear stress
In most cases, the fiercer the agitation, the more impact it exerts on the side wall refractory.This leads to a conflict between the optimization of agitation efficiency and the expected refractory lifespan.In this work, the matte wall shear stress with different tuyere angle combinations is also taken into consideration.
Figure 10 shows a side view of the matte wall shear stress when different tuyere angles are used.As the tuyere structure is simplified in this simulation, the wall shear stress at the inlet region is not discussed.Generally, from the simulated case point of view, the strongest impact on the side wall is made by the matte phase, as the red and yellow region shown in Figure 10 is at the matte level.This means the largest force on the side wall is strongly associated with the surface wave, which may be generated frequently during furnace operation.According to the research of Shui et al., [6] standing waves may disappear in conditions caused by certain combinations of pool depth and gas flow rate.If these ranges are acceptable in industry, a standing wave on the matte surface might be avoidable.Apart from the surface wave, the combination of tuyere angles of 7 and 22 deg also leads to a fierce impact on the side wall shown in Figure 10.At such angles, the impact of the bubble plumes may overlap in some regions as the angle difference is comparatively small.This may be the reason why the wall shear stress rises when the tuyeres are set on only one side of the furnace.
The bottom view of the matte wall shear stress is illustrated in Figure 11, which shows a full picture of the wall shear stress.In Figure 11, the strongest impact on the bottom wall is found in the conditions of tuyere angle combinations of 0 and 14 deg, and 7 and 22 deg.These two angle combinations have a lower tuyere angle difference than the others, and could induce strong interactions between bubble plumes, as discussed in the chapters above.The strong interactions between adjacent bubble plumes expand the high velocity region, as shown in Figure 9, which could be a result of drag force convergence.Such intensification on the velocity field also affects the wall shear stress.It could be inferred that the overlapping of the impact on the bottom wall becomes the dominant factor that influences stress distribution.The angle combination of À 7 and 22 deg shows the best performance on reducing the wall shear stress, which indicates that the distance between two rows of the tuyeres should be wider to alleviate the wall shear stress overlap.However, this arrangement is not recommended because, from the aspect of velocity distribution, as discussed in Section IV-A, it would lead to the appearance of a comparatively quiet low velocity region in both the plume and middle sections.Therefore, to optimize the furnace operation, the tuyere angle difference should be within a feasible range to balance the performance of air agitation efficiency and refractory lifespan.

V. CONCLUSIONS
CFD simulation was carried out for a scaled down SKS furnace operating with different tuyere arrangements using ANSYS Fluent with a Multi-Fluid VOF model.Based on the results, the following conclusions can be drawn: 1.The numerical model adopted in this simulation results in very similar flow patterns to the water model reported in the literature.Therefore, it can be considered that it is feasible to apply the Multi-Fluid VOF model in the investigation of macroscopic flow field in the current work and further CFD simulations of the SKS furnace or other types of furnaces with corresponding geometries and tuyere arrangements.2. With a comparatively high gas flow rate in the range reported in the literature, the simulated bubble plume did not always move in an invariant path.As a conclusion, the bubble plume keeps swinging slightly in the system, resulting in a constantly changing velocity field in the system.3.In the CFD results in this study, some low velocity regions always existed on the far side of the bubble plume.As a conclusion, installing tuyeres on both sides of the furnace is recommended to reduce the size of that low velocity region.The results also suggest that two rows of tuyeres installed with a lower angle and smaller angle difference provide the low velocity region with stronger agitation.4. The liquid flow impact on the side wall mainly comes from the surface wave.On the bottom wall, a strong overlap of the wall shear stress seemed to appear when the distance of two tuyere rows was too narrow.It could be concluded that tuyeres with a bigger angle difference would help to alleviate the impact on the wall refractory. 5.A conclusion for industrial operations is that, although a larger angle difference would help to protect the wall refractory, it would not be beneficial for the agitation in the low velocity region.Consequently, the distance or the angle difference between the two rows of tuyeres should be selected in a limited range, to balance the requirements of agitation efficiency and lining refractory lifespan.6.With the current numerical model, the velocity distribution and wall shear stress at different pool depths, gas flow rates, and slag phase thicknesses merit further simulation.Such systematic investigations would help to understand the air agitation processes inside an SKS furnace better, helping to find more optimized operating parameters for industrial furnaces.

FUNDING
Open access funding provided by Aalto University.This work was supported by the China Scholarship Council, and School of Chemical engineering, Aalto university.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativec ommons.org/licenses/by/4.0/.APPENDIX A preliminary mesh test was conducted to confirm the reliability of the simplified mesh.Two single tuyere models were established with different mesh sizes and inlet geometries.The simplified inlet which has been used in this work was combined with a larger mesh size, of approximately 77,000 cells in total.A circular inlet was adopted with the refined mesh, which raised the total number of cells to around 200,000.A typical example of the mesh test results is shown in Figure A1, which indicates that the general flow field is not  A1.
In this work, according to the comparison results, the mesh was simplified due to the following points: 1.The simplified mesh did not significantly change the macroscopic features of the flow field, especially the area below liquid level compared with the refined mesh with circular inlet.2. To ensure the reliability of the simulation results, the focus points or conclusions in this simulation come from the common features between the tested cases with different mesh arrangements.The points of focus in this work include: (1) Swinging bubble plume.
(2) The effect of an adjacent bubble plume on each plume's velocity distribution.(3) Observation of the low velocity region for 1 to 10 seconds.(4) Impact on the wall mainly comes from a surface wave.(5) Overlap of the wall shear stress.
A schematic diagram of the tuyere angle definition.The positive and negative values of the tuyere angle correspond to tuyeres installed at different sides of the simulated SKS furnace.METALLURGICAL AND MATERIALS TRANSACTIONS B VOLUME 52B, JUNE 2021-1775

Fig. 2 -
Fig. 2-(a) The mesh of the simulated SKS furnace.(b) The mesh arrangement of the tuyere area.The cells in full black are the inlets where the tuyeres are located.

Fig. 3 -
Fig.3-Comparison of the water model images reprinted from Shui et al.[6] (left) and the verification model in current simulation (middle and right).The verification model images are the computed gas-liquid interface with a gas volume fraction of 0.5 for the air-water system (blue) and the air-matte system (orange).Common parameters: nozzle angle (a) = 7 deg, (b) = 7 deg, (c) = 25 deg, pool depth (a) = 120 mm, (b) = 120 mm, (c) = 100 mm.For the water model and computed air-water system, the gas flow rate is (a) = 0.96 m 3 h À1 , (b) = 1.20 m 3 h À1 , (c) = 1.56 m 3 h À1 .For the computed air-matte system, the gas flow rate calculated based on the Fr¢ number is (a) = 1.91 m 3 h À1 , (b) = 2.38 m 3 h À1 , (c) = 3.10 m 3 h À1 (Color figure online).

Fig. 4 -
Fig.4-Evaluation of the surface wave frequency and amplitude.Nozzle angle 7 deg, pool depth 120 mm, gas flow rate 1.20 m 3 h À1 .The estimated amplitude range is set based on the maximum and minimum experimental amplitude from Shui et al.,[6] under the following conditions: pool depth 100 mm and 130 mm, gas flow rate 1.08 m 3 h À1 and 1.26 m 3 h À1 .

Fig. 5 -
Fig. 5-Side view of the simulated SKS furnace.The iso-surface in orange is the air-matte interface with air volume fraction at 0.5 (Color figure online).
7 and 22 deg, À 7 and 14 deg, and À 7 and 22 deg, are shown in Figure 7(a).The bubble plumes are located in the high-speed region where the velocity is over around 0.7 m s À1 .When the tuyeres are installed at one side of the furnace (left part of Figure 7(a)), a large low velocity area (< 1 m s À1

Fig. 7 -
Fig. 7-Matte velocity distribution in Sections I and II of the simulated SKS furnace at 10 s.Gas flow rate = 23.4m 3 h À1 , pool depth = 120 mm.(a) Colored velocity distribution in the range of 0 to 1 m s À1 .(b) Black and white velocity distribution in the range of 0 to 0.1 m s À1 (Color figure online).

Fig. 9 -
Fig. 9-Matte velocity distribution in Sections III through V of the simulated SKS furnace with different tuyere arrangements at 10 s.Gas flow rate = 23.4m 3 h À1 , Pool depth = 120 mm.(a) Colored velocity distribution in the range of 0 to 1 m s À1 .(b) Black and white velocity distribution in the range of 0 to 0.1 m s À1 .

Fig. 10 -
Fig. 10-Side view of the matte wall shear stress of the simulated SKS furnace with different tuyere arrangements at 10 s.Gas flow rate = 23.4m 3 h À1 , Pool depth = 120 mm.

Fig. 11 -
Fig. 11-Bottom view of the matte wall shear stress at different tuyere arrangements at 10 s.Gas flow rate = 23.4m 3 h À1 , Pool depth = 120 mm.

Fig. A1 -
Fig. A1-Air-matte interface of air volume fraction at 0.5 with different mesh arrangement.Left: number of cells is around 200,000 with circular inlet.Right: number of cells is around 77,000 with simplified inlet.

Fig. A2 -
Fig. A2-Velocity distribution in the plume section from 2 to 5 s for different mesh arrangements.