Numerical Simulation of Bubble Behavior before Inclined Solidified Front

Abstract

In order to investigate bubble behavior before inclined solidified front, a numerical simulation model is developed, in which level set method and modified heat transfer equation is applied to simulate the dynamic evolution of gas-liquid interface and solidification process under a fixed grid frame. Meanwhile, an in-situ experiment of bubble behavior before inclined ice solidified front was exacted to validate the numerical model. The effects of bubble diameter, inclined angle of solidified front and cooling rate on its behavior are investigated both numerically and experimentally. The results show that as bubble diameter increasing, inclined angle increasing and cooling rate decreasing, the bubble entrapment possibility is decreasing. Each critical point for entrapment is also studied numerically.

1. Introduction

In the process of continuous casting, argon gas injection via submerged entry nozzle is used to effectively prevent the nozzle from clogging and suppress the secondary oxidation of molten steel.¹⁾ Furthermore, the buoyant rise of argon bubble can stir the molten steel inside the mold and results in the uniformity of thermal field and flow field. Whereas, the entrapment of gas bubbles near the advancing solidified shell causes lots of internal defects to the slab.²⁾ In continuous casting process, these defects due to bubble and inclusion are very harmful to the quality of steel sheet, especially in ultra-low carbon steel CC process.³⁾

Therefore, the research on flow field in the mold and the interaction between gas bubble and solidified front has become a hot topic,^4,5,6) which is regarded as a key factor that determines gas bubble’s entrapment or repulsion. Research is mainly focused on the entrapment mechanism by experimental methods and theoretical analysis. Esaka²⁾ used a transparent substance called Succinonitrile (SCN) to conduct an in-situ observation of bubble behavior in front of inclined solidified front. Through a series of variable controlled experiments and a physical model taking the force balance on a bubble into account, they claimed that with a bigger bubble diameter, higher temperature gradient at the solidified front and a bigger inclined angle of solidified front, the bubble would tend to be entrapped by solidified front rather than repulsed away. Meanwhile, H. Li⁷⁾ built a numerical model to study bubble behavior under different velocity distribution near the solidified front in continuous casting machine. Based on the flow field obtained by solving Navier-Stokes equations, the motion of bubble was studied by using a particle trajectory model accounting a force balance affecting on bubble. But the shortcoming is that the bubble should be regarded as an object with deformable surface rather than a rigid particle in this particle trajectory model.

Moreover, many other studies reveal that the temperature gradient as well as the thermocapillary effect caused by temperature gradient at the solidified front plays a significant role on the bubble behavior. Hadji⁸⁾ claimed that the bubble behavior is mainly determined by thermocapillary effect which is surface tension gradient caused by temperature gradient. And a corresponding physical model shows that the distant between bubble and solidified front is a function of temperature and other parameters. Meanwhile, Min.S.Park⁴⁾ considered the thermocapillary convection caused by temperature gradient and density difference as main influence factors of the bubble behavior before solidified front.

In order to numerically study the interaction between the gas bubble and solidified front, the coupling simulation of free surface evolution of gas bubble and solidification process is required, which is very difficult because the solidified front and gas bubble shape are dynamically varying simultaneously. This is a highly non-linear problem. Previously, Pasandideh-Fard^9,10,11) made numerical simulation of the process of solidification and impinging of tin and water droplet on a hot surface with Volume of Fluid method (VOF), and validated the numerical model with experimental results. The solidification process and thermal distribution are added into Navier-Stokes equation as a body force term. Both the solidified front propulsion and free surface evolution are well tracked, which gives a good example for us.

In this article, a numerical method aiming at simulating the bubble behavior in front of an inclined solidified front is developed, in which level set method^12,13,14) and modified heat transfer equation are applied to track the dynamic evolution of gas-liquid interface and the solidification process under a fixed grid frame. Based on Voller and Prakash’s theory,¹⁵⁾ a body force correction term is added to Navier-Stokes equation to describe the solidification process, in which the body force is a function of temperature-dependent porosity function. Meanwhile the latent heat content is calculated via a modified heat capacity equation which included a smooth Dirac delta function to optimize numerical calculation of solidifying interface. Besides, level set method^12,13,14) is used to describe the shape and behavior evolution of a gas bubble in the liquid phase, and the level set function is also used to distinguish the parameters of the liquid phase from gas phase. As a result, this is a case of three-phase interaction. In order to validate this mathematic model, an in-situ observation experiment of bubble behavior before advancing solidified front is conducted. Analysis and conclusion about experimental and numerical results is given through the comparison.

2. Numerical Model Description

2.1. Model Simplification and Assumptions

Figure 1 shows the process of bubble entrapment by hook structures formed at meniscus in continuous casting mold. During this process, bubble is brought to the hook structure zone at meniscus along with the steel flow field, and its entrapment or repulsion is determined by its physical properties, flow field features, and solidified front characteristics. In order to simplify the practical model, a corresponding geometrical model for numerical simulation is built, as shown in Fig. 2. In this simplified model, the liquid phase is solidifying from an inclined interface where a certain outwards heat flux exists to simulate the water cooling process outside the mold. Meanwhile a gas bubble with a certain diameter in the liquid is floating up under the effects of buoyancy and arriving at the solidified front. In this process, the growth rate of solidified shell is controlled by the heat transfer coefficient at the solidification boundary.

Fig. 1.

Sketch map of bubble behavior with meniscus in continuous casting process.

Fig. 2.

Model schematic diagram.

The assumption of this simplified numerical model is made as below:

1) The influence of steel flow from inject nozzle on bubble behavior, as well as the interaction between bubble and inclusions, is ignored.

2) The flow field near the meniscus is regarded as incompressible viscous laminar flow. The physical properties of both liquid phase and gas phase are considered as independence of temperature and constant.

3) The influences of solidified front structure, hook structure, and oscillation marks (OSM) on bubble behavior are ignored.

4) The effects of thermocapillary force dependent on thermal gradient, added mass force and Basset force which are unsteady forces due to acceleration of the bubble with respect to the fluid, Magnus force which is due to the spinning motion, are ignored.

5) In this model, the meniscus profile is simplified into a straight oblique. According to the Bikerman equation¹⁶⁾ and Takeuchi’s research,¹⁷⁾ length of the slope is chosen as the value of 10 mm which is at the same order of magnitudes with the one in practical examples.

6) The interaction between dendrite growth and bubble behavior is ignored in this model, which involves a meso-scale modeling for dendrite growth and macro-scale modeling for solidification and two-phase flow process. In this article, only the influence of temperature on growth rate of solidified shell is taken into account.

2.2. Incompressible Two-phase Flow: Level Set Method

  

$$\rho \frac{\partial \mathit{u}}{\partial t}+\rho (\mathit{u}\cdot \nabla )\mathit{u}=\nabla \cdot [-p\mathit{I}+\eta (\nabla \mathit{u}+\nabla {\mathit{u}}^{\mathrm{T}})]+\mathit{F}$$(1)

  

$$\nabla \cdot \mathit{u}=0$$(2)

  

$$\frac{\partial \varphi }{\partial t}+\mathit{u}\cdot \nabla \varphi =\gamma \nabla \cdot \left(\mathit{\epsilon }\nabla \varphi -\varphi (1-\varphi )\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right)$$(3)

In this model, the incompressible laminar fluid flow field is described by Navier-Stokes equations for flow as shown in Eqs. (1) and (2), and the gas-liquid interface is tracked by solving level set function^12,13,14) as shown in Eq. (3).

Here the level set function ϕ is used to distinguish the gas and liquid phase. In the Eq. (3), the value of ϕ takes 0 for liquid phase and 1 for gas phase, the remaining 0 to 1 for the transition region from liquid to gas. Smooth Heaviside function is used to specify the transition region to ensure that there’s no singular stiffness matrix existing during numerical simulation. The transition of properties for different phase can be obtained by solving the level set function, leading to a dynamic changing of material properties along with the phase transition. The properties transition for density, viscosity and heat conductivity can be expressed as follows:   

$$\rho ={\rho }_1+({\rho }_2-{\rho }_1)\varphi$$(4)

  

$$\eta ={\eta }_1+({\eta }_2-{\eta }_1)\varphi$$(5)

  

$$\kappa ={\kappa }_1+({\kappa }_2-{\kappa }_1)\varphi$$(6)

The body force term F in Navier-Stokes Eq. (1) includes the surface tension force, gravity force and the correction term for solidification process dependent on temperature (will be described later), in which the surface tension force acting on the bubble surface can be expressed as below:   

$${\mathit{F}}_{st}=\sigma \kappa \delta \mathit{n}$$(7)

Where F_st is a vector field of surface tension force, σ is surface tension coefficient of the liquid-gas interface, n is the unit normal vector to the interface, as defined in Eq. (8), κ is curvature of the interface as defined in Eq. (9), which depends on second derivatives of the level set function ϕ. δ is a Dirac delta function, which is zero everywhere on the real number line except zero point, resulting a concentration on the interface, and the δ-function is approximated by a smooth function according to Eq. (10).   

$$\mathit{n}={\frac{\nabla \varphi }{\left|\nabla \varphi \right|}|}_{\varphi =0.5}$$(8)

  

$$\kappa ={-\nabla \cdot \mathit{n}|}_{\varphi =0.5}$$(9)

  

$$\delta =6\left|\nabla \varphi \right|\left|\varphi (1-\varphi )\right|$$(10)

In addition, the gravity force term is given as below:   

$${\mathit{F}}_g=\rho \mathit{g}$$(11)

As labeled in Fig. 2, the boundary conditions for incompressible two-phase flow level set method are showed in the Table 1 below:

Table 1.Boundary condition settings for two-phase flow level set method.

Boundary number	Boundary type	Equation
①,②,③,④	No slip wall boundary	u = 0
⑤	Initial fluid interface boundary	ϕ = 0.5

2.3. Heat Transfer for Solidification

Similar to the level set method, the temperature distribution and solidus are specified by solving heat transfer equations rather than tracking the solidified front directly. The heat transfer model for solidification is based on Voller and Prakash’s work,¹⁵⁾ where the latent heat fusion is added to the heat capacity term. The general heat transfer equation is shown in Eq. (12).   

$$\rho \cdot c_{pT}\frac{\partial T}{\partial t}+\nabla \cdot (-\kappa \nabla T)=Q$$(12)

Meanwhile, a solid fraction function F_s(T) which depends on temperature distribution is defined as 1 for solid phase, 0 for liquid phase, and 0~1 for transition region which is also known as mushy zone. The specification is as below:   

$$F_s(T)\{\begin{array}{c}0,\\ \frac{T_m+\mathit{\epsilon }-T}{2\mathit{\epsilon }},\\ 1,\end{array}\hspace{1em}\begin{array}{c}T\ge T_m+\mathit{\epsilon }\\ T_m+\mathit{\epsilon }>T\ge T_m-\mathit{\epsilon }\\ T<T_m-\mathit{\epsilon }\end{array}$$(13)

T_m is the solidification temperature. 2ε is the temperature range of transition region from liquid to solid, which stands for the width of mushy zone for alloy solidification process. Here for pure substance, a transition region with a temperature width of 1 K is added as a smooth transition zone for numerical calculation, avoiding the sudden change of properties of materials from solid to liquid. The width mushy zone is chosen as a fairly small value to avoid its effects on bubble behavior.

An assumption is made that the process of solidification is similar to porous media with its porosity varying from1 to 0, where 1 is the liquid phase, 0 is the solid phase and the remaining 1 to 0 is the transition zone. Thus, a porosity function corresponding with solid fraction function is defined as below in Eq. (14).   

$$\lambda =1-F_s(T)$$(14)

  

$$F_s=-\mathrm{A}\mathit{u}$$(15)

  

$$\mathrm{A}=\frac{c{(1-\lambda )}^2}{({\lambda }^3+q)}$$(16)

Afterwards, a body force correction term dependent on the porosity function is added to Navier-Stokes Eq. (1) to dominate the diffusive and convective terms to force the velocity field to be zero as solidifying, which is defined in Eqs. (15) and (16). The term A is corresponding with constants c and q. In order to obtain a desired effect for body force correction term, constant c is arbitrarily chosen as a big enough value, and q is chosen as a small enough value. For the body force correction term, when temperature is higher than the liquidus, the porosity function λ is 1 and solidification body force F_s is zero so that there’s no effect on the Navier-Stokes equation; on the contrary, when the temperature is lower than solidus, the porosity function λ is 0 and solidification body force F_s is a fairly big value so that the velocity field is forced to zero, which is the expected situation for solid phase. For temperature ranging from solidus to liquidus, the transition region is considered as porosity media with a transitional body force correction term. As a result, the full body force form in Navier-Stokes Eq. (1) is as below:   

$${\mathrm{F}}_{\mathrm{x}}=\gamma \kappa \frac{\partial \varphi }{\partial x}\delta (\varphi )+{\mathrm{F}}_{s,x}\varphi$$(17)

  

$${\mathrm{F}}_y=\gamma \kappa \frac{\partial \varphi }{\partial y}\delta (\varphi )+{\mathrm{F}}_{s,y}\varphi +\rho \mathit{g}$$(18)

The influence brought by latent heat release can be included in the total heat capacity term in heat transfer Eq. (12) by redefining the effective heat capacity Eq. (19). Total heat capacity of the system consists of sensible heat h and latent heat ΔH. Sensible heat h can be calculated with Eq. (21). The treatment to latent heat can be expressed as Eq. (23), where δ^*(T) is a smooth Dirac delta function which is treated as a Heaviside smooth function flc2hs with a secondary continuous derivative, which is defined as Eq. (22), ensuring there’s a smooth transition region for latent heat release.   

$$c_{pT}=c_p+c_{pH}$$(19)

  

$$h=c_p\cdot T$$(20)

  

$$c_{pH}=L\cdot \delta (T)=L\cdot flc2hs(T-T_m,\mathrm{\hspace{0.17em} }\mathit{\epsilon })$$(21)

  

$$flc2hs(T-T_m,\mathrm{\hspace{0.17em} }\mathit{\epsilon })=\frac{1}{1+e^{-2\mathit{\epsilon }(T-T_m)}}$$(22)

As a result, the total heat content released between temperature T₁ and T₂ is calculated by integrating the total heat capacity c_pT as shown below:   

$$H={\int }_{T_1}^{T_2}{c}_{pT}\cdot dT={\int }_{T_1}^{T_2}{c}_p\cdot dT+L\cdot {\int }_{T_1}^{T_2}{\delta }_{2\mathit{\epsilon }}^{*}(T)\cdot dT$$(23)

As labeled in Fig. 2, the boundary conditions for heat transfer are showed in the Table 2 below:

Table 2.Boundary condition settings for heat transfer.

Boundary number	Boundary type	Equation
①	Heat flux boundary	$-\mathit{n}\cdot (-\kappa \nabla T)=q_0+h(T_{inf}-T)$
②,③,④	Thermal insulation boundary	$-\mathit{n}\cdot (-\kappa \nabla T)=0$
⑤	Thermal conductivity boundary	$-{\mathit{n}}_u\cdot (-{\kappa }_u\nabla T_u)-{\mathit{n}}_d\cdot (-{\kappa }_d\nabla T_d)=0$

2.4. Parameters and Numerical Solution

The basic parameters of materials and constant terms are shown in Table 3 and Table 4 respectively.

Table 3.Parameters of materials.

Property	Gas	Water	Ice	Description
ρ, kg/m3	1.3	1×103		Density
μ, Pa·s	4.79×10–5	1.79×10–3	—	Viscosity
k, W/(m·K)	0.024	0.6		Thermal conductivity
σ, N/m	0.072		—	Surface tension coefficient
Tm, K	—	273.15		Melting point
2ε, K	—	1		Range of transient zone
L, J/kg	—	3.33×105		Latent heat of melting
cp, J/(kg·K)	1005	4200		Heat capacity

Table 4.Constant terms.

Name	Value	Description
g, m/s2	9.81	Gravity acceleration
C	16000	Parameter for function A
q	1×10–4	Parameter for function A
Tw, K	283.15	Initial temperature of water
h, W/(m2·K))	3000	Heat transfer coefficient
Tinf, K	300	External temperature

This numerical solution is a fully coupled problem including two-phase laminar flow level set method and general heat transfer. In the domain, a mesh grid with 187339 triangle elements and 1315273 degrees of freedom is created from the grids with finer mesh elements on the boundary, where the maximum element size is 0.05 mm growing to coarse one in the center. Meanwhile, several mesh parameters are chosen to validate mesh dependence problem. A series of parametric sweep of mesh size shows that the solution is independent of mesh construction. The time consumption is the only relative factor. The numerical models are built and solved in COMSOL Multiphysics 3.5a. Almost 40 hours are needed to solve this problem with a time-step of 2.0 s with 187339 triangle elements and 1315273 degrees of freedom.

3. Model Validation and Simulation Results

3.1. In-Situ Observation Solidification Experiment

3.1.1. Experiment Setup

In order to validate the numerical model of bubble behavior before inclined solidified front, an in-situ experiment of the interaction between gas bubble and inclined ice solidified front is conducted.

Figure 3 is a sketch map of the in-situ observation experiment system, which consists of low-temperature liquid nitrogen system, solidification system, gas bubble generator and container. The solidification system includes a copper probe with a cavity whose angle with the horizontal line is θ. Liquid nitrogen is circulating with a certain flux in the cavity to constantly chill down the copper probe. As a result, ice is supposed to solidify from the surface of cooper probe. And the solidification rate is controlled by adjusting the liquid nitrogen flux. Meanwhile gas bubble is generated by electrolyzing sodium chloride aqueous solution whose concentration is 10 mol/m³. A pair of stainless wires with diameter of 0.5 mm placed in bottom of the aqueous solution is treated as electrodes. By adjusting the electrolysis intensity, bubbles with diameters ranging from 0.1 mm to 1.0 mm will be generated at the cathode. The diameter and behavior of individual bubble is measured and recorded from the image and video clip recorded by a CCD camera (30 frames per second), as shown in Fig. 3.

Fig. 3.

Diagram for visual solidification experiment system.

3.1.2. Experimental Results

In this experiment, under a certain inclined copper probe angle and liquid nitrogen flux, the influence of the diameter of bubble on its behavior before solidified front is investigated. The gas bubble diameter and its behavior can be observed and recorded from videos and images obtained from the CCD camera.

Figure 4 shows the picture of bubble behavior before inclined solidified front from the sideview recorded by a CCD camera. White arrows are added to point out the bubble in flow field in Fig. 4, where the diameter of bubble emerged from electrode is found to be smaller than 1 mm from the image. It is observed that most of the bubbles are floating up along the solidified shell, whereas a few bubbles are trapped by the advancing inclined solidified shell.

Fig. 4.

Bubble motion before solidified front.

Subsequently, a detailed analysis about the influence of bubble size on its behavior before advancing solidified front is made. Various bubbles with diameter ranging from 0.1 mm to 1.0 mm are observed in this experiment. The experiment ends when the 50% of solidification surface is covered by entrapped bubble, because the entrapped bubbles will affect the behavior of up-coming bubbles. Within the experiment duration, the bubble entrapment or repulsion behaviors for different diameters are observed from recorded video clips and images and their entrapment numbers are also counted. As shown in Fig. 5, it is found that bubble entrapment increases as bubble diameter decreasing. For bubbles with diameter smaller than 0.5 mm, its entrapment number is relatively bigger. On the contrary, bubbles with diameter bigger than 0.8 mm lead to a relatively low entrapment number. Due to its small buoyancy force and relatively high resistance force from the advancing solidified front and thermocapillary force, those bubbles with small diameter are easier to be entrapped by solidified front.

Fig. 5.

Bubble diameter’s effect on its entrapment bubble number.

A series of experiments with increasing liquid nitrogen flux are made to investigate the relationship between bubble entrapment, cooling intensity and solid/liquid interface morphology. The solid/liquid interface morphology is observed to be much coarser with the increasing liquid nitrogen flux and transition from planar to dendrite crystal. More entrapped bubbles are also observed at the dendrite crystal interface than the ones at planar crystal interface. This phenomenon can be explained from the view of solidification theory, which is that the increasing liquid nitrogen flux can result in a larger constitutional supercooling area, and transform the morphology at the solid/liquid interface from planar to dendrite crystal.

3.2. Numerical Simulation Results

3.2.1. Numerical Validation

As a comparison with experiments and a validation to the numerical simulation, the numerical simulation uses the same parameters and properties as experiment materials and circumstances, which are shown in Table 3 and Table 4. Here in the numerical simulation, the inclined solidified front angle is 50° and the initial temperature of water is 10°C, meanwhile the initial solidified boundary has outwards heat transfer coefficient of 3000 W/(m²·K). The effect of bubble diameter on its behavior is studied in different models. The diameter is taking values of 2.0 mm, 1.8 mm, 1.6 mm, 1.4 mm, 1.2 mm, 1.0 mm, 0.8 mm and 0.6 mm respectively.

From the numerical results, there’re three periods for a typical bubble behavior before inclined solidified front, which is floating up, flowing along the solidified front and entrapment or repulsion. A concept of entrapment duration is created to indicate the duration that bubble flows along the solidified front before starting to be entrapped. As shown in Fig. 6, for bubble whose diameter is smaller than 1.2 mm, its entrapment duration is increasing as its diameter increasing. But for those with diameter bigger than 1.2 mm, it will be repulsed by the advancing solidified front. The specified period contrast for different each bubble is listed in Table 5. It can be concluded that the entrapment possibility decreases as bubble diameter increasing. It is indicated that a critical diameter number exist for bubble entrapment. As is shown in Fig. 6 and Table 5, 1.2 mm is considered as the critical number for bubble entrapment in this condition. And this critical diameter number is under the effects of inclined angle of 50° and cooling rate of 3000 W/(m²·K).

Fig. 6.

Distribution of bubble behavior and entrapment duration with the varying diameter.

Table 5.Period duration of interaction between bubble and solidification front.

	Contact	Entrapment	Fully Wrapped	Entrapment Duration
1.8 mm	0.050 s	—	—	—
1.6 mm	0.052 s	—	—	—
1.4 mm	0.055 s	—	—	—
1.2 mm	0.060 s	0.161 s	1.7 s	0.099 s
1.0 mm	0.070 s	0.13 s	1.3 s	0.060 s
0.8 mm	0.082 s	0.15 s	1.2 s	0.058 s
0.6 mm	0.100 s	0.13 s	1.2 s	0.030 s

Meanwhile, two bubbles with diameter of 2.0 mm and 1.0 mm are chosen to demonstrate and analyze the typical behavior of repulsion/entrapment before advancing solidified front. In Fig. 7, the bubble with diameter of 2.0 mm floats up from initial position and touches the solidified front at t = 0.05 s. But it is repulsed and flicked away after a short contact with solidification front and keeps floating up to the upper surface. During its interaction with solidified front and floating up, the deformation of bubble is observed. Due to its small Eötvös number (Eo) and lower Morton number (Mo), shape of the single rising bubble keeps ellipsoidal cap, which is in a good agreement with Bhaga and Weber’s bubble shape regime diagram.¹⁸⁾ Furthermore, during the interaction with the solidified front, the bubble keeps the shape of ellipsoidal cap with the contacting corner much sharper as shown in Fig. 7(c).

Fig. 7.

Behavior of bubble with 2.0 mm diameter before advancing solidified front.

As shown in Fig. 8, the bubble with a diameter of 1.0 mm floats up from initial position and touches the solidified front at t = 0.07 s. Subsequently, it slowly flows upon the solidified front with its speed decreasing gradually. At t = 0.13 s, the bubble halts and starts to be entrapped by the advancing solidified front. It’s fully wrapped up at t = 1.3 s. During the entrapment, due to the high thermal conductivity ratio between water and air, thermal gradient around gas bubble leads to a deformation in the solidified front. As shown in Figs. 8(d) to 8(f), a bulging out is forming in the solidified front as well as the solidification isothermal. This phenomenon is in good agreement with the theoretical study of Hadji⁸⁾ and the X-ray transmission microscopy investigation by S. Sen.¹⁹⁾

Fig. 8.

Behavior of bubble with 1.0 mm diameter before advancing solidified front.

In order to validate the numerical model and its results, the numerical results are compared with the corresponding experimental results. From the numerical model and calculation results, the conditions for the bubble entrapment is that as it touches the solidified front, the balance among buoyancy force, gravity and solidified front resistant force makes it stay still at the solidified front and start to be entrapped by the advancing solidified front gradually. Therefore, for the bubble with bigger diameter, a corresponding higher buoyancy force will drive it floating up along the solidified front. Accordingly, the result is opposite for smaller bubble. Correspondingly, the experimental result is the same as the numerical one. The force balance keeps the bubble stay still at the solidified front and bigger diameter makes it harder for entrapment. And for small bubbles, this phenomenon is for the reason that it’s small enough to enter the spacing between dendrite arms and get entrapped by advancing solidified front. But this explanation for numerical results requires a more detailed numerical model involving dendrite growing and force balance of bubble dependent on dendrite arm space, which will be our future work.

As a result, the numerical simulation model is validated to be authentic based on the comparison with experimental results. Meanwhile, compared with Esaka’s work²⁾ mentioned in the introduction, our numerical results lead to a good agreement and another validation for our numerical model. In addition, our work gives a quantitative result for bubble behavior rather than a qualitative one obtained from the in-situ experiment. So, it’s believed that such a methodology based on Level Set method and modified solidification equation can be used in the simulation of bubble behavior before solidified front.

3.2.2. Effects of Cooling Rate and Inclined Angle on Bubble Behavior

Meanwhile, the effects of cooling rate and inclined angle of solidified front on bubble behavior are also investigated numerically. The cooling rate’s effect is studied through the control of outward heat transfer coefficient h at the inclined solidified boundary with the inclined angle being 45° and bubble diameter being 1 mm. The values of heat transfer coefficient are taking 1000, 2000, 3000, 4000 and 5000 W/(m²·K). Then, inclined angles ranging from 40° to 60° are chosen to observe the bubble behavior with bubble diameter being 1 mm and outwards heat transfer coefficient being 3000 W/(m²·K). The concept of entrapment duration is used to estimate the difficulty for a bubble to be entrapped by advancing solidified front.

The numerical results of cooling rate’s effect which is also known as heat transfer coefficient on bubble behavior given in Fig. 9(a) demonstrate that as cooling rate increasing, less time is needed for a bubble with diameter of 1.0 mm to get entrapped by advancing solidified front. As for the inclined angle of solidified front, 5 different values are taken to investigate its influence on bubble behavior. For inclined angle of 55° and 60°, bubble floats up along the solidified front after contacting with it and flicks away, resulting a repulsion end. Yet, for inclined angle of 40°, 45° and 50°, bubble is successfully entrapped by the advancing solidified front. The entrapment duration comparison for these three angles are shown in Fig. 9(b), which demonstrates that as the inclined angle of solidified front increasing, it leads to a higher entrapment possibility and less time for entrapment duration less time. And the inclined angle of 55° is considered as the critical number for bubble entrapment in this condition. In general, lower inclined angle of solidified front and higher cooling rate at solidified boundary results in a bigger possibility of bubble entrapment, which apparently agrees with work by Esaka.²⁾

Fig. 9.

Distribution of bubble behavior and entrapment duration (a) with the varying heat transfer coefficient, (b) with the varying inclined angle.

4. Discussion

Generally, a numerical method aiming at studying the bubble behavior before solidified front together with experimental validation and analysis is built. The innovation and advancement of this work lies in studying the interaction between gas bubbles and advancing solidified front with an efficient numerical simulation method instead of experimental method or physical model analysis, which leads to a direct and intuitive result to analyze. The simulation results and its comparison with experimental results lead to a consistent result and agreement with traditional theory. In our future work, we plan to improve and perfect our observation method of bubble behavior before solidified shell with the help of art X-ray transmission microscope (XTM)²⁰⁾ or real-time X-ray video microscopy.²¹⁾

However, according to previous researchers’ work, there’s a lot of work left to improve the numerical model in the next stage. In order to obtain the precise value of critical point for bubble entrapment/repulsion, a series of numerical models are needed to gradually narrow down its range. The interaction taken into account in this study is the resistant force effect of advancing solidified front on bubble and the buoyancy force of itself. In the future plans of work, the thermocapillary force dependent on thermal gradient, the unsteady forces due to acceleration of the bubble with respect to the fluid, added mass force and Basset force, Magnus force that is due to the spinning bubble which are acting on bubble will be taken into account.

5. Summary

A finite element method numerical model is presented to simulate the bubble behavior before the advancing solidified front. This methodology is a highly nonlinear coupling mathematic model composed of level set method for two-phase flow interface evolution and a modified heat transfer model for solidified process. Both the solidus interface and gas-liquid interface evolution and deformation are obtained by these two physical fields. The validation for this numerical model is made successfully through comparison between experimental and numerical results of bubble diameters’ effects on its behavior.

In that circumstance, the effects of bubble diameter, inclined angle and cooling rate on bubble behavior before advancing solidified front are investigated numerically and experimentally. Each critical value for bubble entrapment/repulsion is obtained by testing a series of numerical models. Conclusions are made as below:

(1) Due to the high thermal conductivity ratio between water and air, bulging out is observed in the solidified front when bubble is approaching the solidified front.

(2) Entrapment of bubble deceases with increasing diameter of bubble. With the inclined angle of solidified front being 45° and heat transfer coefficient of 3000 W/(m²·K), the critical diameter number for entrapment is 1.2 mm.

(3) Entrapment of bubble deceases with decreasing cooling rate.

(4) Entrapment of bubble decreases with increasing inclined angle. The critical inclined angle number for entrapment is 55° for bubble diameter of 1.0 mm and heat transfer coefficient is 3000 W/(m²·K).

Acknowledgement

The authors wish to express their appreciation to National Natural Science Foundation of China for providing financial support which enabled this study to be successfully carried out [Grant No. 50874133 and No.51274137].

References

• 1)   L.  Zhang and  S.  Taniguchi: Int. Mater. Rev., 45 (2000), 59.
• 2)   H.  Esaka,  Y.  Kuroda,  K.  Shinozuka and  M.  Tamura: ISIJ Int., 44 (2004), 682.
• 3)   B. G.  Thomas and  L.  Zhang: ISIJ Int., 41 (2001), 1181.
• 4)   M. S.  Park,  A. A.  Golovin and  S. H.  Davis: J. Fluid Mech., 560 (2006), 415.
• 5)   Z.  Wang,  K.  Mukai and  I. J.  Lee: ISIJ Int., 39 (1999), 553.
• 6)   A. S.  Dzyuba and  Z. Y.  Yong: Soviet Physics - Crystallography, 30 (1985), 684.
• 7)   H.  Lei,  M.  Zhu and  J.  He: Int. J. Cast Metal. Res., 20 (2007), 98.
• 8)   L.  Hadji: Phys. Rev. E, 75 (2007), 042602.
• 9)   M.  Pasandideh-Fard,  R.  Bhola,  S.  Chandra and  J.  Mostaghimi: Int. J. Heat Mass Transf., 41 (1998), 2929.
• 10)   M.  Pasandideh-Fard,  S.  Aziz,  S.  Chandra and  J.  Mostaghimi: Int. J. Heat Fluid Flow, 22 (2001), 201.
• 11)   M.  Bussmann,  J.  Mostaghimi and  S.  Chandra: Phys. Fluids, 11 (1999), 1406.
• 12)   M.  Sussman,  A. S.  Almgren,  J. B.  Bell,  P.  Colella,  L. H.  Howell and  M. L.  Welcome: J. Comput. Phys., 148 (1999), 81.
• 13)   J. A.  Sethian: Proc. Natl. Acad. Sci. USA, 93 (1996), 1591.
• 14)   S.  Osher and  J. A.  Sethian: J. Comput. Phys., 79 (1988), 12.
• 15)   V.  Voller and  C.  Prakash: Int. J. Heat Mass Transf., 30 (1987), 1709.
• 16)   J. J.  Bikerman: Physical Surfaces, Academic Press, New York, (1970), 20.
• 17)   E.  Takeuchi and  J.  Brimacombe: Metall. Mater. Trans. B, 15 (1984),
• 18)   D.  Bhaga and  M.  Weber: J. Fluid Mech., 105 (1981), 61.
• 19)   S.  Sen,  P.  Curreri,  W.  Kaukler and  D.  Stefanescu: Metall.Mater. Trans. A, 28 (1997), 2129.
• 20)   A. V.  Catalina,  D. M.  Stefanescu,  S.  Sen and  W. F.  Kaukler: Metall. Mater. Trans., 35A (2004), 1525.
• 21)   L.  Arnberg and  R. H.  Mathiesen: JOM, 59 (2007), 20.