Introduction

The production of iron and steel is integral to modern society, but is also a significant contributor to the global production of \({\hbox {CO}}_{2}.\) The steelmaking business accounts for 7 to 9 pct of the direct annual emissions from use of fossil fuel, a majority of which come from the blast furnace.[1] This is not sustainable and must be addressed if steel is to have a future in society. The emissions from the steelmaking process must be reduced without compromising the quality of the steel or the production capacity. Many novel processes utilizing natural gas or hydrogen for reduction have been proposed as alternatives to the blast furnace process and are currently being developed, but their indistrial use is still limited.[1, 2]

The reduction of iron oxide to liquid iron requires energy and a reduction agent to remove the chemically bound oxygen from the iron oxide. Energy is required to heat the material to the reaction temperature as well as to supply the heat consumed by the endothermic reactions associated with reduction. Traditionally, a significant part of the energy is supplied as heat from the combustion of coke, as done in the blast furnace. To reduce the emissions of the blast furnace, this heat can instead be supplied by electric energy from electric arcs or heated gases. If the electricity used for heating is produced in a renewable way the total emissions of the process would be reduced. This is the inspiration for the development of the IronArc process.

By utilizing electric arcs to superheat gas in a plasma generator (PG), energy can be supplied to metallurgical processes without using combustion of coke or hydrocarbons. If this is combined with the addition of hydrocarbons from liquid natural gas (LNG) or liquid petroleum gas (LPG), the reduction of iron ore can be done with significantly lower \({\hbox {CO}}_{2}\) emissions. Previous studies of this process by Bölke et al. by use of computational fluid dynamics simulations (CFD) and water-modeling found promising results.[3]

Currently, there exists a functioning pilot size plant of the IronArc process. However, to further prove the functionality of the process, a larger demonstration-scale reactor is planned. The demonstration-scale plant will feature a two-reactor approach to facilitate a continuous process. The first reactor is a melting reactor where iron ore will be added and heated using several PGs. In this reactor, hydrocarbons are added to reduce the \({\hbox {Fe}}_{2}{\hbox {O}}_{3}\) and \({\hbox {Fe}}_{3}{\hbox {O}}_{4}\) to FeO. Some gangue elements are also expected, which in combination with the remaining FeO will form a slag consisting of approximately 90 pct FeO, 5 pct \({\text {SiO}}_{2},\) and 5 pct CaO. This slag is then transferred to the second reactor where the final reduction takes place with the addition of carbon.

The transfer to the second chamber is done to facilitate the reduction of FeO to Fe, as it requires more reducing conditions. By collecting the CO in the off-gas in the second reactor and recycling it to the first reactor the efficiency of the process is optimized. If all injected carbon is used for reduction and all CO is recycled it is possible to reach a 88 pct use of CO for reduction. This would result in carbon use as low as 0.172 kg C/kg Fe as can be seen in the Rist diagram as shown by J. O. Edström et al. based on the work of Rist and Maysson.[4, 5]

For the IronArc process to be continuous, the reactor walls must be able to withstand contact with the slag during a significant amount of time. This is especially important in the transfer tube between the two reactors where the slag is always high in FeO and in constant motion. FeO is infamous for wearing down most refractory materials very quickly when present in steelmaking slags. In conventional steelmaking slags, the amount of FeO is seldom higher than 30 wt pct.[6,7,8] However, in the IronArc process a slag containing over 90 wt pct FeO is expected, which is a completely different case.

By studying the thermodynamic system of the slag combined with refractory materials it is possible to predict that many systems form a liquid phase at the melting temperature of the slag. Such a combination of refractory and slag would result in the loss of refractory material during the process. The proposed solution to this is to form a frozen layer of slag on the reactor wall to protect it from the molten slag. A freeze-lining is formed by applying an intense water-cooling on the outside of a material with a high heat conductivity. If the interface temperature between the slag and the reactor wall is sufficiently low the contact layer will solidify and protect the reactor wall from corrosive and mechanical wear.[9, 10] A complete freeze-lining will prevent any type of wear on the refractory and enable the demonstration size reactor to be a continuous process. In the pilot plant of the IronArc reactor a freeze-lining is used to protect the reactor wall. Additionally, a thin layer of \({\hbox {Al}}_{2}{\hbox {O}}_{3}\) refractory is used in parts of the reactor to provide extra protection.

To ensure a complete freeze-lining the cooling and the geometry of the reactor must be tailored accordingly. Using the one-dimensional heat equation it is simple to determine the temperature profile and heat loss in a system with a set wall- and freeze-lining thickness, given stationary conditions. All components of the wall will experience an equal heat flux, which can be calculated by their combined thermal resistance. However, in practice, the freeze-lining will also be affected by the convection and movements in the melt. The changing viscosity of the melt near liquidus temperature and the change in thermal conductivity between solid and molten slag will also affect the solidification behavior to some degree.[11, 12] For this reason, numerical simulations are used to determine the solidification behavior with variations of these parameters.

Previous work on numerical simulations of solidifying slags have been done by Guevara and Irons on natural convection in a fayalite slag.[13] Their study concluded that the temperature–viscosity relationship is essential for the accuracy of the simulation. Studies on the thermal and viscous boundary layers in solidification applications have been studied by Voller et al.[14,15,16] These studies use the enthalpy–porosity formulation for the flow in the solidifying region. This model is also used by ANSYS Fluent in simulations of the mushy-zone layer. In the enthalpy–porosity model, the velocity of the flow is modified by inserting a sink term to replicate the slower movement in the solidifying region. By assigning every cell in the solidifying region a porosity value between 0 and 1 depending on the liquid temperature, the velocity of the flow is lowered.[17] The sink term reaches its maximum value at the solidus temperature, effectively causing the velocity to approach zero as the liquid approaches solidification temperature. The sink term also includes the mushy-zone parameter which determines the magnitude of the slowing effect and can greatly affect the results of the simulation.

Recently studies have been done to determine appropriate values of the mushy-zone parameter to produce realistic solidification behaviors when using the enthalpy–porosity model. The experimental results from studies on solidification of the low-temperature phase changing material (PCM) lauric acid by Shokouhman and Kamkari have been used by Fadl et al. and Kherabadi et al. as validation for numerical models.[18, 19] Fadl et al. studied the accuracy of solidification with mushy zone values in the range of \(10^{5}\) to \(10^{7}\) and found that increasing mushy zone parameter decreases the fluid velocity. This, in turn decreases the convection and rate of heat transfer. They concluded that the importance of the mushy zone parameter decreases in regions with high conductive heat transfer, and increases in importance in regions where convection is the dominating heat transfer mechanism. Depending on the type of flow and case setup the mushy zone parameter which fits best to the experimental validation changes. This indicates that there may not be a set value of the mushy zone parameter best fitted to all cases, but that it is case dependent.[20] Kheirabadi et al. found that high values of the mushy zone parameter give slower melting and that too small values give unphysical melt front development. They also concluded that the melting interval for the liquid (\(\Updelta T\)) and the mushy zone parameter were not independent, further indicating that the appropriate mushy zone parameter value is case dependant.[21]

In a paper by Yang et al. the mushy zone parameter is evaluated and a relationship between the mushy zone parameter (\(A_{\text {mush}}\)) and the dynamic viscosity (\(\mu \)) and the length-scale of the microstructure (d) in the mushy zone is proposed.[22] Since the thickness and behavior of the porous region (mushy zone) is dependent on the speed and behavior of the dendrite growth, such a relationship is useful to predict a reasonable mushy zone parameter for the simulation. However, no data on the dendrite formation in the current slag system is available for comparison.

$$ A_{\text {mush}}=\frac{180\mu }{d^{2}}. $$
(1)

Such a relation makes the mushy zone parameter both material dependent and case dependent. Yang et al. also evaluated the impact of the temperature interval between solidus and liquidus used in the enthalpy–porosity approach and found that an isothermal assumption may produce wrong predictions and that a correct temperature interval was a requirement for proper simulation of the solidification front.[22]

The effect of flow rate on the solidification behavior was studied by Crivits et al. using both simulations and physical experiments in the \({\text {CaCl}}{-}{\text {H}}_{2}{\text {O}}\) system. They found that the thickness of the freeze-lining, the interface temperature, and the phases which formed during solidification were sensitive to the parameters that determine the chemical driving force, but not very sensitive to the bulk fluid velocity.[23] In earlier work by Crivits et al. they found that lower viscosity results in higher interface temperature at the freeze-lining. It was also confirmed further that the thinnest areas of the freeze-lining occur in the areas of high turbulence, as previously suggested by Fraser et al.[24, 25]

In the present study, a numerical model was established to aid in the design of a slag transfer system which can form a complete freeze-lining. Two different geometries of slag transfer tubes were tested, a standard runner design called “Slagrunner” (SR), and a slower-moving design called “Stagnated Bath” (SB) which are shown in Figure 1. In these slag transfer tubes a parameter study was performed to determine what parameters of the flow affected the formation of a freeze-lining the most. The thermal and viscous properties of the slag at the given composition are not known as they have not been experimentally determined to the authors knowledge. For this reason a variation of said parameters is of interest to determine the range in which a freeze-lining is applicable. Based on the experimental studies in the \({\text {CaCl}}{-}{\text {H}}_{2}{\text {O}}\) system by Crivits et al. a validation study was performed. This evaluated how accurately the solidification module in ANSYS Fluent replicated solidification with the settings used in the SR and SB parameter studies. The contribution of the current work to the sciences is a focus on mushy zone parameter, viscosity, heat conductivity and varying geometries with higher flow rates than what has previously been evaluated, to the authors knowledge.

Fig. 1
figure 1

The slagrunner (SR) design (a) and the stagnated bath (SB) design (b)

Full size image

Simulation Theory

In finite volume modeling a set of conservation equations are used to calculate the behavior of conduction, convection, diffusion, source terms, and other mechanisms. The conservations equations are solved in each computational cell and then passed on to neighboring cells. The general equation describing this is presented in Eq. [2], where \(\rho \)\(({\text {kg}} \,{\text {m}}^{-3})\) is the density, \(u_{i}\)\(({\text {m}} \, {\text {s}}^{-1}\)) is the velocity, \(\Upphi \) represents the quantity of interest per unit mass, (t) is the time, and \(x_{i}\)\(({\text {m}})\) is the spatial coordinate at index i. For a 2D or 3D system this equation is repeated for each spatial coordinate [26].

$$ \frac{\partial }{\partial t}(\rho \Upphi )+ \frac{\partial }{\partial x_{i}}(\rho \Upphi u_{i})=\frac{\partial }{\partial x_{i}}\left( \Upgamma _{\Upphi } \frac{\partial \Upphi }{\partial x_{i}}\right) +S_{\Upphi }.$$
(2)

The multiphase part of the simulation uses the Volume of Fluid model. Continuity within the Volume of Fluid model is maintained by ensuring that the total volume of each phase in the system is constant, while regarding inflow and outflow of material. The Volume of Fluid model is quite dependent on the quality of the surface mesh where poor surface meshes may cause numerical diffusion of surfaces. However, since the surface behavior is not the primary objective of the study this is seen as an acceptable drawback. The continuity equation used by the Volume of Fluid method is the continuity equation for the volume fraction of a phase in the element divided by the density of that phase, as described by Eq. [3] where it is written for q-number of phases.[17, 27] For the current system the two considered phases are slag and air which then replace the q-subscript in the equation.

$$ \frac{1}{\rho _{\text {q}}} \left[ \frac{\partial }{\partial t}(\alpha _{\text {q}} \rho _{\text {q}}) + \nabla \cdot (\alpha _{\text {q}} \rho _{\text {q}} \mathbf {v_{\text {q}}})=S_{\alpha _{\text {q}}} +\sum _{p=1}^{n}({\dot{m}}_{\text {pq}}-{\dot{m}}_{\text {qp}})\right] , $$
(3)

where \( {\dot{m}} \)\(({\text {kg}} \, {\text {s}}^{-1})\) is the mass transfer from element 1 to 2, \( S_{\alpha _{\text {q}}} \) is the source term, and \(\alpha \) is the fraction of element filled with phase q.[28] The volume of fluid method is combined with the momentum equation, which is solved according to Eq. [4].[17, 27]

$$\begin{aligned} \frac{\partial }{\partial t}(\rho {\mathbf {v}})+\nabla \cdot (\rho {\mathbf {v}} {\mathbf {v}})&= {} -\nabla p+\nabla \cdot \left[ \mu \left( \nabla {\mathbf {v}}+\nabla {\mathbf {v}}^{\text {T}} \right) \right] \\&\quad +\rho {\mathbf {g}}+{\mathbf {F}} \quad \quad {\mathbf {F}}=\begin{Bmatrix} S_{x} \\ S_{y} \\ S_{z} + S_{b} \end{Bmatrix}. \end{aligned}$$
(4)

Here we see the viscosity \(\mu \)\(({\text {Pa}} \, {\text {s}})\) and the introduction of the forcing term (F) which contains different source terms depending on which momentum equation is solved. Additionally, the energy equation for the Volume of Fluid method is implemented as Eq. [5].

$$\begin{aligned}&\frac{\partial }{\partial t}(\rho E)+\nabla \cdot ({\mathbf {v}}(\rho E + p)) \\&\quad = \nabla \cdot \left( k_{\text {eff}} \nabla T - \sum _{J} h_{j} \mathbf {J_{j}} + (\bar{\bar{\tau }}_{\text {eff}} \cdot {\mathbf {v}})\right) +S_{\text {h}}. \end{aligned}$$
(5)

This formulation of the energy equation considers conduction \( k_{\text {eff}}\)\(({\text {W}} \, {\text {m}}^{-1} \, {\text {K}}^{-1})\) of heat, as well as species diffusion flux \(\mathbf {J_{j}}\)\(({\text {kg}} \, {\text {m}}^{-2} \, {\text {s}}^{-1})\) and viscous dissipation from the effective shear stress \( \tau _{\text {eff}}\)\(({\text {Pa}}). \) Since only two discrete phases, and no species are used, the summation term for the species diffusion flux can be disregarded. The source term from the latent heat of solidification (\(S_{\text {h}}\)) is defined further in Eq. [10]. The energy is defined as Eq. [6], where \( h_{\text {q}} \)\(({\text {J}} \, {\text {kg}}^{-1})\) is based on the specific heat of the considered phase and \(\nu \)\(({\text {m}}^{2} \, {\text {s}}^{-1})\) is the kinematic viscosity.

$$ E_{\text {q}}=h_{\text {q}}-\frac{p}{\rho _{\text {q}}}+\frac{\nu ^{2}}{2}.$$
(6)

The solidification model used is based on the work of Voller et al. on the enthalpy–porosity model and is applied using a sink term.[14, 17] The sink term (\(S_{x}\)) is applied to all the conservation equations discussed above, but in separate ways depending on the equation. For the momentum equations the sink term is defined as shown in Eq. [7]:

$$ S_{x,y,z}=C \cdot \frac{(1- \beta )^{2}}{(\beta ^{3}+0.001)} \cdot u, $$
(7)

where (C) is the Mushy-zone parameter (u) is the velocity of the flow in the x, y,  or z direction, and (\(\beta \)) is the liquid fraction which is 0 below solidus and 1 above liquidus. In the mushy zone region \(\beta \) varies depending on the current temperature, as in Eq. [8]:

$$ \beta =\frac{T-T_{\text {solidus}}}{T_{\text {liquidus}}-T_{\text {solidus}}}. $$
(8)

The sink term represents the changing flow behavior when the fluid partially solidifies. In the momentum equations this is done by lowering the velocity in proportion to the liquid fraction. When the liquid fraction reaches zero the material is considered fully solidified. For a 3D case, the momentum sink term is applied for all three momentum equations with an additional term in the z-direction, which is indicated in the forcing term in Eq. [4]. The additional term accounts for the buoyancy force caused by natural convection in the temperature gradient close to the solidification front and is defined in Eq. [9].[15] The buoyancy sink term is calculated based on the density (\(\rho \)), gravity g \(({\text {m}} \, {\text {s}}^{-2}),\) thermal expansion coefficient (\(\gamma \)), sensible heat \(h_{\text {s}}\)\(({\text {w}},\) and the specific heat of the fluid c \(({\text {J}} \, {\text {kg}}^{-1}\, {\text {K}}^{-1}).\)

$$ S_{\text {b}}=\frac{\rho g \gamma (h_{\text {s}}-h_{{\text {s}},{\text {ref}}})}{c}. $$
(9)

The energy equation is modified by the source term described in Eq. [10], which accounts for the energy released by latent heat during solidification, \(\Updelta \)H \(({\text {kJ}} \, {\text {kg}}^{-1})\) as defined by Voller and Prakash. The latent heat is defined as a function of temperature over the mushy zone which is a requirement for a dispersed interface between liquid and solid.[15]

$$S_{\text {h}}=\frac{\partial \rho \Updelta H}{\partial t}+{\text {div}}(\rho u \Updelta H), $$
(10)
$$\Updelta H =f(T). $$
(11)

In addition to the conservation equations the turbulence must be modeled when using the Finite Volume Method. Since a turbulent flow is based on a random process, it is complicated to describe with a numerical function. To directly calculate these motions require enormous computational efforts and is called Direct Numerical Simulation (DNS). Such a solution model is so computationally expensive it can only be used for the simplest of turbulent flows. For more complex flows a slightly less exact method can be used, namely the Large Eddy Simulation (LES) which only calculates the largest of the random swirls, or “Eddies”. However, the most widely utilized methods are the Reynolds Averaged Navier–Stokes (RANS) methods, which are based on the Navier–Stokes continuity equations. Solving a flow problem using the RANS approach yields an averaged and already smoothed interpretation of the turbulence, but requires less computational power than the DNS and LES solutions.

Since near-wall turbulence and the resolution of boundary layers is important for the formation of a mushy zone and solidified layer the choice of turbulence model will affect the appearance of the freeze-lining. The turbulence is calculated using the Baseline Reynolds Stress Model (RSM BSL) which solves more equations than the \(k{-}\varepsilon \) and \(k{-}\varepsilon \) RANS-models, but is not as elaborate as the LES-models. The RSM BSL is known to have superior resolution of boundary layers for some applications since it is based on the \(\omega \)-equation rather than the \(\omega \)-equation.[29] The RSM turbulence model has previously been used for solidification modeling for Czochralski crystal growth in melts by Ristorcelli and Lumley.[30] Additionally, the RSM turbulence model has been shown to be suitable for conjugate heat transfer problems where it performs significantly better than the standard \(k{-}\varepsilon \) turbulence model and equally to the realizable and RNG \(k{-}\varepsilon \) turbulence model.[31]

The RSM BSL turbulence model is based on the calculation of the transportation of Reynolds stresses (\(R_{ij}\)) defined as:

$$ R_{ij}=-\frac{\tau _{ij}}{\rho }. $$
(12)

With the equation for the transportation as shown in Eq. [12]:

$$\begin{aligned} \frac{DR_{ij}}{Dt}&= {} \frac{\partial R_{ij}}{\partial t} + C_{ij} \\&= {} D_{T,ij} + D_{I,ij} + P_{ij} + G_{ij} + \varphi _{ij} - \varepsilon _{ij} + F_{ij} + S_{\text {h}}, \end{aligned}$$
(13)

where C is convection, \(D_{\text {T}}\) is turbulent diffusion, \(D_{\text {I}}\) is molecular diffusion, P is stress production, G is buoyancy production, \(\varphi \) is pressure strain, \(\varepsilon \) is dissipation, F is production by system rotation, and S is a user defined source term.[17, 32,33,34]

Method

The simulations were setup and run using ANSYS Fluent 2019 R1 on a simulation computer running Linux mint 18.3 on Intel Core i9-7940X using 12 cores and 32 Gb of RAM. Also, some supplemental simulations were run on a separate computer with Linux mint 19.1 on AMD EPYC 7301 with 32 cores and 128 Gb of RAM. The geometries were constructed in Spaceclaim in the ANSYS software on windows 10 by using a geometrical symmetry at the \(y{-}z\) plane. This was used for both designs to decrease the computational demand for the simulations. Both geometries have two fluid zones and a solid insulation zone split in three zones lengthwise, all of which share topography for meshing purposes.

The mesh for both designs were constructed using the mosaic mesh function in Fluent meshing. The mosaic mesh was constructed as poly-hexcore with an inflation boundary layer to properly resolve the solidification process and the slag surface. The bulk of the domain is comprised of hexagonal cells significantly larger in size than the polyhedral cells close to the walls. The poly-hexcore approach is useful to minimize the amount of computational cells while maintaining proper resolution in the important areas.

Simulation Setup

The simulations were run with timesteps of 0.01 for 100 s of flow-time for the SR model and 300 s for the SB model to allow the solidification to become stabilized. At this time, the amount of solidified material has stabilized but has not necessarily reached a steady state. Some parameter combinations in the parameter study increase the time required for the simulation to reach a semi-steady state.

The domain for the SB and SR model are both 1 m long and lined with a 20 mm thick solid shell of alumina. The inlet of SB is located 181 mm above the lining and the inlet of the SR is 16 mm above the lining. The schematic geometry of the designs is presented in Figure 2. The boundary conditions for the simulation are set to mimic reality as closely as possible, but are slightly different between the SR and SB designs. Specifically, the boundary conditions were set according to Table I.

Fig. 2
figure 2

Schematic image of SR (a) and SB geometry (b)

Full size image
Table I Boundary Conditions
Full size table

In addition to the standard boundary conditions, the SR was tilted 15 deg by changing the direction of gravity while the SB was kept at a 0 deg incline. In the thermal boundary conditions, a thin steel shell was introduced outside the alumina refractory. The thickness of this steel shell was set to 1 cm for the SB simulation. Furthermore, for the SR a 2 cm thickness was used in the first third of the domain and a 5 cm thickness for the remaining part of the domain. This increased and varying thickness was implemented to control the amount of solidified material to prevent clogging of the runners. The boundary conditions for the inlet mass-flow and the mushy zone parameter were used as variables during the parameter study.

The SR was initialized with no slag in the runner and using a constant temperature of 1650 K in the domain. The SB was patched with slag in the lower fluid region at 1650 K to avoid simulating the initial filling of the domain. The SB was also initialized with a slow movement of 0.05 \({\text {m}} \,{\text {s}}^{-1}\) in the y-direction to prevent solidification of the slag surface by cold back-flowing air before the flow was developed.

To enable a comparison between the designs and the parameter variations, some simplifications were done to the system. The incoming air temperature is constant at 300 K. The cooling water is imitated by constant outer wall temperature of 300 K. The radiation in the system is not considered, since the slagrunner will be covered by an insulating lid. The liquid slag flows into the domain at a constant mass-flow rate and a temperature of 1650 K. The insulating wall on the bottom of the domain is split into three separate insulating zones to allow different boundary conditions in later simulations.

The material properties of slags with FeO content in the range of 90 pct by weight are not well established in literature, since they do not currently have an industrial use. For the simulation the properties were set as general values in the range of values found in literature for similar slags. Since iron oxide slags are ionic melts, which contain a mixture of Fe2+ and Fe3+ ions, it is not possible to define the precise material properties without knowing the distribution of ions. The material properties of alumina refractories also vary depending on porosity and purity of the bricks and were chosen as general values in the range found in literature. The material properties used in the simulations are shown in Table II.

Table II Material Properties [35,36,37]
Full size table

Based on the above listed conditions and properties a parameter study was performed to study the impact of variations in the material properties and process parameters on the simulation results. The focus of the parameter study was on variations in slag viscosity, inlet mass-flow, and mushy-zone parameter. The viscosity was included to account for variations in slag properties depending on temperature and composition. The inlet mass-flow was included to study the possibility to increase the mass-flow in the process, which would allow for an increased productivity of the process. In addition, the mushy-zone parameter was included to study its impact on the solidification behavior to determine what value should be used.

Mesh Analysis

A mesh sensitivity analysis was performed on both designs to act as verification of the accuracy of the simulations. The mesh sensitivity analysis was performed according to the procedure established by Richardson and further developed by Celik et al. where three meshes with significantly different element sizes are compared.[38, 39] By comparing the values of three important parameters from the simulation, a grid convergence index can be determined. The grid convergence index for refinement from a mesh to a finer mesh is calculated using Eq. [13] by using the relative error (\(e_{\text {a}}^{21}\)), the grid refinement factor (\(r_{21}\)), and the apparent order of the calculation (p).

$$ {\text {GCI}}_{\text {fine}}^{21}=\frac{1.25e_{\text {a}}^{21}}{r_{21}^{p}-1}.$$
(14)

The studied parameters were the heat flux through the insulating walls, the volume of solidified material, and the degree of freeze-lining coverage on the walls. The heat flux is calculated as three separate flux values through sections in the outer wall. This results in a negative value as the heat is leaving the domain. The volume of solidified material was determined by calculating the sum of the custom field function in Eq. [14].

$$ {\text {Solidified}}\,{\text {material}} = (1 - {\text {slag}}\,{\text {liquid}}\,{\text {fraction}})\cdot {\text {slag}}\,{\text {VoF}} \cdot {\text {cell}}\,{\text {volume}}. $$
(15)

To determine the degree of freeze-lining coverage on the walls a histogram over the slag liquid fraction in the mesh elements on the walls is constructed. All cells with a liquid fraction of slag below 0.1 are considered to be covered by frozen slag. Since all the elements on the walls have the same size the fraction of cells covered in solidified material is equal to the total freeze-lining coverage of the walls.

The GCI values for the calculated parameters in the SB are in the range of 0.3 to 5 pct and can be seen in Table III. This entails that a refinement from mesh 2 to 3 gives only slightly different values. Thus, a further refinement from mesh 3 will not give significant changes in the calculated values. Therefore, mesh 2 is an acceptable approximation of the solution at an infinitely fine mesh and will be used to calculate values in the parameter study.

Table III Stagnated Bath Mesh Analysis
Full size table

The GCI values for the Slagrunner are in the range of 1 to 6 pct as is shown in Table IV. This allows mesh 2 to be used for all later simulations without compromising the calculation accuracy by using the same logic as explained previously.

Table IV Slagrunner Mesh Analysis
Full size table

Validation Study

A validation study was performed by comparison between the experimental work by Crivits et al. and simulations in ANSYS Fluent.[17, 24] The simulation domain was constructed identical to the experimental setup used by Crivits et al. with a domain of \(300\times 100\times 50\) mm split in two with an insulating layer around all walls and a cooling plate on the side wall. The insulating layers were defined as a mixed thermal boundary condition at 293 K outside of a 3 cm thick polystyrene wall. The cooling plate was defined as a mixed thermal boundary condition at 290 K with a 1 mm copper wall. A mass flow inlet with 0.023 \({\text {kg}} \, {\text {s}}^{-1}\) flow at 323 K was used to match the 854 \({\text {mL}}\, {\text {min}}^{-1}\) at 50 \(^{\circ }\)C used by Crivits et al. In Figure 3 the validation simulation domain is shown from the side with contours of solidified material. The freeze-lining grows from the cooled wall and varies in thickness over the height and width of the domain due to the flow of hot material from the top right. The simulation replicates case 6 in the experimental work which produced a stable freeze-lining of 6.5 mm thickness after approximately 10 hours. The material properties for the liquid \({\text {CaCl}}{-}{\text {H}}_{2}{\text {O}}\) were taken from the work by Guevara et al. on 53 pct \({\text {CaCl}}{-}{\text {H}}_{2}{\text {O}}\) solidification under static conditions.[40] The solid phase was assumed to consist of entirely \({\text {CaCl}} \cdot 6{\text {H}}_{2}{\text {O}}\) with material parameters from Samanta et al.[41] All material parameters are summarized in Table V. To simplify the simulation the density and viscosity variation with phase change was disregarded. All varying material parameters were implemented as piecewise-linear variables with liquid values over liquidus temperature and solid values under solidus temperature and a linear interpolation in-between. The mesh was constructed with 8 cell inflation layer on all surfaces and 0.5 mm hexagonal cells at the cooled surface to capture the solidification front.

Fig. 3
figure 3

Solidification profile after 5 hours of 50 \(^{\circ }\)C flow (side-view)

Full size image
Table V Solidification Validation: Material Parameters
Full size table

It was found that the simulation produced a solidified layer which grew quickly for the first 3 hours of flow-time and then stabilized. However, the thickness of the freeze-lining is not uniform over the width of the domain due to the flow conditions. In the experiments the measurement is reportedly done at half the height of the domain, but no indication is given to the y-coordinate of the measurement. Thus, it is assumed it was done from the side as in Figure 3. In Figure 4 the solidification profile is shown as seen from above with the experimentally measured thickness indicated as a line through the domain at 6.5 mm. If the measurement is done from the side, the observed freeze-lining thickness will be the greatest thickness in the plane, and in such a case the thickness of the freeze-lining in the simulation is very similar to the freeze-lining in the experiments. It should, however, be noted that the time to form a stable freeze-lining was not appropriately replicated in the simulations. Thus, the simulations cannot be used to accurately predict the formation time of a stable freeze-lining in a flowing system. The flow-time in simulations required to form a stable freeze-lining will be dependent on the temporal discretization based on what size of time step is used for the transient simulation.

Fig. 4
figure 4

Solidification profile after 5 hours of 50 \(^{\circ }\)C flow (top-view) with marking for experimental measurement

Full size image

When interpreting the simulation results, the dispersed mushy zone prevents an exact measurement of the freeze-lining thickness and the results will be dependent on where the viewer considers the material solidified. The elongated areas with very high liquid fractions are not physically accurate since solidified material in such areas will break off from the freeze-lining and float away in the bulk melt.

Results and Discussion

The simulations show that both the SB and the SR design are capable of forming freeze-linings to cover the refractory material under certain conditions. The results of the parameter study with variations in viscosity, heat conductivity, mass-flow, and mushy zone parameter show that depending on the chosen material and process parameters the coverage of the freeze-lining and its thickness will vary. The shear stresses on the inner walls and the heat flux through the outer walls is also dependent on the chosen parameters. It is found that the freeze-lining in the SB design is more stable than the SR design, when varying the parameters away from optimal values.

The profile for the solidification in the standard case for both the SR and SB design is shown in Figure 5. The surface is drawn where the liquid fraction (\(\beta \)) is 0.1 which can be considered mostly solid. From this image it is apparent that the wall coverage by the freeze-lining is almost complete in the SB design but lacking in the SR design. This same behavior is represented by the wall temperature as shown in Figure 6. It is clear that the wall temperature in the SR design is on average much higher than in the SB design which is linked to the freeze-lining coverage. The reason for the wall temperature varying along the length in the SR design is that different thermal boundary conditions were used along the length to prevent clogging.

Fig. 5
figure 5

Solidified material (dark/red) and liquid surface (bright/yellow) in the SR design (a) and the SB design (b) (Color figure online)

Full size image
Fig. 6
figure 6

Wall temperature profile in the SR design (a) and the SB design (b)

Full size image

Viscosity

The parameter study shows that the viscosity of the slag has a significant impact on the thickness of the freeze-lining in both the SB and SR designs. The parameter study was performed over a large range of values varying from 0.005 to 10 Pa s, while the expected viscosity of the slag is in the range of 0.05 to 0.3 Pa s. At the outliers in the variation range some simulations experience errors and convergence issues and are thus omitted from the presented data.

In Figure 7 the thickness of the freeze-lining is plotted as a function of the value of viscosity of the slag used in the numerical simulations. The average thickness of the freeze-lining is 2 to 4 times greater in the SB design than in the SR design. The freeze-lining thickness in SB is mainly unaffected by a changed viscosity within the expected viscosity range for the slag, as is shown in Figure 7. Both above and below the expected viscosity range the thickness of the freeze-lining is increasing in the SB, which in itself is interesting. However, it is not significant to the process. The freeze-lining thickness in the SR design continuously decreases with a decreasing viscosity and the SR is thus not appropriate for slags with low viscosities.

Fig. 7
figure 7

Average freeze-lining thickness as a function of slag viscosity

Full size image

In Figure 8 the coverage of freeze-lining in percent of the contact area is shown as a function of the viscosity of the slag. The degree of coverage from the freeze-lining is more stable in the SB design than in the SR design. In the SB design the freeze-lining coverage is over 99.5 pct regardless of the viscosity value, while viscosities below 0.1 Pa s in SR cause the degree of coverage to fall below 92 pct.

Fig. 8
figure 8

Freeze-lining coverage as a function of slag viscosity

Full size image

In Figure 9 the wall shear stress is shown as a function of the viscosity of the slag. This is separated into the maximum measured value on the walls and the average shear stress on all wall elements. Increasing the viscosity increases the shear stresses on the walls in both designs. The increase is similar for both designs but the baseline shear stress is about 10 times larger in the SR design than in the SB design. This result must also be considered in conjunction with the freeze-lining coverage, as it is only the areas not covered in freeze-lining which are subjected to shear stresses.

Fig. 9
figure 9

Maximum and average shear stress as a function of slag viscosity

Full size image

Heat Conductivity

The heat conductivity of the slag is another material parameter which is uncertain for the studied slag and which is expected to have a significant effect on the solidification behavior. The heat conductivity was varied between 0.5 and 3 \({\text {W}}\, {\text {m}}^{-1}\, {\text {K}}^{-1}\) with a baseline value of 2 \({\text {W}}\, {\text {m}}^{-1}\, {\text {K}}^{-1}\) for all other simulations.

It was found that changes in the heat conductivity of the slag do not impact the freeze-lining coverage or the wall shear stresses to a significant degree. However, the increasing heat conductivity does increase the heat flux through the wall. In Figure 10 the heat flux through the different parts of the wall is presented as a function of the heat conductivity of the slag. The increase in heat flux with heat conductivity is more significant in the SB design than in the SR design. This may come as a result of the higher effective conductivity in the liquid slag of the SR design due to the faster movement and turbulence in the flow. If the effective conductivity is already high, a further increase of the heat conductivity will not have a significant influence on the total heat transfer. Increasing the heat conductivity of the liquid slag also increases the heat flux through the solidified slag as the heat conductivity is not dependent on the temperature in the present study. In Figure 11 the average freeze-lining thickness is shown as a function of the heat conductivity of the slag. It was found that increasing the heat conductivity increases the thickness of the freeze-lining. We can see the same tendency for both the SR and SB designs, although the difference is significantly larger in the SB design than in the SR design.

Fig. 10
figure 10

Heat flux as a function of heat conductivity of slag

Full size image
Fig. 11
figure 11

Average freeze-lining thickness as a function of heat conductivity of slag

Full size image

Mass Flow

Variation of the mass flow into the domain helps to determine the maximum production rate possible with each design while still maintaining a stable freeze-lining. It can also help us detect possible problems which may arise if the flow is too low due to issues in the upstream process. The mass flow was varied from 0.5 to 5 \({\text {kg}} \, {\text {s}}^{-1},\) which would allow a production of 15,000 to 150,000 tons per year depending on process uptime.

In Figure 12 the average freeze-lining thickness is shown as a function of the mass flow through the inlet. In the figure we see that an increase in the mass flow rate initially decreases the thickness of the freeze-lining in the SB design, but that the freeze-lining thickness stabilizes at higher mass flows. However, the average thickness of the SR design remains constant, albeit at a very low baseline value. This is linked to that the velocity of the slag near the wall does not increase significantly in the SR design with an increased mass flow. The angled runner already causes high velocities at low mass flows. At a 0.5 \({\text {kg}} \, {\text {s}}^{-1}\) mass flow the maximum velocity is 1 \({\text {m}} \, {\text {s}}^{-1}\) and a quadrupling of the mass flow to 2 \({\text {kg}} \, {\text {s}}^{-1}\) only increases the maximum velocity to 1.5 \({\text {m}} \, {\text {s}}^{-1}.\) For comparison, the velocity is almost 0 close to the walls in the SB design at a 0.5 \({\text {kg}} \, {\text {s}}^{-1}\) mass flow and at a 2 \({\text {kg}} \, {\text {s}}^{-1}\) mass flow the walls below the inlet and at the surface level experience velocity of 0.3 to 0.5 \({\text {m}} \, {\text {s}}^{-1}.\) At higher mass flows the turbulent flow from the inlet reaches the wall and causes high near-wall velocity, at this point the benefit of slow flow near the walls is lost from the SB design.

Fig. 12
figure 12

Average freeze-lining thickness as a function of inlet mass flow

Full size image

In Figure 13 the degree of freeze-lining coverage is shown as a function of the mass flow in the inlet. Despite the significant decrease in freeze-lining thickness for the SB design which was previously discussed, the freeze-lining coverage remains stable at above 99.5 pct until mass flows of over 3 \({\text {kg}} \, {\text {s}}^{-1}.\) However, the freeze-lining coverage in the SR design steadily decreases with increasing mass flow. Similarly to the viscosity, an increased mass flow also increases the shear stress on the walls, although not as significantly. The increase from 0.5 to 5 \({\text {kg}} \, {\text {s}}^{-1}\) multiplies the shear stresses by a factor 5. The ratio of a 10 times higher shear stress in the SR design compared to the SB design remains over the mass flow variation. The heat flux through the outer walls is increased slightly with increasing mass flow for both the SR and SB designs, giving increases of approximately 30 pct from 0.5 to 5 \({\text {kg}} \, {\text {s}}^{-1}.\)

Fig. 13
figure 13

Freeze-lining coverage as a function of inlet mass flow

Full size image

The heat flux is approximately twice as high in the SB design as in the SR design with an average value of 400 \({\text {kW}}\, {\text {m}}^{-2}\) as compared to 200 \({\text {kW}}\, {\text {m}}^{-2}.\) Using the 1D heat transfer equation with the simulation settings and the average thickness of the solidified material in SB and SR return heat flux values of 340 \({\text {kW}}\, {\text {m}}^{-2}\) for a 5 mm skull up to 646 \({\text {kW}}\, {\text {m}}^{-2}\) for a 1 mm skull. The 1D approximation is fairly accurate for the SB design where we calculate a value of 340 \({\text {kW}}\, {\text {m}}^{-2}\) instead of 400 \({\text {kW}}\, {\text {m}}^{-2}.\) However, for the SR design it is not an accurate approximation since the calculation gives 646 \({\text {kW}}\, {\text {m}}^{-2}\) while the measured value from the simulation is only 200 \({\text {kW}}\, {\text {m}}^{-2}.\) This discrepancy is mainly because the slag in the SR design does not cover the entire wall through which the heat flux is measured. The calculation of the average thickness of the freeze-lining from the amount of solidified material divided by the covered surface area also disregards any unevenness in the freeze-lining. Such an unevenness will increase the total heat flux through the wall.

Mushy Zone Parameter

The impact of the mushy zone parameter was investigated for values between 100 and 200,000 with a baseline of 10,000 used for all other simulations. It was found that the amount of solidified material in the SB design decreased when the mushy zone parameter was increased from 100 up to 10,000 where it was stable up to 100,000 and then started to increase again. In the SR design, the amount of solidified material was constant up to a mushy zone parameter of 100,000 where it started to rapidly increase. It was found that values over 100,000 quickly caused increases in the amount of solidified material. When the values approached 1,000,000, the entire studied volume was filled with solidified material.

In Figure 14 the heat flux through the different parts of the wall are shown as a function of the mushy zone parameter. The heat flux through the walls increases with increasing mushy zone parameter. Since increased values of mushy zone parameter decreased the amount of solidified material, this means that the heat flux increases when the amount of solidified material decreases, which is consistent with the previously observed heat flux behavior. In Figure 15 the degree of freeze-lining coverage is shown as a function of the mushy zone parameter. Even though a decreased amount of solidified material is observed with an increased mushy zone parameter, the freeze-lining coverage is unaffected in the SB design and increasing in the SR designs. For the SR design the increase in freeze-lining coverage is linked to the increasing source term (Eq. [7]) in the momentum equation which effectively slows the flow in areas of low temperature. Closer inspection showed that the velocity in the mushy zone in the simulation was lower for a higher mushy-zone parameter.

Fig. 14
figure 14

Heat flux through walls as a function of mushy zone parameter

Full size image
Fig. 15
figure 15

Freeze-lining coverage as a function of mushy zone parameter

Full size image

In the SB design, the freeze-lining coverage is not affected by the mushy zone parameter since the coverage is already very high. The reason the amount of solidified material decreases with increasing mushy zone parameter in the SB design may be because the low-temperature areas solidify more easily. Thus, less cooling of the bulk flow is done. The hotter bulk slag will then prevent the areas where the temperature is just below the solidification temperature from solidifying, which results in less solidified material in total, but a better quality freeze-lining. Increase in the mushy zone parameter appears to increase the stability of the freeze-lining for both designs and can thus be used to fit the simulated solidification performance to verification experiments.

Turbulence Models

For many boundary layer simulations the \(k{-}\omega \) SST model is the most common choice. It performs similarly to the \(k{-}\omega \) BSL but with a better experimental agreement for boundary-layer flows with adverse pressure gradients.[42] Such pressure gradients are not present in the current simulations, due to low flow speeds in metallurgical processes as compared to aerospace applications. Thus, the BSL approach is sufficient. However, for solidification simulations the \(k{-}\varepsilon \) turbulence model is most commonly used.[43]

The initial simulations were only functional with the RSM turbulence model and therefore it was chosen for the simulations. After refining of the geometry designs and meshes the use of the RSM turbulence model was continued for the parameter study. After finishing the parameter study, the turbulence models were again evaluated and it was found that all the tested turbulence models were functional in the updated geometry and mesh. It was found that the \(k{-}\omega \) SST and BSL turbulence models returned very similar values to the RSM turbulence model for all the studied parameters. The \(k{-}\varepsilon \) turbulence model returned heat flux values approximately 40 pct lower than for the RSM and \(k{-}\omega \) turbulence models. The \(k{-}\varepsilon \) model also showed twice the amount of solidified material compared to the RSM turbulence model. The transition k–kl turbulence model was found to give values somewhere in-between \(k{-}\varepsilon \) and \(k{-}\omega \) BSL turbulence models. From this analysis, it was concluded that the RSM turbulence model may not be the most appropriate choice for this study considering the simulation overhead.

However, for the purpose of the current study, it can be concluded that there is a sufficient formation of freeze-lining, regardless of which turbulence model is being used. The RSM turbulence model can be seen as a conservative approximation to the solidified layer thickness.

Conclusions

The study aimed to determine a way to minimize the refractory wear of a slag transfer tube to prolong the refractory life when transporting a highly corrosive slag of 90 pct FeO in the IronArc process. With CFD simulations it was found that a design of the transfer tube which minimized the near-wall turbulence by slowing down the flow of slag experienced the lowest wall shear stress values. By introducing water-cooling a freeze-lining was able to form on the inside of the refractory wall which further reduced the wall shear stress on the refractory.

The fast flowing transfer tube design (SR) and the slow flowing design (SB) were studied in terms of heat flux, shear stress on the walls, freeze-lining coverage, and freeze-lining thickness. The parameter study of the material properties of the flowing slag varied the viscosity, heat conductivity, mushy zone parameter, and mass-flow inlet in the two designs.

It was found that the slow-moving SB design is more affected by changes in parameters than the SR design, but also always maintains a sufficient freeze-lining to protect the refractory wall. The SR design is mostly unaffected by changes. However, in some cases the protective property is entirely lost, which would expose the refractory material to a corrosive and mechanical wear by the slag. This has been attributed to the fact that the freeze-linings are very stable close to the refractory wall, but require favorable conditions to grow thicker.

An increased slag viscosity will reduce the acceleration of the flow due to gravity as well as help to dampen the turbulence from the inlet or near-wall interactions, which will enable an increased solidification of the slag. For the viscosity range of 0.05 to 0.3 Pa s, which is expected in the IronArc slag, the freeze-lining is stable in the SB design and in some regions in the SR design. Overall the properties of the freeze-lining do not change significantly in this range.

The heat conductivity of the slag is found to have a minor impact of the stability of the freeze-lining. However, it is linked to the heat flux in the wall in the SB design, where an increase in heat conductivity increases the heat flux through the wall. The mass flow rate through the inlet is limited to approximately 3 \({\text {kg}} \, {\text {s}}^{-1}\) after which the coverage of the freeze-lining starts to quickly decrease. At lower flow speeds the SB will experience increasing amounts of solidified material, but not to the point of solidifying the entire transfer tube.

The enthalpy porosity model was found to be efficient in slowing the flow in the mushy region and the magnitude of this effect was increased with increasing mushy zone parameter. A modification of the mushy zone parameter can be used to fit the simulation results to verification experiments of solidification to best represent the specific solidification behavior of the studied material and flow case.

The verification simulations achieved good similarity to the experimental values using a mushy zone parameter of 10,000. The validation study also shows that the temporal similarity in the formation of a freeze-lining is unreliable and dependent on the time discretization by the chosen time step. The final thickness of the freeze-lining was similar in the simulation and the validation experiment, but the steady state was reached much faster in the simulations. This should be considered when modeling freeze-lining formation in industrial applications using the enthalpy–porosity model. The RSM turbulence model produced similar results to the \(k{-}\omega \) SST model. However, it is likely not the best choice for solidification simulations due to the higher computational power requirement, instead the \(k{-}\omega \) SST model should be used.

The overall conclusion is that it is possible to use the SB design of the transfer tube with a freeze-lining in order to protect the refractory from the corrosive slag. This will enable continuous production in the IronArc process.