CFD Modeling of the Melting Behavior of a Fayalitic Slag

The melting behavior of a fayalitic slag was investigated with computational fluid dynamics and validated with experimental data. A heating microscope was employed to perform a typical ash melting behavior (AMB) test, in order to gain information about the melting behavior of the slag. The numerical approach is based on using the volume of fluid (VOF) method, in order to track the interface between the slag and the gas phases, and on the enthalpy–porosity model, in order to simulate the slag as a phase change material (PCM). The combined VOF-PCM model was therefore applied to perform a parameter study on the mushy zone constant Amushy\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{\rm mushy}$$\end{document}, which enters the equations of the enthalpy–porosity sub-model and is an a priori unknown. The goal of the investigation was to determine the correct value for Amushy\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{\rm mushy}$$\end{document}, by comparing the numerical results to the AMB test.


INTRODUCTION
The solidification of slag is a physical phenomenon of great importance for high-temperature processes in ferrous and non-ferrous metallurgy. In facts, the formation of a freeze lining at reactors walls is a key design aspect for several operating furnaces, from flash smelters to top submerged lance (TSL) furnaces. [1][2][3] The presence of solid slag on the walls acts as a protective layer for the refractory lining, therefore increasing its life. Nikolic et al. reported an increase in the campaign life of a TSL furnace from 65 to 215 weeks, 4 obtained with a specific development project focusing on the reactor cooling and operating point. One modern type of cooling system involves installing copper cooling elements directly in contact with the liquid slag. 5 The intention is to establish a steady layer of frozen slag, in order to reduce the heat losses and prevent wear of the side walls. 6 In TSL furnaces, slag also solidifies on the outer surface of the injection lance; here, the slag splashed from the bath freezes, cooled by the inner flowing process gas. This also allows the lifetime of the lance to be extended. 7,8 The direct observation of slag lining layering from industrial or even pilot furnaces can still be very challenging and slow in terms of data acquisition. Thus, numerical modeling can be an advantageous tool to support the design and optimization of these processes, and some works in the literature show the growing interest from the scientific and industrial community. Joubert et al. applied finite element analysis to calculate the temperature distributions of copper cooling elements installed in a TSL furnace. 5,6 They compared the use of three different coolants, namely water, mono-ethylene glycol (MEG), and ionic liquid. The numerical results showed that, for both MEG and ionic liquid, the temperatures of the copper element were above the design limit, whereas safer operating conditions were obtained with water. Solnordal et al. developed a detailed 1D model of the heat transfer process for the lance of a pilot TSL furnace. 9 The model included a simplified approach for the slag solidification, which was considered if the slag temperature went below the solidus temperature, ignoring the so-called mushy region between the liquidus and solidus. Nevertheless, they were able to predict the thickness of the solid layer for most of its length over the lance. (Received November 17, 2021; accepted January 25, 2022; published online February 28, 2022) JOM, Vol. 74, No. 4, 2022 https://doi.org/10.1007/s11837-022-05185-4 Ó 2022 The Author(s) Solidification phenomena can also be included in computational fluid dynamics (CFD) models, where the dynamic interactions with the flow and energy fields can be taken into account. Unlike the case of pure substances, the phase transformations of multicomponent materials, such as slags or metallic alloys, are not iso-thermal and involve the presence of a liquid-solid region between the liquidus and solidus points. The solidification of a liquid mixture, l, to produce a solid mixture, a , can therefore be generalized as l ! l þ a ! a. The region between the liquidus and solidus temperatures (l þ a) is defined as the mushy region, due to the co-presence of liquid and solid fractions. 10 These transformations are typically implemented into CFD codes using the enthalpy-porosity approach, a model proposed by Poirier,11 and applied to fluid dynamics by Voller et al. 12,13 The model assumes the mushy region to be a porous medium, whose permeability is a function of the liquid fraction b, and introduces a momentum source S defined as: where u is the fluid velocity, e is a small number to avoid division by 0, and A mushy is the mushy region constant. It must be noted that this last parameter significantly influences the rate of phase transformation and is an a priori unknown variable. From a physical perspective, it represents the degree of penetration of the liquid convection into the porous medium, and is strictly dependent on the morphology of the material. 13 In some works from the literature, this aspect of the model is neglected, making the quality of the results questionable. Ge et al. used CFD to simulate the process of horizontal single belt casting for a metallic alloy. 14 They did not provide any discussion on the determination of the constant and applied a value of 10 5 kg/m 3 s, set by default in the commercial software ANSYS Fluent. Eickhoff et al. re-implemented the enthalpy-porosity model in ANSYS Fluent, 15 in order to take into account the non-linearity of the liquid fraction as a function of the temperature for the material they studied. They reported the importance of A mushy for determining the porous permeability, but did not provide any information on the values used. Tan et al. also did not provide any information on A mushy in their work about solidification in the presence of the Marangoni effect. 16 In some other works, 17 the mushy region constant is estimated as a function of the secondary dendrite arm spacing, indicated here with k, following the work of Poirier. 11 Yang et al. defined A mushy ¼ 180=k 2 and applied the model to calculate the solidification process in continuous casting with electromagnetic stirring. 18 They were able to predict the temporal evolution of the solidified shell thickness. In his doctoral thesis, Giesselmann explains in detail the relationship between the mushy region constant, the porous permeability and the length scale of the dendrite cells. 19 Knowing the geometrical structure of the dendrite for the material examined, he carried out a detailed CFD simulation of the liquid phase flowing through the solidified material, calculating the effective permeability.
Another approach found in the literature is based on conducting a simple phase change experiment with the material under examination, along with a CFD parametric study on A mushy , in order to reproduce the macroscopic physical behavior of the phase transformation. Works can be found for phase change materials (PCM), such as lauric acids and gallium. [20][21][22] Nevertheless, to the authors' best knowledge, no information is available for fayalitic slags, the subject of this study. The values of A mushy strictly depend on the material, and the present work aims to evaluate the phase transformation properties for a fayalitic slag, which has already been the subject of previous works by the authors. 23,24 Correctly modeling the slag solidification behavior is indeed necessary to develop a CFDbased thermal analysis of a TSL furnace, as studied in Ref. 25.

EXPERIMENTAL SETUP
A heating microscope, combined with the image analysis software EMI III (v.3.0.7), was employed to perform an ash melting behavior (AMB) test in line with the DIN 51730. 26 This measurement technique is typically applied to determine the softening and melting temperatures of ash samples by graphically monitoring the evolution of the sample shape as a function of the temperature. Figure 1 shows a schematic of the experimental apparatus. The ash is pressed into a cylindrical shape, with a diameter and height of 3 mm, and inserted in the furnace with a support tube made of alumina. The sample lies on a small alumina plate (10 mm 9 10 mm), placed on the support tube, and is inserted at t ¼ 0 s, when the temperature inside the microscope is about 300°C. A heating rate of 60 K/min is then applied until a temperature of 800°C is reached, after which it is reduced to 10 K/min. During the heating phase, an external camera records changes of the shape of the ash cylinder every 5 s. A platinum-rhodium-platinum thermocouple is located inside the support tube to measure the temperature directly beneath the sample. 26 This temperature will be referred here as the ''support temperature'', since it measures the temperature inside the support tube. Another thermocouple of the same type measures the temperature of the furnace, referred here as the ''furnace temperature''. This is integrated inside the heating conductor and is not reported in Fig. 1 for the purpose of simplification. A mixture of CO and CO 2 is employed as a carrier gas, which flows continuously inside the furnace. The test is considered to have been completed, after the ash has totally melted. Figure 2 shows the profile of the furnace and support temperatures over time. It can be seen that, after 15 min, the furnace temperature and the support temperature differ by approximately 50 K. From that time on, both temperatures increase with approximately the same slope until the end of the experiment.
The change in the sample shape during the test is shown in Supplementary Figure S1 (refer to online supplementary material). During the first 40 min, no modifications can be observed. Afterwards, the cylinder begins to deform, first slightly shrinking and then growing in diameter and height. Hence, the slag assumes a semi-spherical shape until, after around 49 min, it is completely liquefied. In the first pictures recorded by the camera, the profile of the slag is seen to be rectangular. Since the diameter and the height of the slag cylinder are equal, additional slag material must have been on top of the cylinder. This fragment probably remained after the stamping process, and is not considered in the geometry used for the numerical investigations.

General Description
The commercial software ANSYS Fluent (v. 19.4) was employed to perform the CFD calculations. The approach used to simulate the AMB test is based on using the volume of fluid (VOF) method, as a multiphase model, and the enthalpy-porosity model for the phase change. Hence, the numerical model considers a set of equations including the continuity, momentum, and energy balance, void fraction advection, the enthalpy-porosity equations and the radiation model. The flow develops in a laminar regime, therefore turbulence modeling is not included.
In the VOF formulation, only one set of the Navier-Stokes equations is applied: where u is the velocity vector, p is the pressure, f g is the gravitational force, and f r is the interfacial force due to the surface tension. Here, a is the liquid volume fraction, defined as:  At the interface, where its value varies between 0 and 1, a geometrical reconstruction scheme is applied to reconstruct the interface using the piecewise-linear interface construction approach. 27 After the reconstruction, a is advected in the given velocity field following the advection equation: The formulation of the enthalpy-porosity model that is the subject of this study is worth closer explanation. When a material undergoes a phase change, its enthalpy can be computed as the sum of the sensible enthalpy, h, and the contribution of latent heat, DH: and Here, L is the heat of fusion of the material and b is the liquid fraction of the liquid-solid mixture at a specific temperature between the liquidus and the solidus temperatures. It is defined as: For a combined VOF and enthalpy-porosity model, the energy equation can therefore be written as: where a is the VOF volume fraction, and H includes the information about the latent heat. 28 As mentioned above, the enthalpy-porosity model considers the mushy region to be a ''pseudo'' porous medium, whose permeability is a function of the liquid fraction b. The porosity varies between 0 and 1. In the case of material solidification, when the fluid in a cell is completely solidified, the porosity goes to 0, as well as the velocity field. 28 The consequent momentum sink in the mushy region is taken into account by adding a source term S to the momentum equation, expressed as: The weight of the mushy region constant A mushy and the need for an appropriate evaluation is evident, since it proportionally affects the source term. The P-1 model is applied to include radiation in the heat transfer. The mathematical formulation can be found in Ref. 28. A thermodynamic calculation was performed in FactSage, to determine the liquidus and solidus temperatures, the liquid fraction-temperature curve and the heat of fusion of the slag examined. The slag composition is reported in Table I. This originates from a secondary copper slag converter, and has already been the object of several studies. 23,24 FToxid and FTmisc databases were applied, allowing the formation of all possible phases and two-and three-phase immiscibility. A stepwise calculation of equilibrium points was performed by varying the temperature and evaluating the liquid fraction at each of these points. The resulting liquid fraction-temperature curve is shown by the blue line in Fig. 3, determining the liquidus and solidus points, where b was respectively 1 and 0. However, from 750°C to around 1000°C the liquid fraction has a constant value of about 0.01, and is practically negligible. The temperature-dependency of b was therefore simplified by a linear approximation in the range with significant variation (red line), identifying the effective solidus temperature at 1003°C. The liquidus temperature is equal to 1185°C and the heat of fusion L is 657.3 kJ/kg. The linear approximation of b is also consistent with the definition of bðTÞ in the enthalpy-porosity model in ANSYS Fluent (see Eq. 9). This linear dependency should always be verified for the material under examination, and, if not fulfilled, the implementation of a new function bðTÞ is necessary. 15 The temperature-dependency of the slag's thermo-physical properties was taken into account by implementing user-defined functions. The density, viscosity, and surface tension were measured experimentally, while the specific heat and thermal conductivity have been modeled as proposed by Mills. 29,30 More details can be found in Ref. 23, where the authors focused on the temperaturedependency of the slag's physical properties and their influence on the modeling of slag flows. The contact angle formed between the liquid slag, the gas atmosphere, and the plate is an important parameter to take into account, because it directly affects the shape of the droplet at the contact point. This was measured during the AMB test, and the function of its temperature-dependency was implemented in the CFD model with a UDF. Figure 4a and b shows the definition and the thermal behavior monitored during the test. Figure 5 shows the three-dimensional computational domain, with the geometry of the heating microscope accurately reproduced as during the AMB test. At the inlet, a gas mixture of 55 vol% CO and 45 vol% CO 2 enters with a velocity of 0.001943 m/s and leaves the furnace from the outlet. The inner furnace wall is modeled with a wall boundary condition, and the temperature is set according to the information resulting from the test, see ''Furnace temperature'' in Fig. 2. The support and the plate made of Al 2 O 3 are modeled as solid regions in order to calculate the conjugate heat transfer at their boundaries, which is otherwise unknown.

Modeling Strategy
The full 3D simulation of the AMB test is, however, still too demanding in terms of computational effort. The reasons are related to the different geometrical scales of the furnace and the probe. In fact, considering that the diameter and height of the probe are 3 mm, a grid resolution of at least 10 cells per diameter leads to a mesh size of 0.3 mm. The time step that results if the convergence is kept stable and the Courant-Freidrichs-Lewy condition (Courant number < 1) is satisfied, would make it unfeasible to simulate 50 min of real time, which is the duration of the test. To overcome this inconvenience, the investigation is broken down into two steps: 1. 3D simulation of the heat transfer inside the experimental furnace 2. 2D axi-symmetric simulation of the slag melting behavior.   Since the slag probe is cylindrical, a 2D axi-symmetric approach can be applied, restricting the domain just to the region where the sample is located. This expedient drastically reduces the computational effort, allowing the AMB test to be reproduced. Nevertheless, the thermal boundary conditions that apply to this sub-domain are unknown. Hence, a detailed 3D simulation of the heat transfer inside the heating microscope is first carried out, in order to determine the temperature distributions in the furnace that must be applied as an input for the second step. These are validated by direct comparison with the experimental data from the AMB test. The modeling strategy is summarized here, with specific information for each of the two steps.

3D Simulation of the Heat Transfer
The VOF equation is not applied in the first modeling step, and the slag sample is modeled as a solid material. Given the size of the sample, the heat of fusion is assumed to be a negligible contribution to the determination of the temperature distribution in the furnace. The goals of this analysis are to validate the temperature results with experimental data and to extract the transient development of the temperature inside the slag sample. Figure 6a and b offers an insight into the computational grid used for this analysis, generated with the software ANSA Beta CAE (v.18.01.1). The mesh consists of 105,454 polyhedral cells, and has a refinement region around the sample with a resolution of 10 cells per diameter of the sample. This was set as a trade-off between the level of detail and the feasibility of simulating the elapsed times. As mentioned, the heat transfer inside the furnace is calculated by including the conjugate heat transfer with the solid regions, and the radiative heat transfer using the P-1 model.
The numerical setup is summarized in the Table -II. During the simulation, the temperature inside the alumina support is monitored at the same location where the thermo-couple is placed in the experiment, and compared to the ''Support temperature'', shown in Fig. 2. Alongside this verification, the temperature is also monitored at the center of the slag sample, and is used as input data for the 2D axi-symmetric simulation.

2D Axi-symmetric Model of Slag Melting
At this step, it is assumed that the temperature is homogeneously distributed inside the slag sample, and that the carrier flow has no effect on the melting process of the sample. The aim of the 2D model is to simulate the melting behavior with the VOF-PCM method, and to perform the parametric study on A mushy .
The computational domain for this simulation step is shown in Fig. 7. Here, a structured grid of quad elements is employed, with a level of refinement of 10 cells per sample diameter. The temperature of the slag sample, recorded during the detailed 3D simulation of the heat transfer, is superimposed in the domain, according to the assumption of homogeneous thermal distribution in the sample. Eventually, the melting phenomenon is computed with the combined VOF-PCM model, described above. The setup of the 2D axi-symmetric approach is summarized in Table III.

3D Analysis of the Heat Transfer
In order to validate the calculation of the heat transfer inside the heating microscope, a parametric study of the setup was required. In particular, the study focused on the emissivity coefficient e of  alumina, the material from which the support and the plate under the sample are made, and which was found to have a major influence on the temperature distributions. The values of e for the specific Al 2 O 3 material used in the heating microscope are not known. Jones et al. reported that emissivity coefficients for oxidized aluminum surfaces vary in the range from 0.03 to 0.2. 31 Several simulations were therefore performed varying the emissivity e of the support and the plate in this range. Figure 8a shows the calculated time-dependent support temperature profiles in comparison with the experimental data. The influence of the alumina emissivity is clear, as also expected. Nevertheless, the choice for a correct prediction demands further considerations. While the curve with e ¼ 0:2 correlates well with the experimental curve during the first 10 min of the simulation, after 30 min, the curves with e ¼ 0:03, 0.05 have the strongest correlation. It should be noted that, for the purpose of the melting analysis, the initial phase of the experiment is not relevant. In fact, as shown in Fig. 8b, the relevant phase of the experiment starts from temperatures around 850°C, corresponding to a time of 30 min, after which the slag sample starts to deform, as indicated by the changes in the projected surface area and the form factor. This allows the relevant observation window to be reduced to 850-1200°C, or 30-45 min. As shown by the changes in the surface area, the sample shrinks until reaching the softening temperature, then expands to the semi-spherical shape, and is eventually flattened after complete melting.
The selection of the proper value for e is based on the evaluation of a relative error with respect to the experimental curve, weighted over the selected time range from 30 min to 45 min. The relative error is calculated as: The configuration with e ¼ 0:05 proved the most accurate with T err ¼ 0:98% and was therefore selected for the final setup of the 3D simulation to monitor the temperature inside the slag sample. This result is shown in Fig. 9, where the CFD   ''Support temperature'' (blue) is compared to that measured experimentally (red). The curve related to the ''Sample temperature'' (orange) represents the time-dependent temperature inside the slag sample, which is then superimposed when performing the 2D axi-symmetric VOF-PCM simulation. As expected, this lies in between the ''Furnace temperature'' (black) and the ''Support temperature'' (blue).

Parametric Study on A mushy
Once the time-dependent temperature in the sample had been determined, the melting behavior of the slag was simulated with the VOF-PCM-based approach in the 2D axi-symmetric domain. This model configuration was therefore applied to perform a parametric study on the parameter A mushy , which was the objective of the study.
After preliminary investigations, the correct value of A mushy was expected to be in the range of 10 10 to 10 13 kg/m 3 s. In order to check the influence of the melting behavior in this specific range, the first simulations were performed with A mushy = 10 10 , 10 11 , 10 12 and 10 13 kg/m 3 s. In each case, the cross-section profiles were captured at different times and compared to the corresponding experimental snapshots. The comparison of the temporal sample shapes is in fact sufficient to prove the correct replication of the macroscopic behavior of the phase change. In Fig. 10, it can be seen that the correct value for A mushy must be close to 10 12 kg/ m 3 s. Thus, four more simulations were performed, narrowing the range of investigation in that area. Comparing the shapes of the melting slag, it becomes evident that the configurations with A mushy from 7:5 Á 10 11 to 2:5 Á 10 12 kg/m 3 s follow the melting behavior throughout the duration of the experiments. Indeed, the slag melts quickly at lower values and too slowly at the highest values. Slight differences between the three configurations are still observed at the time of 47 min. Figure 11 shows a more detailed comparison, where the gas-slag interfaces from the three cases overlap with the experimental profile. Here, the CFD interface is extracted as an iso-contour line at a void fraction equal to 0.5. The results indicate that the appropriate value of A mushy for the fayalitic slag under investigation is equal to 10 12 kg/m 3 s.

Verification with Other Sample Shapes
The CFD simulation of the melting slag was repeated with a different slag sample geometry, in order to verify that the value of A mushy ¼ 10 12 kg/ m 3 s was geometry-independent and only a characteristic of the material. A new AMB test was performed with a slag sample pressed in the shape of a truncated cylindrical cone. The result of this verification study is reported in Fig. 12, showing good agreement between the experiment and the simulation. This confirms the reproducibility of the analysis and the value of the mushy region constant that was previously identified.

CONCLUSION
The melting behavior of a fayalitic multi-compound slag was investigated by means of CFD modeling. The numerical results were validated by experimental measurements carried out with an AMB test in a heating microscope, in line with the DIN 51730. 26 The main focus of the analysis was on the determination of an appropriate value for the mushy zone constant A mushy , as this is a governing material-specific parameter for modeling a phase changing material following the enthalpy-porosity approach. 28 Due to the high computational requirements, the simulation of the AMB test was split into two steps: a detailed simulation of the heat transfer in the microscope in 3D and a 2D axi-symmetric simulation of the melting slag. This was necessary to determine the time-dependent temperature profile inside the slag sample, which could then be applied as input data to simulate the phase change.
A parametric study was performed on the mushy zone constant A mushy , varying its value in a range from 10 10 to 10 13 kg/m 3 s, and comparing the CFDcalculated shape of the melting sample with experimental snapshots taken at different time instants during the test. A mushy region constant of A mushy = 10 12 kg/m 3 s was found to correctly predict the shape changing behavior for this specific fayalitic slag. The consistency of the result was verified by reproducing the AMB measurement and the simulation with a different slag sample geometry, namely a truncated cylindrical cone. A mushy = 10 12 kg/m 3 s is hence an intrinsic value for the specific slag under examination.
The reader should note how much the final result differs from the default values found in commercial software, and often used in works published in the literature. It is important to stress once again that the value of A mushy is strongly connected to the physical properties of the material, and therefore evaluations such as the one presented in this paper should be carried out for each material of interest. Similar values are, however, expected for fayalitic slags and can be applied to include the enthalpyporosity model in the CFD simulation of non-ferrous metallurgical processes, where slag solidification plays a major role in terms of heat transfer and the development of the fluid flow.

ACKNOWLEDGMENTS
The authors would like to thank the German Federal Ministry of Education and Research (BMBF) which has funded this research within the framework of the Center for Innovation Competence

FUNDING
Open Access funding enabled and organized by Projekt DEAL.

CONFLICT OF INTEREST
On behalf of all authors, the corresponding author states that there is no conflict of interest.

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 h ttp://creativecommons.org/licenses/by/4.0/.