Numerical Investigation and Experimental Validation of Motion and Distribution of Nonmetallic Inclusions in Argon Protection Electroslag Remelting Process

A transient coupled mathematical model has been proposed to predict the motion and distribution of nonmetallic inclusions in the electroslag remelting process. The electromagnetic field, temperature distribution, flow pattern, solidification, and inclusion behavior were solved simultaneously. The Euler-Lagrange approach was used to describe the interaction between continuous phase and inclusions. The discrete phase model (DPM) was adopted to track the particles acted by a number of forces. An experiment under argon protection was implemented, and it agrees well with the simulation. The effect of inclusion diameter was also discussed. The experiment result shows that the typical inclusions in final ingot are Al2O3 core surrounded by an outer sulfide layer and pure MnS inclusions. The total area of inclusions increases from center toward ingot surface, while the average equivalent diameter decreases. The simulated result indicates a denser distribution of inclusions locations at the region ranging from 0.7 to 0.9 radii. Few inclusions are entrapped near the mold due to the strong descending flow. Moreover, the inclusion amount increases with height. The total amounts of inclusions with diameter of 2 μm, 5 μm, and 10 μm in ingot are 2199, 1871, and 1704, respectively. The imposed buoyancy and floatation contribute to this phenomenon.


Introduction
Nonmetallic inclusions in steel cause serious defects and affect the mechanical properties of final ingots negatively [1].Electroslag remelting (ESR) process, which could effectively remove the nonmetallic inclusions and improve the steel cleanliness, has been widely used to produce high quality steels and alloys for applications, such as energy, aeronautics, and tooling [2,3].With the increasing demand of high performance alloys, a deeper understanding of the mechanism of inclusion behavior during ESR process must be done.
As shown in Figure 1, the molten slag is electrically heated by an alternating current (AC) to melt the consumable electrode in the ESR process under argon protection.The molten metal passes through the slag and feeds a liquid metal pool, which solidified directionally in a water-cooled mold.There are two main sources of nonmetallic inclusions in ESR ingots, one is the original inclusion of the electrode and the other is the newly precipitated inclusion at the solidification front during the solidification.The former is supposed to be mainly removed at the liquid metal film and the growing droplet on the electrode tip through washing out mechanically and dissolution [4].A further reduction of original inclusions is found during the falling of metal droplet, despite of the very short time about 0.1-0.15s.Few original inclusions from electrode get into the liquid metal pool, but precipitation and growth of new formed nonmetallic inclusions occur in two-phase region between liquidus and solidus temperature [5].The motion, removal, precipitation, and final distribution of nonmetallic inclusions are generally affected by the thermodynamic and kinetic factors.time about 0.1-0.15s.Few original inclusions from electrode get into the liquid metal pool, but precipitation and growth of new formed nonmetallic inclusions occur in two-phase region between liquidus and solidus temperature [5].The motion, removal, precipitation, and final distribution of nonmetallic inclusions are generally affected by the thermodynamic and kinetic factors.Many metallurgists experimentally investigated the refining parameters and processes on the characteristics of inclusions in ESR ingots.Dong et al. focused on the effect of fluoride containing slag on the inclusions in ESR ingots [6].The results determine that the inclusions of ESR ingot using multicomponent fluoride containing slag are MgO-Al2O3 inclusions, but Al2O3 inclusions when binary fluoride containing slag is adopted.Yang et al. studied the generation mechanism of TiN inclusion for GCr15SiMn during the ESR process [7].The results show that TiN inclusion that is found in final ingot is regenerated during solidification of metal, not the one from consumable electrode.A higher content of SiO2 in slag leads to a lower TiN content in ESR ingot due to the Si/Ti exchange reaction.Six ingots were remelted by Paar et al. to analyze the effect of electric power supply on the type and content of non-metallic inclusions in ESR ingots [8].They noticed that the control of oxygen and oxide inclusions is coherent.The results indicate that the one with 1 Hz AC power supply has a better cleaning effect than the ones with both direct current (DC) polarities.The differences of size and type of inclusions under different electric power supply can be explained by related electrolytic reactions.
Due to the limited information provided by opaque reactor of experiments, the behavior of nonmetallic inclusions during the refining process is still not clear.Numerical simulation thereby has been widely used to investigate the complicated phenomenon [9].Kelkar et al. numerical tracked the inclusions with different size and density during the ESR process for the first time [10].A temperature-dependent dissolution rate of inclusions is considered.The simulated results show that large inclusions are mainly affected by buoyancy, while smaller ones move with the fluid flow.However, the inclusions were added after the flow field reached steady state, which is not in accord with the fact.In other words, the flow field at a certain time determines the whole trajectory of inclusions from the addition to the static state in solidified ingot.Kharicha et al. established a twodimensional (2D) numerical model to study the origin of radial distribution of nonmetallic inclusions in ESR ingot [11].The frequency of the inclusions was found to increase from the centre toward the ingot surface.But the periodic formation and dripping of molten droplet were simplified, which affects the flow field significantly.
As mentioned above, few attempts have been made to numerically investigate the motion and distribution of nonmetallic inclusion in ESR process coupled with transient multi-physical fields, which restricts the wider application of ESR technology.In this case, a transient 2D axisymmetric mathematical model was developed to analyze the mechanism of inclusion behavior during the ESR process.The behavior of original inclusions from the electrode was simplified due to an extremely high removal ratio.The newly generated inclusions at the solidification front were carefully discussed.The Euler-Lagrange approach was used to describe the interaction between continuous  [8].They noticed that the control of oxygen and oxide inclusions is coherent.The results indicate that the one with 1 Hz AC power supply has a better cleaning effect than the ones with both direct current (DC) polarities.The differences of size and type of inclusions under different electric power supply can be explained by related electrolytic reactions.
Due to the limited information provided by opaque reactor of experiments, the behavior of nonmetallic inclusions during the refining process is still not clear.Numerical simulation thereby has been widely used to investigate the complicated phenomenon [9].Kelkar et al. numerical tracked the inclusions with different size and density during the ESR process for the first time [10].A temperature-dependent dissolution rate of inclusions is considered.The simulated results show that large inclusions are mainly affected by buoyancy, while smaller ones move with the fluid flow.However, the inclusions were added after the flow field reached steady state, which is not in accord with the fact.In other words, the flow field at a certain time determines the whole trajectory of inclusions from the addition to the static state in solidified ingot.Kharicha et al. established a two-dimensional (2D) numerical model to study the origin of radial distribution of nonmetallic inclusions in ESR ingot [11].The frequency of the inclusions was found to increase from the centre toward the ingot surface.But the periodic formation and dripping of molten droplet were simplified, which affects the flow field significantly.
As mentioned above, few attempts have been made to numerically investigate the motion and distribution of nonmetallic inclusion in ESR process coupled with transient multi-physical fields, which restricts the wider application of ESR technology.In this case, a transient 2D axisymmetric mathematical model was developed to analyze the mechanism of inclusion behavior during the ESR process.The behavior of original inclusions from the electrode was simplified due to an extremely high removal ratio.The newly generated inclusions at the solidification front were carefully discussed.The Euler-Lagrange approach was used to describe the interaction between continuous Metals 2018, 8, 392 3 of 20 phase and inclusions.The electromagnetism, flow pattern, heat transfer, solidification, and inclusion behavior were solved simultaneously.Besides, a laboratory experiment under argon protection was implemented, and both the electrode and ingot were analyzed to validate the model.

Model Description
In order to establish the model computationally, the following assumptions are invoked in the current study: (a) The simulated domain includes the molten slag, liquid metal pool and solidified ingot; (b) The leakage current toward the mold is not considered.The slag and metal are assumed to be electrically insulated from the mold [12]; (c) The densities of both phases and the electrical conductivity of slag are temperature dependent, while the other properties are supposed to be constant [13]; (d) The start-up stage and the feeding stage of ESR process are neglected; and, (e) The behavior of inclusions in slag is not considered.

Electromagnetism
Maxwell's equations are used to describe the electromagnetic field created by an AC current: The current displacement which is lower than the electric conduction can be ignored [14].The constitutive equations for the electromagnetic fields are expressed by: The electric potential method is adopted to solve the electromagnetic field, while the magnetic potential vector is introduced to solve the magnetic field [15]: Finally, the time-average Lorentz force and Joule heating are calculated, and are incorporated into momentum and energy conservation equations, respectively:

Fluid Flow
The two-phase incompressible fluid flow in ESR mold is modeled by the continuity equation and Navier-Stokes equation: ε , where C µ = 0.0845 [16].→ F st is the surface tension described by the continuum surface force model that was developed by Brackbill et al. [17].→ F b is the buoyancy force that was determined by the Boussinesq approximation [18].→ F damp represents the damping force in the mushy zone.→ S is a source term accounting for the momentum exchange between the melt and inclusions.
The volume of fluid (VOF) approach is employed to track the slag/metal interface [19]: where α represents the volume fraction of the metal phase.

Heat Transfer and Solidification
The energy conservation equation is described by [20,21]: where k e f f is the effective thermal conductivity.L is the latent heat of fusion of metal.C Joule is the corrected coefficient of Joule heating accounting for the loss of leakage current toward the mold.E is the internal energy of the mixture phase, which is related to the specific heat and temperature.
The second term on the right side represents the latent heat released during solidification.The third term on the right side represents the Joule heating that maintains the heat balance in the slag layer.
The liquid fraction f l in Equation ( 12) is estimated as [22,23]: where T l and T s are the liquidus and solidus temperatures related to local solute concentrations.The permeability in the metal mushy zone determines the flow resistance and can be modeled with Blake-Kozeny equation [24,25]: where λ represents the characteristic length, and the default is the secondary dendrite arm spacing of the ingot.The damp force in the mushy zone mentioned in Equation (10) therefore can be written, as follows:

Inclusion Motion
The nonmetallic inclusions are spatially discontinuous, and can be simplified to discrete spherical geometries, which could be tracked by the discrete phase model (DPM) [26].The flow of melt is described by the Euler equation, as mentioned above, while the motion of inclusions is described by the equation of Lagrange form.As illustrated in Figure 2, the inclusion distributed in the melt is acted upon by a number of forces.The motion of inclusions is written by the Newton's second law [27]: Metals 2018, 8, 392 5 of 20 The terms on the right side of Equation ( 16) represent in the following order: drag force, buoyancy force, gravity force, lift force, virtual mass force, electromagnetic pressure force, pressure gradient force, Magnus force, and Marangoni force.
The drag force, → F Drag , imposed on the particle by the melt makes the particle follow the fluid flow, and it can be calculated by [28]: where C D is drag force coefficient.d p is particle diameter.( ) 3   1 6 The Saffman's lift force, Lift F  , which is induced by the local velocity gradient around the particle, can be deduced by [30]: The virtual mass force, Vm F  , was caused by relative acceleration between the particle and the surrounding fluid [31]: where Vm C is virtual mass force model constant.
The electromagnetic pressure force, Emf F  , imposed on the particle is generated due to the pressure gradient caused by the Lorentz force in the melt [32].The inclusions are assumed to be electric insulated, and this force thus is written, as follows: The pressure gradient force, P F  , is important for the low-density particles within a large velocity gradient of the melt: The Magnus force, M F  , caused by a pressure difference will 'suck' the particle and curve its trajectory, and it is defined as [33]: where ω  is relative fluid-particle angular velocity.
The Marangoni force, Ma F  , is related to the gradient of surface tension [34]: The buoyancy force, → F B and gravity force, → F G act upward and downward along the gravitational acceleration, respectively.The two forces are related to the difference between the melt and the particle densities [29]: The Saffman's lift force, → F Li f t , which is induced by the local velocity gradient around the particle, can be deduced by [30]: The virtual mass force, → F Vm , was caused by relative acceleration between the particle and the surrounding fluid [31]: where C Vm is virtual mass force model constant.
The electromagnetic pressure force, → F Em f , imposed on the particle is generated due to the pressure gradient caused by the Lorentz force in the melt [32].The inclusions are assumed to be electric insulated, and this force thus is written, as follows: The pressure gradient force, → F P , is important for the low-density particles within a large velocity gradient of the melt: Metals 2018, 8, 392 6 of 20 The Magnus force, → F M , caused by a pressure difference will 'suck' the particle and curve its trajectory, and it is defined as [33]: where → ω is relative fluid-particle angular velocity.The Marangoni force, → F Ma , is related to the gradient of surface tension [34]: The random walk model is adopted to include the chaotic effect of the turbulence on the inclusion trajectories in the melt.The instantaneous velocity is composed of a mean velocity and a random value that is proportional to the local turbulent kinetic energy [27]: where ς is a Gaussian-distributed random number with the mean of zero and standard deviation of unity.
During the solidification, the segregated elements that were rejected by the solidifying dendrites react among them and form new various inclusions, such as oxides and sulphides, in the mushy zone.The entrapment of inclusions is significantly affected by the dendrite arm spacing and the fluid flow condition in the mushy zone, as shown in Figure 3.It is firstly assumed that the particles are entirely wrapped by the liquid.The inclusions near the tip of primary dendrite can move freely.With the decrease of the liquid fraction and the growth of secondary dendrite, the fluidity of interdendritic liquid becomes lower.According to the study of Yamazaki et al. [35], the inclusions would be entrapped if the local liquid fraction is less than 0.6.
where ς is a Gaussian-distributed random number with the mean of zero and standard deviation of unity.
During the solidification, the segregated elements that were rejected by the solidifying dendrites react among them and form new various inclusions, such as oxides and sulphides, in the mushy zone.The entrapment of inclusions is significantly affected by the dendrite arm spacing and the fluid flow condition in the mushy zone, as shown in Figure 3.It is firstly assumed that the particles are entirely wrapped by the liquid.The inclusions near the tip of primary dendrite can move freely.With the decrease of the liquid fraction and the growth of secondary dendrite, the fluidity of interdendritic liquid becomes lower.According to the study of Yamazaki et al. [35], the inclusions would be entrapped if the local liquid fraction is less than 0.6.

Boudnary Conditions
An axisymmetric physical geometry is established to study the inclusion behavior, especially the radial distribution of inclusions.The boundary conditions of the calculated domain are shown in Figure 4.A potential gradient is imposed at the electrode tip interface, accounting for the skin effect [36].Since the leakage current is neglected, the potential gradient at the mold lateral wall is supposed to be zero.Besides, the electrode melt rate is related to the Joule heating near the electrode tip [14]:

Boudnary Conditions
An axisymmetric physical geometry is established to study the inclusion behavior, especially the radial distribution of inclusions.The boundary conditions of the calculated domain are shown in Figure 4.A potential gradient is imposed at the electrode tip interface, accounting for the skin effect [36].Since the leakage current is neglected, the potential gradient at the mold lateral wall is supposed to be zero.Besides, the electrode melt rate is related to the Joule heating near the electrode tip [14]: where η represents the utilization rate of Joule heating for the melting of electrode, and is adjusted according to experiments and literature [37,38].
The mold wall and the bottom are supposed to be no-slip boundary, while the top surface is set to be zero shear stress.A velocity inlet boundary with a 30 K superheat than the liquidus of steel is enforced at the electrode tip.A suppositional pull velocity that is coherent with electrode melt rate is patched on the whole domain so that the ingot growth can be considered.Because of the air gap and the solidified slag skin, it is hard to determine the heat transfer at mold wall.In this study, only heat convection is taken into account, and an equivalent heat transfer coefficient with a value of 300~400 W/m 2 •K is used [39,40].The detailed parameters are listed in Table 1 [13,40].
Only few inclusions from electrode enter into the liquid metal pool, and if this happens, then they are subject to an inclusion modification for a simple model in the current work [4,5].Due to the interaction and rejection of segregated elements, the new inclusions are assumed to nucleate uniformly from 0.9 to 0.99 liquid fraction.The diameters of inclusions are 2, 5, and 10 µm, respectively.The inclusions will be removed from the calculated domain when they pass through the slag/metal interface.Moreover, the inclusions would be entrapped by the solidified slag skin when they reach the mold wall, and escape from the calculated domain once touching the bottom.uniformly from 0.9 to 0.99 liquid fraction.The diameters of inclusions are 2, 5, and 10 μm, respectively.The inclusions will be removed from the calculated domain when they pass through the slag/metal interface.Moreover, the inclusions would be entrapped by the solidified slag skin when they reach the mold wall, and escape from the calculated domain once touching the bottom.

Calculation Procedure
The CFD software FLUENT (Version 12.0, ANSYS Inc., Canonsburg, PA, USA) with user defined functions (UDFs) was employed to calculate the coupled transient axisymmetric model in the current work.Complicated phenomena involving electromagnetism, fluid flow, heat transfer, phase change, and inclusion behaviors were solved numerically and simultaneously.A mesh independence study was done.Three cases with structured grid sizes of 0.5 mm, 1 mm, and 2 mm were simulated, respectively.The averaged deviation of the temperature distribution along the half radius between the former two cases is about 4%, but 9% between the latter two cases.Finally, the axisymmetric domain is divided into structured grids with the size of 1 mm, as shown in Figure 4 when considering the reasonable computational time and adequate accuracy.The step size was reduced to be 0.005 s to fulfill the convergence criteria (10 −6 ).The transient calculation of one case needs approximately 230 h using a multi-CPU cluster (8 cores, 4.0 GHz, parallel computing).

Experiment Details and Model Validation
An experiment was carried out using a 200 kg scale ESR furnace under argon protection, as displayed in Figure 5.The inner diameter and height of mold are 200 mm and 1150 mm, respectively.The electrode was AISI (American Iron and Steel Institute) 304 stainless steel with a diameter of 100 mm.The original premelting slag was (wt %): CaF 2 70.0, Al 2 O 3 30.0.The dried argon gas was blown into the airtight cover with a flowrate of 30 Nl/min.The current was 3100 A during the steady melting stage.The whole ESR process, including the arcing and feeding stages lasted about 130 min, and neither deoxidizer nor desulphurizer was added.
displayed in Figure 5.The inner diameter and height of mold are 200 mm and 1150 mm, respectively.The electrode was AISI (American Iron and Steel Institute) 304 stainless steel with a diameter of 100 mm.The original premelting slag was (wt %): CaF2 70.0, Al2O3 30.0.The dried argon gas was blown into the airtight cover with a flowrate of 30 Nl/min.The current was 3100 A during the steady melting stage.The whole ESR process, including the arcing and feeding stages lasted about 130 min, and neither deoxidizer nor desulphurizer was added.After remelting, a longitudinal section with a thickness of 11 mm was obtained by wire-electrode cutting along the axis.The section was firstly etched by mixture acid (75% HCl, 25% HNO3) for the macrostructure.Twelve samples (10 mm × 10 mm × 10 mm) were taken from the section along radial direction at different heights.These samples were carefully ground, polished, and analysed by the SEM/EDS (JSM-7800F, JEOL Ltd., Tokyo, Japan) for the morphology and composition of the inclusions.The amount and size of nonmetallic inclusions were counted by Image J software (Version 1.48, National Institutes of Health (NIH), Bethesda, MD, USA) through 40 random SEM fields of view of each sample.The total oxygen content was determined by oxygen/nitrogen/hydrogen analyser (ONH836, LECO, St. Joseph, MI, USA).
Figure 6 displayed the comparison of experimental and simulated metal pool profiles.The growth of columnar crystals are dominant in final ingot, while no equiaxed crystal is observed.The depth of metal pool increases with the sequent solidification due to a decreasing cooling effect of the After remelting, a longitudinal section with a thickness of 11 mm was obtained by wire-electrode cutting along the axis.The section was firstly etched by mixture acid (75% HCl, 25% HNO 3 ) for the macrostructure.Twelve samples (10 mm × 10 mm × 10 mm) were taken from the section along radial direction at different heights.These samples were carefully ground, polished, and analysed by the SEM/EDS (JSM-7800F, JEOL Ltd., Tokyo, Japan) for the morphology and composition of the inclusions.The amount and size of nonmetallic inclusions were counted by Image J software (Version 1.48, National Institutes of Health (NIH), Bethesda, MD, USA) through 40 random SEM fields of view of each sample.The total oxygen content was determined by oxygen/nitrogen/hydrogen analyser (ONH836, LECO, St. Joseph, MI, USA).
Figure 6 displayed the comparison of experimental and simulated metal pool profiles.The growth of columnar crystals are dominant in final ingot, while no equiaxed crystal is observed.The depth of metal pool increases with the sequent solidification due to a decreasing cooling effect of the water tank at bottom.The depths are 37 mm, 43 mm, and 45 mm at the bottom, middle, and upper part of ingot, respectively.The simulated quasi-steady metal pool profile, of which depth is 44 mm, agrees well with the upper and middle ones.The discrepancy probably results from the uncertainty of material properties and boundary condition, the electrical insulation hypothesis of the mold, and the measuring errors.
Figure 7 shows the morphology and composition of typical inclusions in samples of electrode and ESR ingot.The inclusions are mainly MnS, TiN, and (Al, Ti, Mn) O-MnS with irregular structure in electrode, as illustrated in Figure 7a-c.After the ESR process, the inclusions are Al 2 O 3 core surrounded by an outer sulfide layer and pure MnS inclusions, as shown in Figure 7d-h.Most of the oxide inclusions are near-spherical with a smaller size than the ones in electrode.The large amount of existence of sulfide inclusions results from the absence of desulphurization treatment.Few inclusions including TiN, as shown in Figure 7i, is observed.
water tank at bottom.The depths are 37 mm, 43 mm, and 45 mm at the bottom, middle, and upper part of ingot, respectively.The simulated quasi-steady metal pool profile, of which depth is 44 mm, agrees well with the upper and middle ones.The discrepancy probably results from the uncertainty of material properties and boundary condition, the electrical insulation hypothesis of the mold, and the measuring errors.The size of almost all the nonmetallic inclusions is 1-5 µm, and no inclusion larger than 10 µm is observed.The evolutions of total area of nonmetallic inclusions along radial direction at different heights of ESR ingot are shown in Figure 8.The experimental results indicate that the total area of inclusions increases from center toward ingot surface, while the average equivalent diameter decreases over radial direction.Meanwhile, the total area increases with the ingot growth due to the accumulation of inclusions in a liquid metal pool.The total oxygen content, which is closely related to oxide inclusions, also shows an upward trend with height.
Given that the total number of inclusions in simulation differs from that in experiment, the distribution ratio along radius direction thus is used to validate the model.The ratios at three different radiuses are assessed, as demonstrated in Figure 9.The good agreement between the calculation and experiment indicates that the model is able to give a further prediction of inclusion behaviour in ESR process.
The size of almost all the nonmetallic inclusions is 1-5 μm, and no inclusion larger than 10 μm is observed.The evolutions of total area of nonmetallic inclusions along radial direction at different heights of ESR ingot are shown in Figure 8.The experimental results indicate that the total area of inclusions increases from center toward ingot surface, while the average equivalent diameter decreases over radial direction.Meanwhile, the total area increases with the ingot growth due to the accumulation of inclusions in a liquid metal pool.The total oxygen content, which is closely related to oxide inclusions, also shows an upward trend with height.
Given that the total number of inclusions in simulation differs from that in experiment, the distribution ratio along radius direction thus is used to validate the model.The ratios at three different radiuses are assessed, as demonstrated in Figure 9.The good agreement between the calculation and experiment indicates that the model is able to give a further prediction of inclusion behaviour in ESR process.The size of almost all the nonmetallic inclusions is 1-5 μm, and no inclusion larger than 10 μm is observed.The evolutions of total area of nonmetallic inclusions along radial direction at different heights of ESR ingot are shown in Figure 8.The experimental results indicate that the total area of inclusions increases from center toward ingot surface, while the average equivalent diameter decreases over radial direction.Meanwhile, the total area increases with the ingot growth due to the accumulation of inclusions in a liquid metal pool.The total oxygen content, which is closely related to oxide inclusions, also shows an upward trend with height.
Given that the total number of inclusions in simulation differs from that in experiment, the distribution ratio along radius direction thus is used to validate the model.The ratios at three different radiuses are assessed, as demonstrated in Figure 9.The good agreement between the calculation and experiment indicates that the model is able to give a further prediction of inclusion behaviour in ESR process.

MHD (Magnetohydrodynamics) Fluid Flow and Temperature Distribution
The Joule heating is the unique heat source to melt the electrode and to maintain the slag temperature in ESR process.Figure 10 demonstrates the predicted distributions of current density and Joule heating at 2425.15 s.The molten steel that is driven by gravity and electromagnetic force converges to the center of electrode tip, and then forms molten droplets periodically.The current tends to pass through the molten droplet that has a higher conductivity once flowing into the calculated domain.A higher current density is thus observed in the droplet, as well as the periphery of electrode tip due to the skin effect.The Joule heating, which is closely related to the current density and electrical conductivity, is constantly changing with the phase distribution, and is mainly produced in the slag layer.The maximum Joule heating locates under the periphery of electrode tip due to the skin effect.Few Joule heating is generated in the outer layer of slag due to the nonconductive mold.Additionally, a higher region of Joule heating is also observed in the slag around the droplet tip with the growth of droplet.
temperature in ESR process.Figure 10 demonstrates the predicted distributions of current density and Joule heating at 2425.15 s.The molten steel that is driven by gravity and electromagnetic force converges to the center of electrode tip, and then forms molten droplets periodically.The current tends to pass through the molten droplet that has a higher conductivity once flowing into the calculated domain.A higher current density is thus observed in the droplet, as well as the periphery of electrode tip due to the skin effect.The Joule heating, which is closely related to the current density and electrical conductivity, is constantly changing with the phase distribution, and is mainly produced in the slag layer.The maximum Joule heating locates under the periphery of electrode tip due to the skin effect.Few Joule heating is generated in the outer layer of slag due to the nonconductive mold.Additionally, a higher region of Joule heating is also observed in the slag around the droplet tip with the growth of droplet.The distribution of temperature, as well as the flow pattern, has been displayed in Figure 11.A hotter region locates in upper part of slag layer.The maximum temperature, which can be seen under the periphery of electrode tip, has reached 2021.48K.The strong cooling effect near the mold leads to a great temperature gradient.Moreover, the droplets also transfer heat from slag to liquid metal pool.The temperature distribution differs from the Joule heating distribution due to the fluid flow.Two vortices are clearly observed in the slag layer.The strong cooling effect of mold results in a descending flow near the mold, while the falling droplets and the electromagnetic force dominate the flow in the center.The falling droplets also cause fluctuations of slag/metal pool interface.On the other hand, the thermal buoyancy force plays a major role in the liquid metal pool, and a bigger vortex can be seen near the mold.The sheer stress that is imposed on the molten metal by slag also contributes to a local shear flow near the interface.
The flow pattern would significantly affect the motion of nonmetallic inclusions.The flow pattern in liquid metal pool keeps a dynamic stability except when the droplet falling into the metal pool, as shown in Figure 12.During the formation period of droplet as shown in Figure 12a, a wellmarked clockwise flow is observed near the mold.The descending flow induced by the cooling effect, which has a higher velocity, reaches at the mushy zone, and then sweeps through the solidification front from outside toward center.Conversely, a counterclockwise flow affected by the Lorentz force with a slower velocity pushes the inclusions outward in the center.Once the droplet falls into the liquid metal pool, as displayed in Figure 12b, the flow pattern in the center changes with a temporary flow with the direction of falling droplet.It also should be noted that the velocity decreases with the reduction of liquid fraction.The distribution of temperature, as well as the flow pattern, has been displayed in Figure 11.A hotter region locates in upper part of slag layer.The maximum temperature, which can be seen under the periphery of electrode tip, has reached 2021.48K.The strong cooling effect near the mold leads to a great temperature gradient.Moreover, the droplets also transfer heat from slag to liquid metal pool.The temperature distribution differs from the Joule heating distribution due to the fluid flow.Two vortices are clearly observed in the slag layer.The strong cooling effect of mold results in a descending flow near the mold, while the falling droplets and the electromagnetic force dominate the flow in the center.The falling droplets also cause fluctuations of slag/metal pool interface.On the other hand, the thermal buoyancy force plays a major role in the liquid metal pool, and a bigger vortex can be seen near the mold.The sheer stress that is imposed on the molten metal by slag also contributes to a local shear flow near the interface.
The flow pattern would significantly affect the motion of nonmetallic inclusions.The flow pattern in liquid metal pool keeps a dynamic stability except when the droplet falling into the metal pool, as shown in Figure 12.During the formation period of droplet as shown in Figure 12a, a well-marked clockwise flow is observed near the mold.The descending flow induced by the cooling effect, which has a higher velocity, reaches at the mushy zone, and then sweeps through the solidification front from outside toward center.Conversely, a counterclockwise flow affected by the Lorentz force with a slower velocity pushes the inclusions outward in the center.Once the droplet falls into the liquid metal pool, as displayed in Figure 12b, the flow pattern in the center changes with a temporary flow with the direction of falling droplet.It also should be noted that the velocity decreases with the reduction of liquid fraction.

Inclusion Motion and Distribution
Figure 13 shows the motion of newly generated inclusions in liquid metal pool.The inclusions at the half radius are firstly uplifted by the flow, as displayed in Figure 13a.Subsequently, the inclusions moves with the vortexes, as presented in Figure 13b,c.It should be noted that more inclusions are entrapped into the outer vortex.Part of the inclusions rise to the slag/metal interface under the buoyancy, and a small amount enter into the slag through flotation.The descending flow near the mold, wrapping most of the inclusions, sweeps through the solidification front continuously.Finally, an aggregation is easily observed near the mold, as shown in Figure 13d.

Inclusion Motion and Distribution
Figure 13 shows the motion of newly generated inclusions in liquid metal pool.The inclusions at the half radius are firstly uplifted by the flow, as displayed in Figure 13a.Subsequently, the inclusions moves with the vortexes, as presented in Figure 13b,c.It should be noted that more inclusions are entrapped into the outer vortex.Part of the inclusions rise to the slag/metal interface under the buoyancy, and a small amount enter into the slag through flotation.The descending flow near the mold, wrapping most of the inclusions, sweeps through the solidification front continuously.Finally, an aggregation is easily observed near the mold, as shown in Figure 13d.

Inclusion Motion and Distribution
Figure 13 shows the motion of newly generated inclusions in liquid metal pool.The inclusions at the half radius are firstly uplifted by the flow, as displayed in Figure 13a.Subsequently, the inclusions moves with the vortexes, as presented in Figure 13b,c.It should be noted that more inclusions are entrapped into the outer vortex.Part of the inclusions rise to the slag/metal interface under the buoyancy, and a small amount enter into the slag through flotation.The descending flow near the mold, wrapping most of the inclusions, sweeps through the solidification front continuously.Finally, an aggregation is easily observed near the mold, as shown in Figure 13d.Figure 14 demonstrates the inclusion distribution in the metal pool and ingot.The color of particles represents the residence time of particles.Few inclusions are observed near the ingot surface due to the strong descending flow.Most of the inclusions are carried inward and precipitate with solidification sequence.A clear denser distribution is observed at the region ranging from 0.7 to 0.9 radius.Since the metal pool profile does not vary with the solidification in this quasi-steady simulation, the radial distributions of inclusions at different height of ingot are similar.The earlier generated inclusions are mainly distributed in the lower part of ingot, but a small proportion still locates in the upper part even the metal pool.This phenomenon indicates that the newly generated inclusions would not be solidified entirely in a short time, and a gradual accumulation of inclusions in metal pool occurs with the growth of ingot.In this case, a denser distribution of inclusions would be calculated in the upper part of ingot.
Figure 15 shows the predicted distribution of inclusions with different diameters.The inclusions with different diameters are easier to concentrate in the region ranging from 0.7 to 0.9 radius.A further analysis of the distribution of inclusions with different diameters along radius is conducted in Figure 16.The inclusion amount increases along radius, but slightly decreases when ranging from 0.3 to 0.4 radius because of a local rising flow.The inclusion amount reaches a maximum at 0.8 radiuses, and it rapidly falls down at the ingot surface.There are no significant differences among the inclusion amount with different particle diameters in the center.However, the inclusion amount with smaller diameter is obviously more than the one with larger diameter at the region ranging from 0.7 to 0.9 radius.The ratios of the inclusion amount in this range to the total with the diameter of 2 μm, 5 μm, and 10 μm are 48.76%,44.23%, and 42.48%, respectively.The origin of this phenomenon is that small inclusions easily follow the liquid flow.The simulated result can account for the experimental result reasonably, that is, the total area of inclusions increases along radius but the average equivalent diameter decreases.Figure 14 demonstrates the inclusion distribution in the metal pool and ingot.The color of particles represents the residence time of particles.Few inclusions are observed near the ingot surface due to the strong descending flow.Most of the inclusions are carried inward and precipitate with solidification sequence.A clear denser distribution is observed at the region ranging from 0.7 to 0.9 radius.Since the metal pool profile does not vary with the solidification in this quasi-steady simulation, the radial distributions of inclusions at different height of ingot are similar.The earlier generated inclusions are mainly distributed in the lower part of ingot, but a small proportion still locates in the upper part even the metal pool.This phenomenon indicates that the newly generated inclusions would not be solidified entirely in a short time, and a gradual accumulation of inclusions in metal pool occurs with the growth of ingot.In this case, a denser distribution of inclusions would be calculated in the upper part of ingot.
Figure 15 shows the predicted distribution of inclusions with different diameters.The inclusions with different diameters are easier to concentrate in the region ranging from 0.7 to 0.9 radius.A further analysis of the distribution of inclusions with different diameters along radius is conducted in Figure 16.The inclusion amount increases along radius, but slightly decreases when ranging from 0.3 to 0.4 radius because of a local rising flow.The inclusion amount reaches a maximum at 0.8 radiuses, and it rapidly falls down at the ingot surface.There are no significant differences among the inclusion amount with different particle diameters in the center.However, the inclusion amount with smaller diameter is obviously more than the one with larger diameter at the region ranging from 0.7 to 0.9 radius.The ratios of the inclusion amount in this range to the total with the diameter of 2 µm, 5 µm, and 10 µm are 48.76%,44.23%, and 42.48%, respectively.The origin of this phenomenon is that small inclusions easily follow the liquid flow.The simulated result can account for the experimental result reasonably, that is, the total area of inclusions increases along radius but the average equivalent diameter decreases.Figure 17 displayed the effect of particle diameter on the inclusion distribution along height in ingot.The solidified domain (liquid fraction of metal less than 0.6) is evenly divided into three parts along height.It can be found that the inclusion amount increases with height.The ratios of inclusions with the diameter of 5 µm to total inclusion amount in the lower part, middle part, and upper part are 32.19%,32.21%, and 32.78%, respectively.Similarly, the ratios of inclusions with the diameter of 10 µm to total inclusion amount in the lower part, middle part, and upper part are 29.27%,29.58%, and 29.67%, respectively.This increasing trend results from that the large inclusions are imposed by a greater buoyancy force.A bigger average equivalent diameter of inclusions in the upper part is therefore deduced.Moreover, the total amounts of inclusion with the diameter of 2 µm, 5 µm, and 10 µm in ingot are 2199, 1871, and 1704, respectively.This is because the large inclusions are easier to enter into the slag by floatation.Figure 17 displayed the effect of particle diameter on the inclusion distribution along height in ingot.The solidified domain (liquid fraction of metal less than 0.6) is evenly divided into three parts along height.It can be found that the inclusion amount increases with height.The ratios of inclusions with the diameter of 5 μm to total inclusion amount in the lower part, middle part, and upper part are 32.19%,32.21%, and 32.78%, respectively.Similarly, the ratios of inclusions with the diameter of 10 μm to total inclusion amount in the lower part, middle part, and upper part are 29.27%,29.58%, and 29.67%, respectively.This increasing trend results from that the large inclusions are imposed by a greater buoyancy force.A bigger average equivalent diameter of inclusions in the upper part is therefore deduced.Moreover, the total amounts of inclusion with the diameter of 2 μm, 5 μm, and 10 μm in ingot are 2199, 1871, and 1704, respectively.This is because the large inclusions are easier to enter into the slag by floatation.

Conclusions
A transient 2D axisymmetric model has been developed to predict the motion and distribution of nonmetallic inclusions in ESR process.The Euler-Lagrange approach was adopted to describe the interaction between the continuous phase and nonmetallic inclusions.The VOF approach was employed to track the two-phase flow, while the DPM model was used to track the inclusions.The simulated result agrees well with the experimental one with an acceptable discrepancy.
1.The experiment shows that the growths of columnar crystals are dominant in the ESR ingot.The typical inclusions in ESR ingot are Al2O3 core surrounded by an outer sulfide layer and pure MnS inclusions.Most of the oxide inclusions are near-spherical with a smaller size than that in electrode.The total area of inclusions increases from center toward ingot surface, while the average equivalent diameter decreases.2. The thermal buoyancy force plays a major role in the metal pool, and it drives a descending flow near the mold.The descending flow sweeps through the solidification front inward carrying the inclusions.A reverse smaller vortex is observed in the center.A rising flow is hence produced at 0.5 radius.The velocity decreases with the reduction of liquid fraction in the mushy zone.3. The simulated result overall shows that the inclusion amount increases along radius.A denser distribution of inclusions is observed at the region ranging from 0.7 to 0.9 radius.A slight decrease when ranging from 0.3 to 0.4 radius is observed because of a local rising flow.Few inclusions are entrapped near the mold due to the strong descending flow.The earlier generated inclusions are mainly the distribution in the lower part of ingot, but a small proportion is still located in the upper part.

Conclusions
A transient 2D axisymmetric model has been developed to predict the motion and distribution of nonmetallic inclusions in ESR process.The Euler-Lagrange approach was adopted to describe the interaction between the continuous phase and nonmetallic inclusions.The VOF approach was employed to track the two-phase flow, while the DPM model was used to track the inclusions.The simulated result agrees well with the experimental one with an acceptable discrepancy.

1.
The experiment shows that the growths of columnar crystals are dominant in the ESR ingot.The typical inclusions in ESR ingot are Al 2 O 3 core surrounded by an outer sulfide layer and pure MnS inclusions.Most of the oxide inclusions are near-spherical with a smaller size than that in electrode.The total area of inclusions increases from center toward ingot surface, while the average equivalent diameter decreases.

2.
The thermal buoyancy force plays a major role in the metal pool, and it drives a descending flow near the mold.The descending flow sweeps through the solidification front inward carrying the inclusions.A reverse smaller vortex is observed in the center.A rising flow is hence produced at 0.5 radius.The velocity decreases with the reduction of liquid fraction in the mushy zone.

3.
The simulated result overall shows that the inclusion amount increases along radius.A denser distribution of inclusions is observed at the region ranging from 0.7 to 0.9 radius.A slight decrease when ranging from 0.3 to 0.4 radius is observed because of a local rising flow.Few inclusions are entrapped near the mold due to the strong descending flow.The earlier generated inclusions are mainly the distribution in the lower part of ingot, but a small proportion is still located in the upper part.4.
There are no significant difference among the inclusion amount with different particle diameters in the center of ingot, but the smaller inclusions are clearly more than the larger ones at the region from 0.7 to 0.9 radius.

5.
The inclusion amount increases with height.The total amounts of inclusions with diameter of 2 µm, 5 µm and 10 µm in ingot are 2199, 1871, and 1704, respectively.The imposed buoyancy and floatation contribute to this phenomenon.

Figure 2 .
Figure 2. Force balance on a typical particle in the liquid metal pool.

Figure 2 .
Figure 2. Force balance on a typical particle in the liquid metal pool.

Figure 3 .
Figure 3. Schematic illustration of inclusion behavior during the solidification.

Figure 3 .
Figure 3. Schematic illustration of inclusion behavior during the solidification.

Figure 5 .
Figure 5. Photos of the experiment: (a) 200 kg argon protection ESR furnace; (b) ESR ingot; (c) Illustration of wire cutting on ESR ingot, the longitudinal section labelled by gray color is used for the current study.

Figure 5 .
Figure 5. Photos of the experiment: (a) 200 kg argon protection ESR furnace; (b) ESR ingot; (c) Illustration of wire cutting on ESR ingot, the longitudinal section labelled by gray color is used for the current study.

Figure 6 .
Figure 6.Comparison of metal pool profiles between experiment and simulation.

Figure 7
Figure 7 shows the morphology and composition of typical inclusions in samples of electrode and ESR ingot.The inclusions are mainly MnS, TiN, and (Al, Ti, Mn) O-MnS with irregular structure in electrode, as illustrated in Figure 7a-c.After the ESR process, the inclusions are Al2O3 core surrounded by an outer sulfide layer and pure MnS inclusions, as shown in Figure 7d-h.Most of the oxide inclusions are near-spherical with a smaller size than the ones in electrode.The large amount of existence of sulfide inclusions results from the absence of desulphurization treatment.Few inclusions including TiN, as shown in Figure 7i, is observed.

Figure 6 .
Figure 6.Comparison of metal pool profiles between experiment and simulation.

Figure 6 .
Figure 6.Comparison of metal pool profiles between experiment and simulation.

Figure 7
Figure 7 shows the morphology and composition of typical inclusions in samples of electrode and ESR ingot.The inclusions are mainly MnS, TiN, and (Al, Ti, Mn) O-MnS with irregular structure in electrode, as illustrated in Figure 7a-c.After the ESR process, the inclusions are Al2O3 core surrounded by an outer sulfide layer and pure MnS inclusions, as shown in Figure 7d-h.Most of the oxide inclusions are near-spherical with a smaller size than the ones in electrode.The large amount of existence of sulfide inclusions results from the absence of desulphurization treatment.Few inclusions including TiN, as shown in Figure 7i, is observed.

Figure 7 .
Figure 7. Morphology and composition of typical inclusions in samples of electrode and ESR ingot: (a-c) electrode and (d-i) ESR ingot.

Figure 8 .
Figure 8. Evolution of total are of nonmetallic inclusions along radial direction at different heights of ESR ingot.Only 1-5 μm inclusions are counted.

Figure 9 .
Figure 9.Comparison of distribution ratio of inclusions along radius direction between simulation and experiment.

Figure 8 .
Figure 8. Evolution of total are of nonmetallic inclusions along radial direction at different heights of ESR ingot.Only 1-5 µm inclusions are counted.

Figure 8 .
Figure 8. Evolution of total are of nonmetallic inclusions along radial direction at different heights of ESR ingot.Only 1-5 μm inclusions are counted.

Figure 9 .
Figure 9.Comparison of distribution ratio of inclusions along radius direction between simulation and experiment.Figure 9. Comparison of distribution ratio of inclusions along radius direction between simulation and experiment.

Figure 9 .
Figure 9.Comparison of distribution ratio of inclusions along radius direction between simulation and experiment.Figure 9. Comparison of distribution ratio of inclusions along radius direction between simulation and experiment.

Figure 10 .
Figure 10.Predicted distributions of current density (a) and Joule heating (b) at 2425.15 s.The red lines represent the slag/metal interface.The yellow lines represent the current streamlines.

Figure 10 .
Figure 10.Predicted distributions of current density (a) and Joule heating (b) at 2425.15 s.The red lines represent the slag/metal interface.The yellow lines represent the current streamlines.

Figure 11 .
Figure 11.Predicted distribution of temperature at 2485.65 s.The red lines represent the slag/metal interface, and the black ones are streamlines.

Figure 12 .
Figure 12.Typical map of velocity vector in liquid metal pool at 2485.65 s (a) and 2621.56 (b).The solidified ingot and mushy zone are marked in gray color.

Figure 11 . 19 Figure 11 .
Figure 11.Predicted distribution of temperature at 2485.65 s.The red lines represent the slag/metal interface, and the black ones are streamlines.

Figure 12 .
Figure 12.Typical map of velocity vector in liquid metal pool at 2485.65 s (a) and 2621.56 (b).The solidified ingot and mushy zone are marked in gray color.

Figure 12 .
Figure 12.Typical map of velocity vector in liquid metal pool at 2485.65 s (a) and 2621.56 (b).The solidified ingot and mushy zone are marked in gray color.

Figure 14 .
Figure 14.Predicted distribution of inclusions in the metal pool and ingot.Only the inclusions with diameter of 2 μm are displayed.

Figure 15 .
Figure 15.Predicted distribution of inclusions with different diameters.Only one-sixth of the total inclusions are displayed.

Figure 14 . 20 Figure 14 .
Figure 14.Predicted distribution of inclusions in the metal pool and ingot.Only the inclusions with diameter of 2 µm are displayed.

Figure 15 .
Figure 15.Predicted distribution of inclusions with different diameters.Only one-sixth of the total inclusions are displayed.

Figure 15 .
Figure 15.Predicted distribution of inclusions with different diameters.Only one-sixth of the total inclusions are displayed.

Figure 15 .
Figure 15.Predicted distribution of inclusions with different diameters.Only one-sixth of the total inclusions are displayed.

Figure 16 .
Figure 16.Effect of particle diameter on the inclusion distribution along radius in ingot.Figure 16.Effect of particle diameter on the inclusion distribution along radius in ingot.

Figure 16 .
Figure 16.Effect of particle diameter on the inclusion distribution along radius in ingot.Figure 16.Effect of particle diameter on the inclusion distribution along radius in ingot.

Figure 17 .
Figure 17.Effect of particle diameter on the inclusion distribution along height in ingot.

Figure 17 .
Figure 17.Effect of particle diameter on the inclusion distribution along height in ingot.

F
damp damping force in the mushy zone (N/m 3 ) Marangoni force (N/m 3 ) → F P pressure gradient force (N/m 3 ) → F st surface tension (N/m 3 ) → F Vm virtual mass force (N/m 3 ) particle angular velocity (rad/s) is the effective viscosity composed of the molecular and turbulent viscosities, i.e., µ e f f = µ + µ t .The RNG (Re-Normalized Group) k − ε turbulence model that performs a wider application of turbulent number is used to calculate the turbulent one, i.e., µ t = ρC µ

Table 1 .
Physical properties and operation conditions.

Table 1 .
Physical properties and operation conditions.