Modelling the Nucleation, Growth and Agglomeration of Alumina Inclusions in Molten Steel by Combining Kampmann–Wagner Numerical Model with Particle Size Grouping Method

Recent inclusion models are mainly focused on the compositional evolution of inclusion, steel and slag. Due to the importance of inclusion size distribution to steel properties, the evolution of inclusion size distributions should also be accounted for. As the first step to establish a model to predict the evolution of inclusion size distribution, the nucleation, growth and removal of alumina inclusions in molten steel were modeled by combining Kampmann and Wagner numerical model for nucleation, growth and coarsening with particle size grouping method. The model could simulate the time evolution of the size distribution of alumina inclusions after aluminum de-oxidation. The model was validated by using the experimental size distribution data of alumina inclusions available in the literature. The model calculation results were also compared with previous simulation results. The influences of interfacial tension between steel and inclusion and diffusion coefficient on the calculated inclusion size distribution were investigated. As interfacial tension between steel and alumina increases, the maximum number density decreases and the peak value of radius increases. As diffusion coefficient increases, the maximum number density decreases and the peak-value radius increases. The calculated size distribution curves showed a change from log normal to fractal, which is due to the change of dominating mechanisms for crystal growth and agglomeration.


I. INTRODUCTION
DEOXIDATION of steel during secondary metallurgy generates enormous non-metallic inclusions which are detrimental to mechanic properties of steels. [1] Therefore, the inclusions should be removed to the top slag layer by the flotation from ladle treatment to tundish. The size distributions of inclusions are essential to the properties of steel due to the fact that the inclusions with large sizes are particularly harmful to properties (e.g., fatigue) of final products. [2] In steel refining practice, inclusions with larger size could be removed more easily due to the larger terminal velocity of the flotation.
There were some modelling work to calculate the evolution of size distribution during nucleation, growth and agglomeration of inclusions. [3][4][5][6][7][8] Nakaoka et al. [3] established a so-called particle size grouping (PSG) method to enable a simple calculation of the agglomeration by a small number of size groups with complete conservation in total particle volume. Zhang and Pluschkell [4] proposed a numerical model to simulate the nucleation, Ostwald ripening and agglomeration due to Brownian and turbulent collisions of inclusions. They employed the method by Kampmann and Kahlweit [5] to simulate the nucleation and Ostwald ripening process of inclusions. An improved model was later put forward by Zhang and Lee. [6] They divided the total size spectrum into discrete regime and sectional regime to solve the population balance model (PBM) efficiently. Later, this model was combined with a Computational Fluid Dynamics (CFD) model to predict the size distribution of alumina at any position in ladle. [7] Lei et al. [8] also developed a similar numerical model to simulate the nucleation and growth of inclusions. It should be mentioned that Zhang-Pluschkell, Zhang-Lee and Lei's model all employed the method by Kampmann and Kahlweit to simulate the nucleation and Ostwald ripening of inclusions. It was found that particle size grouping method could not easily be combined with Kampmann and Kahlweit method since monomer addition in particular group can never support enough particle growth in order for it to be listed in the neighboring larger size group. [6,7] Later than Kampmann and Kahlweit method, Kampmann and Wagner [9] treated nucleation, growth and coarsening as couple processes by using a framework suitable for the numerical calculation based on theory by Langer and Schwartz. [10] Kampmann and Wagner numerical model (KWN) has been extensively applied in the phase precipitation during heating. [11][12][13][14] Recently Spring [15] attempted to employ KWN to model the inclusion size distribution in Al deoxidation, however, the algorithm to treat the agglomeration in their work is not clear. Besides, the model was not validated and agreements between their calculation results and previous simulation results were not very good.
There are increasing efforts in recent years to model the inclusion evolution during ladle treatment. [16][17][18][19] However, these models are mainly focused on the evolution of compositions of inclusion, steel and slag. Due to the importance of inclusion size distribution, the evolutions of inclusion size distributions should be also included in the future model.
In this work, a new method combining Kampmann and Wagner numerical (KWN) model with Particle Size Grouping (PSG) method to simulate the nucleation, growth and agglomeration of alumina inclusions in molten steel was proposed. It is shown that the KWN model can be intimately compatible with the PSG method, and there is no need to separate the total size spectrum into the discrete regime and the sectional regime.

A. Thermodynamics of Nucleation
The de-oxidation reaction by Al in molten steel could be represented as follows: The standard Gibbs-free energy for reaction [1] DG could be related to the equilibrium constant K where R is gas constant, T is temperature in Kelvin. a Al ½ and a O ½ are activities of Al and O in steel, respectively.
Since liquid steel is normally dilute solution, it is convenient to use 1 weight-percentage solute obeying Henrian law as standard state and calculate activity of solute as follows: where e j i is a first-order mass-based interaction parameter. The interaction parameters [20] used in this study are shown in Table 1.
According to classic nucleation theory, the formation of nuclei should overcome the barrier due to bulk Gibbs free energy change of de-oxidation reaction and the interfacial energy between nuclei and its parent phase. Then the total Gibbs energy change for the formation of spherical nuclei could be expressed as the following equation: where r is the interfacial tension between nuclei and the parent phase; r is the radius of nuclei;DG v is the bulk Gibbs free energy for de-oxidation reaction per unit volume.
where P is supersaturation degree represented by the following equation: V m is the molar volume of nucleation phase. The nuclei start to form when @ DG v ð Þ @r ¼ 0, then we get the critical size for nucleation B. Kampmann and Wagner Numerical Model

Nucleation
It was assumed that the nucleation of inclusions in molten steel proceeds in a way of homogenous nucleation. This is reasonable considering large supersaturation degree after additions of alloys. Then, the steady state nucleation rate I could be expressed as where I A is a pre-exponential factor, I A can be calculated by I A ¼ N A V m i cr , N A is Avogadro number, V m is molar volume, i cr is the critical number of molecules to form a nucleus which can be obtained as hom is the energy barrier for homogeneous nucleation and can be calculated as: where r is the interfacial tension between inclusion and residue steel liquid; DG V is Gibbs free energy change per molar volume for inclusion formation.

Crystal growth or dissolution
Considering a spherical inclusion in supersaturated steel melts, the concentration at the inclusion-steel interface C i could be calculated by Gibbs-Thomson equation: where Ce is the equilibrium concentration; r is the interfacial tension between inclusion and residue steel liquid; V m is the molar volume of a inclusion. It can be seen directly from the equation that the inclusion with smaller size would have a higher concentration at the interface. Depending on whether the concentration at the inclusion-steel interface is larger than the concentration in residue liquid steel, the crystal will undergo growth or dissolution. The rate of crystal growth or dissolution could be expressed as the following equation.
where C is the concentration in residue liquid steel, C p is the concentration in inclusion, C i is the concentration at the interface, D is the diffusion coefficient of oxygen in liquid steel, r is the radius of the inclusion.
The critical radius of a particle where growth or dissolution rate is zero could be obtained by combining Eq. [11] with Eq. [12]:

Numerical method
To numerically simulate the nucleation, growth and coarsening process, the particle size distribution could be discretized into a series of size groups with widths of Dri. The radius and volume of particles in each group increases with constant ratios of R r and R v respectively (R v = R r 3 ). A finite volume method (FVM) which is widely applied in computational fluid dynamics could be employed to establish the discretized equation for the evolution of size distribution. [11] A grid configuration compatible to FVM is shown in Figure 1. Following ''upwind scheme'', discretization equations could be obtained as follows [11] : where N P is the number density at P point after one iteration; N 0 P , N 0 W and N 0 E are the current number density at P, W and E grids before iteration, respectively; Dt is the time step; Dr W , Dr P and Dr E are the radius width for W, P and E group, respectively. v W and v E are rate of growth or dissolution for grid W and E.
Based on the conservation of mass, we can calculate the matrix solute concentration as follows: where C p is solute concentration in inclusions in residue liquid steel; C 0 is the initial solute concentration in residue liquid steel, q st and q P are density of liquid steel and inclusion, respectively. C. Particle Size Grouping (PSG) Method The agglomeration of particles due to collision could be well described by population balance equation derived by Smoluchowski [21] : where n k is the number density of k-mer. collision frequency function b C ij comes from the contributions of Brownian collisions b B ij , stokes collisions b S ij and turbu- where l is the dynamic viscosity of steel, k B is Boltzmann constant, e is turbulent energy dissipation rate, g is the acceleration of gravity, q st and q st is the density of steel and alumina, respectively, r i and r j is the radius of particles i and j, respectively. a t is turbulent coagulation coefficient and can be calculate by the following equation: where A as is the Hamaker constant of alumina in steel and here we use 2.3 9 10 À20 J, a 1 is the radius of primary particle. Nakaoka et al. [3] presented a particle-size-grouping method to efficiently solve the population balance equation. The particles are divided into M groups. The volume of particles of each group increases from one to the next group with a constant ratio: R v . The procedure can be described by the following equation: where d ik is the Dirac delta function; i c;k is the critical size number which was determined by the R v ratio, when R v = 2, i c;k ¼ k À 1. n i;kÀ1 and n i;k are the correction factor of the particle number density which are given by: In the present work, the particles were divided into a series of groups with volume ratio of 2 (radius ratio of 1.26). The particle size grouping method could be easily incorporated into the model by directly using Eq. [26].

D. The Removal of Inclusions
The previous studies on model of nucleation and growth of inclusions did not properly account for the removal of inclusions. In model by Zhang and Pluschkell, [4] inclusions larger than 36 lm were assumed to be removed instantly to top slag. In Zhang and Lee's work, [6] it is supposed that inclusions larger than 50 lm in diameter are immediately removed from molten steel. Lei et al. [8] assumed that the inclusions larger than 10 lm are removed from molten steel at once.
There are three mechanisms including Stokes flotation and absorbed by top slag layer, adherence to the wall and bubble-inclusion attachment for removal of inclusions proposed in literatures. [7,22,23] Due to lack of reliable parameters, only stokes flotation mechanism is considered in the present work.
According to Stokes floatation mechanism, the number density change due to the Stokes flotation in one step of simulation can be represented by the following equation: where v t is the terminal velocity ; n is the number density of group i; Dt is the time step size; H is the depth of the inclusions in steel vessel. In each iteration, Dn st was calculated for all size groups.
Considering the inclusion removal, the equation (19) should be modified as follows to keep the mass balance: The density of liquid steel can be varied with the composition of steel. For AISI 1008 (composition 0.06 pct C, 0.38 pct Mn, 0.01 pct Si) carbon steel, the density value is given as 7000 kg/m 3 at 1823 K. [24] We adopted 7000 kg/m 3 as the density value of liquid steel.
The diffusion of oxygen in steel was assessed and compiled in a handbook by institute of iron and steel Japan as D O = 2.7 9 10 À9 m 2 /s (at 1873 K). [25] The model parameters employed in the calculation are shown in Table II.

F. Model Strategy
The time step for calculation of nucleation and growth of inclusions was set to be 1 9 10 À7 seconds. The time step for solving PBE using PSG method is 1 9 10 À3 seconds. Due to the extreme small time scale for nucleation and growth calculation, the simulation of nucleation and growth for a long timescale is a very heavy job for computer. According to Zhang et al. [29] the control mechanism for evolution of inclusion size distribution should be agglomeration due to collisions after several seconds. There is a research trend to build the online model to predict the evolution of composition and size of inclusions during ladle treatment. To incorporate the present model into future complicated model, a simplified model (model S), in which nucleation and growth of inclusions is omitted after 5 seconds, was proposed to reduce the computation time. For model S, the calculation time for obtaining a size distribution at 1000 seconds after de-oxidation is around 360 seconds (6 minutes) when the program was run on a personal computer with intelÒ core i7-4770 3.4 GHz and 8G RAM. This guarantees that model S can be run in an online mode. In comparison, the calculation time for full model (F) to obtain a size distribution at 1000 seconds after de-oxidation is around 28800 seconds (480 minutes).The calculation results by model S were compared with the full model (model F) which considered all phenomena, including nucleation, growth, agglomeration and removal by floatation.

III. RESULTS AND DISCUSSION
A. Comparisons with Previous Experimental and Modelling Results

Comparison with experimental work by Nakajima et al.
The present model was applied to calculate the size distribution of alumina in Al deoxidation experiments performed by Nakajima et al. [30] and the calculation results were compared with the experimental three-dimensional size distribution of alumina inclusions after Al de-oxidation.
In their experiments, Fe-10mass pct Al deoxidizer was added into the 70 g Fe-10 mass pct Ni melts with initial oxygen content between 100 and 120 ppm, then the melt was stirred for 10 s with an alumina rod. The particle size distribution of deoxidation products was measured by scanning electron microscope (SEM) and electron probe micro-analyzer (EPMA) on filter after electrolytic extraction.
The initial oxygen content was set to be 100 ppm. The calculated equilibrium oxygen content is 2.8 ppm. Due to the lack of exact stirred intensity condition, different values (e = 0.2, 0.11 and 0.012 m 2 s À3 ) of dissipation rate of turbulence kinetic energy were employed in the present simulation. No specific information on the geometry of crucibles employed in Nakajima's experiments was found. It is assumed that the diameter of crucible is 0.03 m, and then a depth of liquid steel pool can be calculated by assuming a density value of 7000 kg/m 3 . Therefore, we assumed a depth of 0.014 m for inclusion flotation. Figure 2 showed the calculated size distributions by model S for alumina inclusions employing different dissipation rate of turbulence kinetic energy values compared with experimental data by Nakajima et al. The calculated size distribution by model F with Density of Alumina Inclusion 3823 kg/m 3 [ 26] Molar Mass of Alumina 0.10196 kg/mol [26] Dynamic Viscosity of Liquid Steel 0.0067 Pa.s [27] Diffusion Coefficient of Oxygen in Steel Melts Do at 1873 K 2.7 9 10 À9 m 2 /s [25] Interfacial À 1202000 + 386.3T J/mol [28] Temperature 1873 K -e = 0.0122 m 2 s À3 is also showed in Figure 2. As can be seen from Figure 2, the calculated size distribution for alumina inclusions using e = 0.0122 m 2 s À3 fits very well with the experimental size distribution by Nakajima et al. [30] in the part with radius larger than 0.07 lm.
There are some deviations between calculated and experimental size distribution for radius less than 0.07 lm, which should be due to the uncertainty of electrolytic extraction. As dissipation rate of turbulence kinetic energy increases, the maximum number density decreases, and the size distribution curve becomes wider. The peak-value radius (the radius corresponding to the maximum number density) only slightly varied with the dissipation rate of turbulence kinetic energy. The comparison between the calculated size distribution by model S and distribution by model F shows that the calculated results by both model are in good agreements for inclusions with radius larger than peak-value. This indicates that the agglomeration of inclusions dominates the size distribution of large inclusions.

Comparison with simulation by Zhang and Lee
Zhang and Lee [6] proposed a nucleation-growth model and simulated the evolution of size distribution of alumina inclusions. It was supposed that the initial oxygen is about 300 ppm and that the final oxygen content is roughly 3 ppm. The value of turbulent energy dissipation rate employed by them was unknown.
In the present simulation, we employed the same parameters in Table II. The initial oxygen content is assumed to be 300 ppm. The final equilibrium oxygen was calculated to be 2.77 ppm with Al content being 0.3 mass pct. The chemical composition of investigated steel is listed in Table III. The turbulent energy dissipation rate is assumed to be 0.0122 m 2 s À3 , by which the present model will give a size distribution approaching that by Zhang and Lee. The depth for inclusion flotation is set to be 1 m considering the industrial conditions in ladle.
The size distribution of alumina inclusions at 60 and 600 seconds was calculated using model S and model F and compared with the simulation results by Zhang and Lee. Figure 3 showed the calculated size distribution at 60 and 600 seconds and size distribution by Zhang and Lee. It could be seen from the figure that the calculation results by model S are in good agreements with the results by Zhang and Lee. A sudden decrease of number density could be found in inclusion size distribution at 60 seconds calculated by Zhang and Lee, which could be due to the assumption of immediate removal of inclusion larger than 50 lm. The calculation results by model F fit well with results by Zhang and Lee for inclusions with radius larger than peak-value radius. However, the peak number density by model F is less and peak-value radius is larger than results by Zhang and Lee. It should be mentioned that Zhang and Lee employ a simplified method to divide the total spectrum into discrete regime and sectional regime and the nucleation and growth of inclusions may not be taken into account in full accuracy.

Comparison with simulation by Zhang and Pluschkell
Zhang and Pluschkell [4] simulated the nucleation and growth kinetics using a numerical method. In their simulation, the total oxygen content before Al deoxidation is 300 ppm, and the final oxygen content is around 3 ppm at equilibrium. The turbulent energy dissipation rate was 0.0122 m 2 s À3 . It should be mentioned that they did not consider the Stokes collisions in their simulation. Therefore, Stokes collision is also not considered in the present calculations to maintain the consistency.
In the present calculation, the initial oxygen and Al content were set to be 300 ppm and 0.2 pct respectively. The calculated equilibrium oxygen content is around 2.8 ppm. The composition of other elements are as same as those in Table III. The turbulent energy dissipation rate was also set to be 0.0122 m 2 s À3 . The depth for inclusion flotation is also assumed to be 1 m. As for the size distribution, a comparison between the present calculation and results by Zhang and Pluschkell are shown in Figure 4. The number density values calculated by the model S and model F are generally in good agreements with that by Zhang and Pluschkell. Both model S and model F give close results on inclusions with radius larger than peak-value radius. However, the model F predict larger peak-value radius and smaller peak values of number density than model S, indicating the effect of nucleation and growth.

Effect of interfacial tension between steel and alumina
Interfacial tension between steel and alumina inclusion is a key parameter for the nucleation of inclusions. It could be seen from Eqs. [8] and [10] that interfacial tension determines the critical size for nucleation and the Gibbs energy change for nucleation, and thereby has an impact on nucleation rates. It is well known that  [29] oxygen is a surface active element in liquid iron. With the decrease of oxygen content in liquid iron, the surface tension and interfacial tension between iron and alumina would be enhanced. [30] This indicates that with the progress of de-oxidation, the interfacial tension value between alumina and steel could undergo a significant change. However, there is a lack of equations for calculating the interfacial tension between alumina and steel as a function of oxygen concentration. In simulation work by Zhang-Pluschkell [4] and Zhang-Lee, [6] the interfacial tension value between steel and alumina was set to be 0.5 N/m. Lei et al. [8] employed a quite different values (1.6 N/m) of interfacial tension. They tested different values of interfacial coefficients but failed to find a relationship between size distribution and the interfacial coefficients.
The effect of interfacial tension value between alumina and steel on the size distribution of alumina inclusion was investigated in the present work by employing four values of interfacial tension: 0.5, 1.1, 1.3 and 1.5 N/m. The initial oxygen content was set to be 300 ppm. The composition of steel is shown in Table III. The calculated size distributions of alumina inclusion at 1 seconds using different interfacial values were shown in Figure 5. As seen in the figure, only slight change of the size distribution could be found as interfacial tension increases from 0.5 to 1.1 N/m. However, further increase of interfacial tension would lead to significant change of the calculated size distribution of alumina inclusions. On the whole, as interfacial tension between steel and alumina increases, the maximum number density decreases and the peak value  [4] of radius increases. This behavior could be explained by considering both nucleation and coarsening of the inclusions. As interfacial tension between steel and alumina increases, the energy barrier for nucleation rises, as shown in Eq. [10]. The critical size for nucleation also increases with increasing interface tension, leading to the decrease of I A value in Eq. [9]. Accordingly, the nucleation rate will be suppressed by the increasing of interfacial tension. This would lead to the decrease of maximum number density of inclusion as shown in Figure 5. Lifshitz-Sloyzov-Wagner (LSW) theory [32,33] was often employed to interpret the coarsening (Ostwald ripening) of particles. The kinetics of Ostwald ripening could be controlled by the diffusion in liquid between large and small crystals, or the growth or dissolution at interface. If we assume that Ostwald ripening during de-oxidation is controlled by the diffusion in liquid, we could have where d is the mean crystal size at time t; d 0 is the initial mean crystal size; D is effective diffusion coefficient; r SL is solid-liquid interfacial tension; V S is the molar volume of crystal; c 0 is the mass concentration of mobile species in liquid equilibrated with a crystal with infinite large size. According to Eq. [31], as interfacial tension increases, the coarsening rate k would be raised. The enhanced coarsening could be the main reason for the increase of peak-value radius with increasing interfacial tension.

Effect of diffusion coefficient
Diffusion coefficient is one of the most important parameters for crystal growth, especially for the coarsening of alumina inclusions. In previous simulations, various values of diffusion coefficients were employed. Zhang and Pluschkell [4] employed a value of 3 9 10 À9 m 2 /s for oxygen diffusion coefficient. Zhang and Lee [6] used a value of 2.7 9 10 À9 m 2 /s which is close to that by Zhang and Pluschkell. However, Lei et al. employed Stokes-Einstein relation to obtain a much lower value of 9 9 10 À10 m 2 /s compared with previous simulations.
In the present work, the effect of diffusion coefficient values on the size distribution of alumina inclusions was investigated by employing two diffusion coefficient values of 9 9 10 À10 and 3 9 10 À9 m 2 /s. The composition of steel is shown in Table III. The calculated size distributions were shown in Figure 6. As seen in Figure 6, the maximum number density decreases and the peak-value radius moves to higher values, as the diffusion coefficient increases from 9 9 10 À10 to 3 9 10 À9 m 2 /s. The number density values before peak calculated with D = 9 9 10 À10 m 2 /s are one order of magnitude higher than that with D = 3 9 10 À9 m 2 /s . This indicates that the coarsening of inclusion particles was enhanced by increasing diffusion coefficient. This could be also interpreted by LSW theory. As shown in Eq. [31], the increase of diffusion coefficient would raise the coarsening rate. Accordingly, the maximum number density decreases, and the peak-value radius increases, as the interfacial tension between steel and inclusion increases. It should be mentioned that Lei et al. [8] also investigated the effect of diffusion coefficient on size distribution using their model, but they found very weak effect.

Effect of initial oxygen content
The effect of initial oxygen content on the size distribution of alumina inclusion at 60 s was investigated by employing the model S. Four levels of initial oxygen content (100, 200, 300,400 ppm) were considered. The aluminum content is kept to be constant as 0.1 pct. The compositions of other elements in steel are as same as those in Table III. The energy dissipation rate for all calculation is set to be 0.0122 m 2 s À3 and the depth for inclusion flotation is 1 m. The calculation results were shown in Figure 7. As seen in Figure 7, the number density of inclusions is decreased with increasing the initial oxygen content. However, the calculated largest inclusion increases in size at enhanced initial oxygen content. For the steel with 100 ppm oxygen, only inclusions smaller than 10 lm can be found. In contrast, inclusion larger than 100 lm can be observed in steels with initial oxygen content larger than 300 ppm. This reflects the enhanced coarsening and agglomeration at higher oxygen content. The enhanced coarsening with rising oxygen content can be well understood from Eq. [31]. As concentration of oxygen increases, the coarsening rate increases. The enhanced coarsening would lead to higher number density of large inclusions, which further promote the agglomeration of inclusions. Accordingly, larger inclusions were observed in steels with higher initial oxygen content.

Effect of the depth for inclusion flotation
The influence of the depth for inclusion flotation, H, on the size distribution of inclusions was also studied. The initial aluminum and oxygen content is assumed to be 0.1 pct and 300 ppm respectively. The compositions of other elements in steel are as same as those in Table III. The energy dissipation rate is also set to be 0.0122 m 2 s À3 . The size distributions of inclusions at 60 seconds were calculated by assuming H = 0.014 and 1 m and shown in Figure 8. As can be seen in Figure 8, the number density of the inclusion which is smaller than 40 lm is hardly affected by the depth for inclusion flotation. As H decreases, the number density of inclusions which are larger than 40 lm decreases due to the stokes flotation. The present study indicates that the Stokes flotation only has a significant effect on large inclusions.

C. The Evolution of Inclusion Size Distribution with Time
The evolution of supersaturation degree and inclusion size distribution with time after deoxidation is simulated on a steel grade. The initial oxygen content is 400 ppm. The aluminum content is 0.2 mass pct. Then the calculated equilibrium oxygen content is 3.16 ppm. The composition of other elements in steel is as same as those in Table III.
The evolution of supersaturation degree with time was calculated using the present model and shown in Figure 9. It could be seen from the figure that the supersaturation degree only drops slightly during the first 10 ls because the number of nuclei precipitated is still small. After 2 9 10 À5 seconds, the supersaturation degree drops rapidly and approaches 1 at 3 9 10 À4 seconds. After then, the supersaturation changes slightly again. The present calculation result indicates the nucleation process of alumina inclusion is very fast and can be nearly completed within several hundred microseconds. Zhang and Lee [4] also simulated the evolution of supersaturation degree with time and find that nucleation process could be accomplished within ten microseconds. It is probable that they overestimated the nucleation rate due to a simple nucleation model they employed. The size distributions of inclusions at 1, 60, 100, 600 and 1000 second after de-oxidation were calculated using model S and model F. The dissipation rate of turbulent energy of 0.0122 m 2 s À3 was employed in calculation. The calculated size distributions of inclusions using model S and model F were shown in Figure 10. The removal of inclusion due to Stokes flotation was also considered. The depth of inclusions in steel vessel was set to be 1 m.
As seen in Figure 10, the inclusion size distribution at 1 second is close to a lognormal distribution. It has been reported by many researchers that the nucleation and growth process would generate a log normal type particle size distribution. [6,12] As time approaches 60 s, the size distribution curves become more flat. The left part of the inclusion size distribution curve still obeys lognormal distribution. In comparison, the middle part of the inclusion size distribution curve shows a straight line, indicating that the part of distribution has changed from log normal to fractal. This behavior has already been observed in previous simulations and experiments and could be due to the agglomeration of inclusion. [29][30][31] The middle part of the inclusion size distribution can be further analyzed by fitting the distribution with the fractal distribution function: f r ð Þ ¼ Cr D ; where r is the radius of inclusions, C and D are the intercept and the slope for the straight line shown in logarithmic coordinates. In order to compare with previous results better, the particle density function (PDF) was employed in fitting. Zinngrebe et al. [34,35] proposed that while the size distribution results can be variable depending on bin size setting, PDF results will show a constant distribution regardless of bin size setting. PDF can be calculated from particle size distribution using the following equation: where the n(r) is the number density of inclusions for bin size of Dr.
The middle parts of calculated PDF are shown in Figure 11. All PDFs show linear functions with regard to radius of inclusions in logarithmic coordinates. Accordingly, all straight lines can be fitted using the fractal distribution function. The fitting results are shown as follows: PDF r ð Þ ¼ 0:0084r À3:24 t ¼ 60 s ½33 PDF r ð Þ ¼ 0:000504r À3:42 t ¼ 100 s ½34 PDF r ð Þ ¼ 0:000234r À3:45 t ¼ 600 s ½35 PDF r ð Þ ¼ 0:000128r À3:48 t ¼ 1000 s ½36 Zinngrebe et al. [34] characterized the inclusion populations in 120 Ti-alloyed Al-killed steel samples using automated scanning electron microscope and applied PDF in analyzing the inclusion populations of samples during secondary metallurgy. They reported an exponent value of À 3.5 for all PDF plots, which is very close to the present value of À 3.48 for t = 1000 s. Seo et al. [36] investigated the size distribution of inclusions in  Ti-alloyed Al-killed low-carbon steel samples. It was found that the exponents in fractal distribution were À 3.54 near the end of RH process at 7-13 minutes after deoxidization, À 2.76 in tundish and À 2.38 in cast slab, respectively. The exponent value in RH process reported by them is very close to the present results. They also calculated the PDF values for simulation results by Zhang and Lee [6] and found the exponent values of À 3.16(60 s), À 3.55(120 s) and À 3.60(300 s). The similar change tendency of exponent value in fractal distribution with treating time, i.e., decreasing with time, is also found in the present study. The opposite number of exponent value in fractal distribution can be understood as the fractal dimension value of the system. According to the theory of fragmentation, [37,38] fractal dimension is a measure of the fracture resistance of the material relative to the process causing fragmentation.
With the progress of agglomeration, the inclusion aggregates become more complicated and stable, and the fracture resistance of inclusion aggregates should increase, leading to higher value of fractal dimension. The final part of the curve drops quickly as radius increases, which is attributed to the removal of large inclusions. The maximum and total number density significantly decreases compared with that at 1 seconds. As time further increases, the number density decreases, which is mainly due to the agglomeration by Brownian, turbulent and Stokes collisions. The number of large inclusion continuously decreases due to the removal of inclusion by the flotation towards top slag layer.
Compared with model S, model F gives very close results on inclusions with radius higher than peak-value radius. For model F, the calculated maximum peak number density by model F is lower and the calculated peak-value radius given by model F is higher than those by model S. It is well known that larger inclusions are more harmful to the mechanic properties of steels compared with tiny inclusions. In spite of the difference between model S and model F, model S is still applicable for calculating inclusion size distribution due to its good performance in calculating the size distributions of large inclusions.
IV. SUMMARY 1. The nucleation, growth and removal of inclusions in molten steel were modeled by combining Kampmann-Wagner numerical (KWN) model with particle size grouping method. The model could simulate the time evolution of size distribution of alumina inclusions after aluminum de-oxidation. 2. The model was validated by using the experimental size distribution data of alumina inclusions by Nakajima et al. The model calculation results were also compared with previous simulation results. 3. The influences of interfacial tension between steel and inclusion, diffusion coefficient and initial oxygen content on the calculated inclusion size distribution were investigated. As diffusion coefficient increases, the maximum number density decreases and the peak-value radius increases. As interfacial tension between steel and alumina increases, the maximum number density decreases and the peak value of radius increases. The calculated largest inclusion increases in size with increasing initial oxygen content. The Stokes flotation of inclusions only has a significant effect on large inclusions. 4. The calculation results showed that the nucleation process of alumina inclusions is very fast and can be accomplished within several hundred microseconds. 5. The calculated size distribution curves showed a change from log normal to fractal, which is due to the change of dominating mechanisms for crystal growth and agglomeration. The exponent value in calculated fractal distribution decreases as treating time increases, which can be explained by the theory of fragmentation.

ACKNOWLEDGMENTS
Financial supports from the Academy of Finland for Genome of Steel Grant (No. 311934) is gratefully acknowledged. FUNDING Open access funding provided by University of Oulu including Oulu University Hospital.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.