Numerical Analysis of E ﬀ ect of Initial Bubble Size on Captured Bubble Distribution in Steel Continuous Casting Using Euler-Lagrange Approach Considering Bubble Coalescence and Breakup

: A mathematic model considering the bubble coalescence and breakup using the Euler-Lagrange approach has been developed to study the e ﬀ ect of the initial bubble size on the distribution of bubbles captured by the solidiﬁcation shell. A hard sphere model was applied for dealing with the bubble collision. Advanced bubble coalescence and breakup models suitable for the continuous casting system and an advanced bubble captured criteria have been identiﬁed established with the help of user-deﬁned functions of FLUENT. The predictions of bubble behavior and captured bubble distribution agree with the water model and plant measurements well respectively. The results show that the number of small bubbles captured by solidiﬁcation shell is much higher than that of large bubbles. What is more, the number of captured bubbles at the sidewalls decreases with the distance from the meniscus. For the case of large gas ﬂow rate (gas ﬂow fraction of 8.2%), the initial size of bubbles has little e ﬀ ect on bubble captured distribution under various casting speeds. When the gas ﬂow rate is small (gas ﬂow fraction of 4.1%), the number density of captured bubbles increases as the initial bubble size increases, and the e ﬀ ect of initial bubbles size on captured bubble number density is ampliﬁed when the casting speed decreases. The average captured bubble diameter is about 0.12–0.14 mm. Additionally, for all cases, the initial bubble size hardly a ﬀ ects the average size of captured bubbles.


Introduction
The bubbles captured by the solidifying steel shell during the continuous casting process are the main cause of defects of blisters and slivers in the slab. Argon gas is usually injected from the submerged entry nozzle (SEN) to withstand nozzle clogging [1,2]. The transport of argon bubbles has a great effect on the product quality. The bubbles escaping at the slag layer, while not casting the slag entrapment, are helpful for removing the inclusions by attaching them to the slag, which is good for the slab quality. However, the bubbles captured near the meniscus cause surface defects, and those captured deep in the mold result internal defects, which is harmful for the slab quality [3].
Argon gas is injected through the porous refractory at the upper part of the SEN wall. The initial size of argon bubble is affected by the refractory properties, gas flow rate and fluid velocity [4]. After entering the mold, both the evolution and redistribution of bubbles have great impact on the fluid flow and product quality. Due to the high temperature and invisibility, it is difficult to study the phase is considered in the model, which is neglected in the DPM (Discrete Phase Model). The collision of bubbles is modeled by means of the hard sphere model. The liquid phase is solved with the commercial package FLUENT 18.1 (ANSYS Inc., Canonsburg, PA, USA, www.ansys.com). The coupling between two phases is realized by computing the liquid phase volume of each grid and adding interfacial forces as a source term, with the help of extensive user-defined functions of FLUENT. The governing equations of two-phase are summarized in Table 1. Table 1. Governing equations for the model.

Bubble Coalescence Model
In the present work, the bubble collision was resolved by the hard sphere model. By integrating the works of Yang et al. [29,30], a bubble coalescence model considering the bubbles size and off-center degree was used in the present work to determine if the bubble collision will result in coalescence or not. With the equivalent diameter of two collision bubbles smaller than 2.3 mm, bubbles coalesce when We 1 − B 2 < 0.16 + B/24 where We is the relative Weber number; r e = 2r 1 r 2 /(r 1 + r 2 ) is the equivalent radius of two interacting bubbles and v r is the relative velocity of bubbles.
where B is the dimensionless number characterizing the off-center degree. Where r 1 and r 2 are the radii of the two bubbles and b is the distance from the center of the one bubble to the relative velocity vector originating from the center of the other bubble at contact. With the equivalent diameter of two collision bubbles larger than 2.3 mm, two collision bubbles result in coalescence if the relative approach velocity exceeds 0.11 m/s [31,32].

Bubble Breakup Model
The bubble breakup model used in the present work is the model proposed by Yang et al. [30]. In the turbulent two-phase flow, a bubble breaks up into two daughter bubbles if the uneven pressure acting on the bubble surface larger than the surface force [33]. The breakup of the bubble occurs if the bubbles exceed a maximum stable size, d max .
where We crit is the critical Weber number for bubble break-up with a value of 0.53; ε is the turbulent kinetic energy dissipation rate and σ is the gas-liquid surface tension. Besides the bubble breakup criteria, the volume fraction of the daughter bubble is the other essential parameter. According to the work of Yang et al. [30], the volume fraction of daughter bubble suit for the tanh function, following the U-shape of daughter size distribution, in the steel continuous casting system.
where f bv is the volume fraction of daughter bubble and x is a random variable from 0 to 1.

Bubble Capture Model and Predict Sample Observations
The bubble capture model as proposed by Yuan and Thomas [25,26,28] was used in the present work. This model considers the effect of PDAS (primary dendrite arm spacing), bubble size and local flow field. If the diameters of bubbles are smaller than the PDAS, the bubbles will be captured when touching the solidification shell. When the bubbles are larger than PDAS, the force balance is considered to determine if bubbles are captured by the solidification shell when the local flow velocity changes. Besides the forces introduced in the previous section, three additional forces, including the lubrication force F lub , surface tension gradient force F Grad and Van der Waals force F IV , were exerted on the captured result.
where R b is the bubble radius; u sol is the solidification front moving velocity and the details about it can be found elsewhere; h t is the thickness between the bubble and dendrite tip and r d is the radius of dendrite tip.
where α = 1 + nC 0 , β = (C * − C 0 )nr d , ξ = R b+ + r d + h o and m, n, C 0 and C * are the model parameters, which are given in the previous work.
a o is the steel atomic diameter, σ sb , σ sl and σ bl represent surface tensions for solid-bubble, solid-liquid and bubble-liquid respectively.
For comparing with plant measurement, it is essential to determine which layer of product steel the bubbles appear on. Whether the captured bubble q with radius R q at the distance D q beneath the casting product surface can be found in the examined layer of the sample w at the distance D w beneath the surface of test sample is shown in Figure 1. The solidification shell thickness is determined by t [28]. The solidification time t is calculated as t = L/v c = 40L, where L is the length beneath the meniscus and v c is the casting speed. So, D q = S| L q = 3 40L q , where L q is the distance between the entrapped bubble q and the meniscus.
solid-liquid and bubble-liquid respectively. For comparing with plant measurement, it is essential to determine which layer of product steel the bubbles appear on. Whether the captured bubble q with radius Rq at the distance Dq beneath the casting product surface can be found in the examined layer of the sample w at the distance Dw beneath the surface of test sample is shown in Figure 1. The solidification shell thickness is determined by 3 s t = [28]. The solidification time t is calculated as

Numerical Implementation
To quantitatively validate the mathematic model, the benchmark case in the paper of Jin et al. [28] was built. The computational domain and boundary conditions are shown in Figure 2. Table 2 presents the geometry parameters and operation conditions. The liquid phase flow field was solved with the commercial package FLUENT. The model of bubble was written by C. The couple between two phases was realized by computing the liquid phase volume of each grid and adding interfacial forces as a source term, with the help of extensive user-defined functions of FLUENT. The time step was set as 0.001 s, and the total simulation time was 60 s. The bubbles were assumed as spheres with the initial diameter ranging from 1 to 2 mm at this operation condition, based on Lee et al. [4].

Numerical Implementation
To quantitatively validate the mathematic model, the benchmark case in the paper of Jin et al. [28] was built. The computational domain and boundary conditions are shown in Figure 2. Table 2 presents the geometry parameters and operation conditions. The liquid phase flow field was solved with the commercial package FLUENT. The model of bubble was written by C. The couple between two phases was realized by computing the liquid phase volume of each grid and adding interfacial forces as a source term, with the help of extensive user-defined functions of FLUENT. The time step was set as 0.001 s, and the total simulation time was 60 s. The bubbles were assumed as spheres with the initial diameter ranging from 1 to 2 mm at this operation condition, based on Lee et al. [4].

Bubble Distribution
The bubble distribution in the mold from the computation model was compared with a simple water model results in Figure 3. From Figure 3, it can be noticed that the entire bubble distribution from the computation model agreed qualitatively well with the water model. Most of the bubbles floated up near the SEN due to the buoyancy, while smaller bubbles got to the further location following with the liquid flow. Due to the introduction of bubble interaction, the bubble size distribution was controlled by the bubble coalescence and breakup. In the submerged entry nozzle, the bubble was small at the initial size and the number of bubble was large, so a lot of bubbles coalesced to larger ones, as shown in Figure 3. At the location near the SEN ports, due to the high turbulent stress, and the large bubble size causing by amounts of coalescence of bubbles in the SEN, a lot of bubbles broke up, as shown in Figure 3. It can be noticed that the bubble distribution and bubble interaction agreed well in the mold with the experiments results. However, there were more bubbles at the top surface in the water model than the simulation results. The bubbles that got to the gas-liquid interface would stay for a while in the water model resulting in the aggregation of bubbles, but the bubbles that get to the interface will be removed immediately by the program setting in the mathematic simulation. Metals 2020, 10, x FOR PEER REVIEW 7 of 15

Comparison of Bubbles Captured in Samples
For validating the computation result, the measurement of samples with six layers examined on each by Jin et al. [4], shown in Figure 4a, were compared with the computation results. For the plant measurements, the part above the examining layer was milled away, and then an optical microscope was used to record the bubble diameter and number of the examining layer [4]. Figure 4b,c shows comparison of bubble number density and mean diameter between the simulation and measurement at the wide face. It can be seen that the computation simulation predicted 0.1-0.2 bubble per cm 2 on the first four layers at the wide face, which matched the measurements of the sample well. On the fifth and sixth layers, there were less than 0.05 bubble per cm 2 . The bubbles in these two layers were entrapped at the downward recirculation region, where only a few small bubbles could reach. Figure 4c shows the average bubble diameter at the wide face. The simulation result shows that the average bubble diameter slightly decreased and the distance beneath the surface was 0.1-0.15 mm, which also matched the measurement. Figure 4d,e shows the bubble number density and the bubble diameter of simulation and measurement at the narrow face. The size and number of bubbles captured at the narrow was similar with the wide face, and the simulation results matched the measurements.

Comparison of Bubbles Captured in Samples
For validating the computation result, the measurement of samples with six layers examined on each by Jin et al. [4], shown in Figure 4a, were compared with the computation results. For the plant measurements, the part above the examining layer was milled away, and then an optical microscope was used to record the bubble diameter and number of the examining layer [4]. Figure 4b,c shows comparison of bubble number density and mean diameter between the simulation and measurement at the wide face. It can be seen that the computation simulation predicted 0.1-0.2 bubble per cm 2 on the first four layers at the wide face, which matched the measurements of the sample well. On the fifth and sixth layers, there were less than 0.05 bubble per cm 2 . The bubbles in these two layers were entrapped at the downward recirculation region, where only a few small bubbles could reach. Figure 4c shows the average bubble diameter at the wide face. The simulation result shows that the average bubble diameter slightly decreased and the distance beneath the surface was 0.1-0.15 mm, which also matched the measurement.    Figure 5 shows the fluid velocity counter and vector in the mold at the center cross-sectional plane. It can be seen that two jet flows were in the upper recirculation region. One jet flow impacted into the narrow face, and then it branched into an upper and a down flows. The other jet flow got to the top surface after out of the port, which was driven by the bubble buoyancy.

Two-Phase Flow in the Mold
Metals 2020, 10, x FOR PEER REVIEW 9 of 15 Figure 5 shows the fluid velocity counter and vector in the mold at the center cross-sectional plane. It can be seen that two jet flows were in the upper recirculation region. One jet flow impacted into the narrow face, and then it branched into an upper and a down flows. The other jet flow got to the top surface after out of the port, which was driven by the bubble buoyancy. Figure 6 shows the whole instantaneous bubble distribution in the continuous casting mold after 25 s. In Figure 6, the bubble size is presented with color for a better visualization of the small bubbles. It can be noticed that the large bubbles (green color) floated up to the top surface close to the SEN. Additionally, only the bubbles (blue color) smaller than 1 mm could get to the narrow face and the deep of the mold, with the reason that the effect of drag force on the small bubbles was greater than that of buoyancy. Most of the bubbles captured by the solidification shell were these small bubbles.    Figure 6 shows the whole instantaneous bubble distribution in the continuous casting mold after 25 s. In Figure 6, the bubble size is presented with color for a better visualization of the small bubbles. It can be noticed that the large bubbles (green color) floated up to the top surface close to the SEN. Additionally, only the bubbles (blue color) smaller than 1 mm could get to the narrow face and the deep of the mold, with the reason that the effect of drag force on the small bubbles was greater than that of buoyancy. Most of the bubbles captured by the solidification shell were these small bubbles.

Distribution of Bubble Captured Location
Metals 2020, 10, x FOR PEER REVIEW 9 of 15 Figure 5 shows the fluid velocity counter and vector in the mold at the center cross-sectional plane. It can be seen that two jet flows were in the upper recirculation region. One jet flow impacted into the narrow face, and then it branched into an upper and a down flows. The other jet flow got to the top surface after out of the port, which was driven by the bubble buoyancy. Figure 6 shows the whole instantaneous bubble distribution in the continuous casting mold after 25 s. In Figure 6, the bubble size is presented with color for a better visualization of the small bubbles. It can be noticed that the large bubbles (green color) floated up to the top surface close to the SEN. Additionally, only the bubbles (blue color) smaller than 1 mm could get to the narrow face and the deep of the mold, with the reason that the effect of drag force on the small bubbles was greater than that of buoyancy. Most of the bubbles captured by the solidification shell were these small bubbles.

Distribution of Bubble Captured Location
The bubble captured locations at the narrow and wide faces are discussed in this section. Different from the previous papers with a constant bubble size, the number of small bubbles captured by the solidification shell was determined by the actual flow field and the bubble interaction. The distribution of the bubble entrapped location under 1.5 m/min casting speed and a bubble fraction of 8.2% is shown in Figure 7. For a better analysis of the bubble entrapped location, the bubble was divided into three levels: small bubble (diameter smaller than 0.2 mm), medium-size bubble (diameter ranges from 0.2 to 1 mm) and large bubble (diameter larger than 1 mm). From Figure 7, it can be seen that many small bubbles were captured by the solidifying shell, and the number of bubbles decreased with the distance from the top surface. In the narrow face, a number peak of captured bubbles appeared at the region about 0.5 mm below the top surface, where the mainstream impacts. In the deep part of the wide face, more bubbles were captured near the narrow face. From Figure 7, it can be noticed that the number of entrapment bubbles with medium-size bubbles and large size was far less than the small bubbles. There was a small amount of medium-size bubbles captured in the deep, but almost no large bubbles were captured in the deep of the mold. Additionally, the most of the large bubbles were entrapped close to the top surface. The bubble captured locations at the narrow and wide faces are discussed in this section. Different from the previous papers with a constant bubble size, the number of small bubbles captured by the solidification shell was determined by the actual flow field and the bubble interaction. The distribution of the bubble entrapped location under 1.5 m/min casting speed and a bubble fraction of 8.2% is shown in Figure 7. For a better analysis of the bubble entrapped location, the bubble was divided into three levels: small bubble (diameter smaller than 0.2 mm), medium-size bubble (diameter ranges from 0.2 to 1 mm) and large bubble (diameter larger than 1 mm). From Figure 7, it can be seen that many small bubbles were captured by the solidifying shell, and the number of bubbles decreased with the distance from the top surface. In the narrow face, a number peak of captured bubbles appeared at the region about 0.5 mm below the top surface, where the mainstream impacts. In the deep part of the wide face, more bubbles were captured near the narrow face. From Figure 7, it can be noticed that the number of entrapment bubbles with medium-size bubbles and large size was far less than the small bubbles. There was a small amount of medium-size bubbles captured in the deep, but almost no large bubbles were captured in the deep of the mold. Additionally, the most of the large bubbles were entrapped close to the top surface.

Effect of the Initial Bubble Size on the Distribution of Captured Bubbles
The initial size of bubbles may affect the bubble secondary size distribution after the interaction of bubbles and the complex effect by the flow field, which in turn affects the capture of bubbles by the solidification shell. The effect of the initial size of bubbles on the distribution of captured bubbles may vary with different casting speeds and gas flow rates. For analyzing the distribution of captured bubbles, number density and average size of captured bubbles are discussed in this section. The number density of captured bubbles was calculated as the ratio of the number of captured bubbles to the area during the casting time. The area during the casting time was calculated as: S k = v c × t c × L k , where S k is the area during the casting time, v c is the casting speed, k is the wide or narrow face, t c is the casting time and L k is the length of wide face or narrow face. From the work of Lee et al. [4], few initial bubbles generated from the porous refractory are smaller than 1 mm. So, to investigate the effect of initial bubble size on the captured bubble distribution, three different initial bubble sizes of 1 mm, 1.5 mm and 2 mm were compared in this section. Table 3 shows the distribution of captured bubbles at the wide and narrow faces with three different initial bubble sizes under the operation of the casting speed of 1.5 m/min and gas flow rate of 8.2% gas volume fraction. The number density of captured bubbles on the wide face and narrow face was about 3.6 and 5.2 per cm 2 , which clearly indicates that there were more bubbles per unit area on the narrow face. Since the main fluid stream rushed to the narrow face, more bubbles were carried to the narrow face. The average size of captured bubbles was about 0.130 and 0.126 mm respectively. Due to the greater fluid velocity towards the narrow face, a larger drag force acted on the bubbles facing the direction to the narrow face, resulting in larger bubbles brought to the narrow face. With the different initial bubble sizes, there were no obvious differences in the number density and size of captured bubbles. In the operation of casting speed of 1.5 m/min and gas flow rate of 8.2% gas volume fraction, the gas volume was large enough and the flow field was sufficiently turbulent, so the bubbles in the SEN with different initial bubble sizes could fully collide. The number and size of large bubbles, which could break up to the small bubbles that could be captured by the solidified shell, had little difference under different initial bubble sizes. So, there were no obvious differences at the number density and size of captured bubbles with different initial bubble sizes. The number density and size of captured bubbles with the casting speed of 1.7 m/min and 1.3 m/min and gas flow rate of 8.2% gas volume fraction are present in Tables 4 and 5. It can be noticed that more bubbles were captured and the average bubble size decreased with increasing casting speed. It also can be found that when the gas flow rate was 8.2%, the initial bubble size had no obvious effect on the captured bubble distribution at the operation of the casting speed of 1.7 and 1.3 m/min. The reason was same as the casting speed of 1.5 m/min.  Table 6 shows the number density and average size of captured bubbles with different initial bubble sizes. For all cases, the gas flow rate was 4.1% argon volume fraction and the casting speed was 1.5 m/min. The results indicate that the bubble number density decreased comparing with the gas flow rate of 8.2% argon volume fraction. For different initial bubble sizes, it can be found that the captured bubble number density decreased with the decrease of initial bubble size. The captured bubble number density with the initial bubble size of 2 mm was about 1.21 times of that of 1 mm. For the small gas flow rate (argon volume fraction of 4.1%), the probability of collision between bubbles was lower than that of a large gas volume fraction. Therefore, in the SEN, the bubbles with a large initial size could coalesce into more large bubbles, which could break up into the small bubbles that were captured by the solidified shell, than small initial size bubbles. So, the captured bubble number density with the large initial bubble size was larger than that with a small initial bubble size. The effect of the initial bubbles size on captured bubble distribution under the condition of casting speed of 1.7 m/min and gas flow rate of 4.1% gas volume fraction is presented in Table 7. It can be found that the captured bubble number density with the initial bubble size of 2 mm was about 1.05 times of that of 1 mm. For the case of a gas flow rate of 4.1% argon volume fraction, the effect of the initial bubble size on the captured bubble number density decreased when the casting speed increased. By increasing the casting speed, the maximum velocity that bubbles can get will increase, and the bubble motion will be more chaotic. Therefore, the probability of collision between bubbles increases with increasing casting speed, and the difference in the number of large bubbles under different initial bubble sizes will be reduced. This is why the effect of the initial bubbles size on captured bubble distribution decreased with the increase of casting speed. From Table 8, for the case of casting speed of 1.3 m/min, it can be found that the captured bubble number density with the initial bubble size of 2 mm was about 1.55 times of that with an initial bubble size of 1 mm. Similar to the previous explanation, for the case of gas flow rate of 4.1% argon volume fraction and casting speed of 1.3 m/min, the gas volume fraction was small, and the bubbles motion was more stable with a lower casting speed. Therefore, the probability of collision between bubbles was lower with a lower casting speed. So, the effect of initial bubble size on the captured bubble number density was higher. From the simulation results, it can be noticed that the initial bubble size had no significant effect on the captured bubble size. The reason is that the average size of bubbles that can reach the advancing solidified shell was mainly determined by the casting speed.

Conclusions
The Euler-Lagrange model was used to simulate the bubble transport in the continuous casting mold and bubble captured by the solidification shell. The bubble coalescence and breakup were considered in the model, and an advanced bubble captured model was established. The location, number density and size of captured bubbles were analyzed in the present work. Besides, the effect of the initial bubble size on the captured bubble distribution under different operation conditions was discussed. The conclusions are the following: 1. Most of large bubbles floated up to the top surface close to the SEN, and only some small bubbles could get to the narrow face and the deep of the mold.
2. The number of small bubbles (diameter smaller than 0.2 mm) captured by the solidification shell was far more than medium-size (diameter ranges from 0.2 to 1 mm) and large bubble (diameter larger than 1 mm). The number of small and medium-size bubbles captured by the solidification shell decreased with the distance from the top surface. The large bubbles were captured only near the top surface.
3. For the case of a large gas flow rate, the initial bubble size had no obvious effect on the number density and average size of captured bubbles under different casting speeds. When the gas flow rate was small, the number density of captured bubbles by the solidification shell decreased with the decrease of the initial bubble size, and the effect of initial bubble size on captured bubble number density increased when the casting speed decreased. For all cases, the initial bubble size hardly affected the average size of captured bubbles.
Author Contributions: Z.Z., Z.L. and W.Y. conceived and designed the study; W.Y. and N.Z. did the numerical simulation; all the authors discussed the results; W.Y. wrote the paper; Z.Z. and Z.L. edited the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding:
No funding has supported this work.