An Analysis of Bubble Migration in Horizontal Thermo-Capillarity Using the VOF Modeling

Abstract

Due to various engineering applications, spontaneous bubble movement on the heated surface has brought huge attention. This work numerically studied the bubble migration driven by the thermo-capillary force under the temperature gradients perpendicular to the gravity direction. This problem is constructed in a two-dimensional domain, and the volume of fluid (VOF) method is adopted to capture the properties of the bubble interface between the vapor and the liquid. One still vapor bubble is initially positioned at the center of the liquid domain, and the temperature gradient is applied to two side walls. The results show that the bubble with a size greater than the capillary length can only oscillate near the initial position even with a larger temperature gradient. The deformation of the bubble such as spheroid and spherical cap can be found around this regime. However, the movement of the bubble with a size smaller than the capillary length is significant under a higher temperature gradient, and it remains a spherical shape. The coefficient of thermo-capillary force (C_Th) is defined within this work, and it is found that a larger Weber number ($W_e$) accomplishes a larger C_Th. This work may provide more precise guidance for smart bubble manipulation and critical heat flux estimation for future nuclear reactor design.

1. Introduction

Bubbles and drops can be found everywhere, and they also appear in various applications in engineering and material processing. To actively manipulate the bubble migration is essential and with a sound potential. Surface tension has been considered constant in most cases under the mesoscale, while on the microscale, surface tension as a function of temperature progressively decreases with increasing temperature, which cannot be ignored. Due to this temperature gradient, the non-uniform distribution of surface tension causes the thermocapillary effect. The non-uniform surface tension at the fluid interface contributes to viscous forces acting on the outer fluid, causing fluid particle motion in the thermal gradient direction. The bubbles suspended in a temperature gradient fluid will then move towards a hot zone. Thermo-capillary force can be crucial in the processing of materials in a reduced gravity environment, such as in a space shuttle as well as in separation processes for recycling life-sustaining water and oxygen in space travel.

The studies of bubble/droplet migration due to thermo-capillarity can be tracked back to 1959 by Young et al. [1]. They first derived an analytical solution of the bubble migration velocity under thermo-capillary and buoyant forces. Balasubramaniam et al. investigated the bubble movement for a small Marangoni number [2]. They also studied the bubble migration at reduced gravity with different Marangoni numbers [3,4]. Nas and Tryggvason showed the numerical simulations of the thermo-capillary motion of a pair of two- and three-dimensional fully deformable bubbles and drops [5]. Zhang et al. numerically investigated the bubble/droplet migration due to thermo-capillary force under uniform thermal radiation with the level set method [6]. Oliver and DeWitt demonstrated the analytical solution of surface tension-driven flows for a droplet in a micro-gravity environment [7]. Ma et al. numerically investigated the bubble migration driven by the thermo-capillary effect (temperature gradient parallel to the gravity direction) under variable buoyancy [8]. Haj-Hariri et al. showed the numerical simulations of the 3D thermo-capillary motion of deformable viscous drops under the influence of a constant temperature gradient within the liquid medium [9]. Using a finite difference method, Chang analyzed the effect of bubble deformation while predicting thermo-capillary migration [10]. They found that decreased surface tension contributes to a higher deformation rate at higher temperatures. Hadland et al. showed the results from the experiments on the thermo-capillary migration of air bubbles and fluorinert drops in a Dow-Corning silicone oil aboard a NASA Space Shuttle mission [11]. The capillary-driven force also benefits the development of high-performance heat transfer devices, such as microscale heat pipes [12].

Other than spontaneously manipulating the bubble movement, bubble behavior is also important to determine the heat flux. In a nuclear reactor, boiling can be found on the heating surface of the fuel rods. For safety, knowing the critical heat flux (CHF) of the fuel rods is vital for designing a nuclear reactor. Bubbles generated on the fuel rods will depart from the heating surface and migrate into the flow channel under lateral forces such as wall lubrication, dispersion, thermo-capillarity, etc. It is intuitive to imagine a more significant void fraction near the heating surface due to the heat and mass transfer. The CHF will decrease since the conductivity of vapor is smaller than liquid. Understanding the bubble dynamics near the heated surface is essential to determine the critical heat flux. Previous research mostly focused on the thermo-capillarity parallel to the gravity direction. However, this work investigates how a single bubble migrates under the temperature gradient perpendicular to gravity using the volume of fluid method (VOF). The primary goal of this work is to get a better understanding of the physics in a thermo-capillary bubble to pursue more precise critical heat flux estimation in the near future.

2. Governing Equations

Various methodologies have been proposed to model the vapor-liquid two-phase problems, such as the molecular dynamic (MD) simulation [13], lattice Boltzmann method [14,15], immersed boundary methods [16], level set (LS) [17] methods, and volume-of-fluid method [18,19]. The VOF method proposed by Hirt and Nichols in 1981 tracks the interface between two phases [20], and it has been used widely in studying various two-phase flow systems. The VOF method is adopted here since the fluid mass can be appropriately conserved, and it can be applied on a larger scale compared with the LBM and MD methods. Also, it’s computationally and timely cheap. Therefore, the VOF model is chosen in this work to characterize the liquid/vapor interface to illustrate the bubble dynamics. In the VOF model, α is the volume fraction of two phases with subscripts υ and l denoting vapor and liquid. The volume fractions of the two phases in the VOF model can be defined as:

$${\alpha }_{\upsilon }+{\alpha }_l=1,$$

(1)

The continuity equation describes the conservation of mass as:

$$\{\begin{array}{c}\frac{\partial }{\partial t}\left({\alpha }_{\upsilon }{\rho }_{\upsilon }\right)+\nabla ·\left(u{\alpha }_{\upsilon }{\rho }_{\upsilon }\right)=0,\\ \frac{\partial }{\partial t}\left({\alpha }_l{\rho }_l\right)+\nabla ·\left(u{\alpha }_l{\rho }_l\right)=0,\end{array}$$

(2)

where $t$ is time and u is the average fluid velocity. The momentum equation is expressed as:

$$\frac{\partial }{\partial t}\left(\rho u\right)+\nabla ·\left(\rho uu\right)=-\nabla p+\nabla ·\left(\mu \nabla u\right)+\rho g+F,$$

(3)

where p is the pressure, g is the gravitational acceleration, and F is the external force. The non-dimensional fluid properties of the mixture density $\rho$, dynamic viscosity μ, specific heat $C_p$, and thermal conductivity k can be expressed as:

$$\{\begin{array}{c}\begin{array}{c}\rho ={\rho }_l{\alpha }_l+{\rho }_v{\alpha }_v,\\ \mu ={\mu }_l{\alpha }_l+{\mu }_v{\alpha }_v,\end{array}\\ C_p=C_{pl}{\alpha }_l+C_{pv}{\alpha }_v\\ k=k_l{\alpha }_l+k_v{\alpha }_v\end{array}$$

(4)

The conservation of energy is expressed as:

$$\frac{\partial }{\partial t}\left(\rho C_{p}T\right)+\nabla ·\left(\rho u{C}_{p}T\right)=-\nabla ·\left(k\nabla T\right)+S_h$$

(5)

where T, $C_p$, k, and $S_h$ are the temperature, the specific heat, the mixture thermal conductivity, and the heat source, respectively. Brackbill [21] proposed the continuum surface force (CSF) model, which assumes that surface tension is constant at the same temperature, where the forces are typical to the two-phase interface. In the CSF model, the surface tension is considered the source term in the momentum equation in the VOF model. With the surface tension σ and two radii orthogonal to the interface $R_1$ and $R_2$, the equation for the pressure drop, $p_2$ − $p_1$, is expressed as:

$$p_2-p_1=\sigma \left(\frac{1}{R_1}+\frac{1}{R_2}\right)$$

(6)

In previous research, surface tension σ is considered a constant. However, the surface tension σ is a function of temperature in this work. Hence, in the momentum equation, F is the surface tension force between the two phases that are expressed as a volume force density:

$$F=\sigma \left(T\right)\frac{{\alpha }_l{\rho }_l{\kappa }_l\nabla {\alpha }_l+{\alpha }_{\upsilon }{\rho }_{\upsilon }{\kappa }_{\upsilon }\nabla {\alpha }_{\upsilon }}{\frac{1}{2}\left({\rho }_l+{\rho }_{\upsilon }\right)},$$

(7)

where κ is the interface curvature denoted as:

$${\kappa }_l=-{\kappa }_{\upsilon }=-\nabla ·\left(\frac{\nabla {\alpha }_l}{\left|\nabla {\alpha }_l\right|}\right).$$

(8)

3. Problem Definition and Numerical Method

In Figure 1, the schematic diagram illustrates the problem definition and boundary conditions in this work. The computational domain is a 0.12 m × 0.04 m rectangle, and the time step size is chosen as $\Delta t=2.5\times {10}^{-4}$ s. One spherical vapor bubble is initially positioned near the bottom wall at the center region of the liquid. The distance of the bubble centroid from the bottom wall is R_bubble + 0.005 m. All solid walls are set to the no-slip boundary conditions. The surface tension is defined as a function of temperature σ = f(T). The temperature of the left wall is fixed at T_L = 373 K, and that of the right wall is varied from T_R = 273 K to 373 K. The computations are performed for two-dimensional, unsteady, and incompressible flow with no mass change. The governing equations are solved using the finite volume method performed on ANSYS Fluent 2020 R1. The PISO scheme is chosen for pressure-velocity coupling. The 2nd order upwind scheme is used for the conservation of both momentum and energy. The PRESTO! and the Geo-Reconstruct schemes are adopted for the pressure and the volume fraction interpolation, respectively. The Courant number is set at 0.25. The solution is converged well below the residuals of ${10}^{-5}$. The numerical simulation settings are listed in

Table 1. Numerical simulation settings.

Item	Content
Multiphase model	VOF explicit
The surface tension force model	Continuum Surface Force
Viscous model	Laminar
Pressure-velocity coupling	PISO
Momentum and energy	2nd order upwind
Volume fraction	Geo-Reconstruct
Maximum iteration	100

.

To ensure numerical stability, the density ratio of liquid to vapor can’t be large in the VOF method [14]. The working fluid properties used in the simulations are listed in

Table 2. Properties of the working fluid in simulations.

Parameter	Liquid	Vapor
$\rho$ Density ($\mathrm{kg}/{\mathrm{m}}^3$)	$2.0\times {10}^2$	$5.0\times {10}^0$
$\mu$ Viscosity ($\mathrm{Pa}·\mathrm{s}$)	$1.0\times {10}^{-1}$	$5.0\times {10}^{-3}$
$k$ Thermal conductivity ($\mathrm{W}/\mathrm{m}·\mathrm{K}$)	$4.0\times {10}^1$	1.0
$C_p$ Specific heat ($\mathrm{J}/\mathrm{kg}·\mathrm{K}$)	$4.0\times {10}^2$	$2.0\times {10}^2$
$h_{lv}$ Latent heat ($\mathrm{J}/\mathrm{kg}$)	$1.0\times {10}^4$
$\sigma$ Surface tension ($\mathrm{N}/\mathrm{m}$)	(−0.001 × Temperature + 0.473) × 15

. Further comparison between the properties of working fluid and water can be done based on non-dimensional parameters. The primary non-dimensional parameters are listed below:

$$Re=\frac{{\rho }_l\sqrt{g_y{\lambda }_0}{\lambda }_0}{{\mu }_l}, Pr=\frac{C_{pl}}{{\mu }_l{k}_l}, We=\frac{g_y{\lambda }_0{}^2{\rho }_l}{\sigma }, Pe=\frac{C_{pl}{\lambda }_0\sqrt{g_y{\lambda }_0}{\rho }_l}{k_l}, Ja=\frac{C_{pl}\Delta T}{h_{lv}}.$$

where ${\lambda }_0$ is the capillary length, ${\rho }_l$ is the liquid density, $g_y$ is the gravity acceleration, ${\mu }_l$ is the liquid viscosity, $k_l$ is the liquid conductivity, $C_{pl}$ is the liquid specific heat, $\sigma$ is the surface tension, $h_{lv}$ is the latent heat, and the $\Delta T$ is the temperature difference. The calculated non-dimensional settings are $R_e=3.85$, $P_r=1$, $W_e=1.023$, $P_e=3.849$, and $J_a=8\times {10}^{-2} \mathrm{to} 4$. The Jakob number ($J_a$) of the working fluid in the simulations is larger than that of real water, which represents a better ability to phase change. The Peclet number ($P_e$) of the working fluid is smaller than that of water, representing lower convection. The Reynolds numbers ($R_e$) in both fluids are chosen in the laminar regime. The Weber numbers ($W_e$), which depend on surface tension, are set close to unity in both fluids. The Prandtl numbers ($P_r$) in both fluids are within the same order of magnitude. The parameters chosen for the working fluid in our simulations are based on previous work [22,23], as listed in Table 2. The only difference is the surface tension now is chosen according to the temperature and is in the same order as the previous one.

4. Results and Discussion

This section discusses the model verification and simulation results with varying temperature gradients and bubble radius sizes.

4.1. Numerical Verification

First, the numerical results of our model are benchmarked against the theoretical solution for verification. The capillary length $L_c$ is essential for the surface tension effect and is defined as follows:

$$L_c=\sqrt{\frac{\sigma }{g\left({\rho }_l-{\rho }_v\right)} }=3.23\times {10}^{-2} \mathrm{m}$$

(9)

where $\sigma$ is the surface tension, $g$ is the gravity acceleration, and ${\rho }_l$ and ${\rho }_v$ represent the density of liquid and vapor. The theoretical thermos-capillary migration velocity $V$ of a bubble/drop in a vertical temperature gradient was first given by Young [1]

$$V=\frac{2}{3}\left[\mu {\gamma }^{\prime }R{T}_C-\left(\rho -{\rho }^{\prime }\right)g{R}^2\left(\mu +{\mu }^{\prime }\right)\right]{\left(3\mu +2{\mu }^{\prime }\right)}^{-1},$$

(10)

where R is the bubble radius, $T_c=3T_1/\left(2+(h^{\prime }/h\right))$, $\rho$ and ${\rho }^{\prime }$ are the densities and viscosities of fluid outside the bubble, respectively, ${\rho }^{\prime }$ and ${\mu }^{\prime }$ are the densities and viscosities inside the bubble, and ${\gamma }^{\prime }$ is the temperature coefficient of surface tension. Balasubramaniam also performed experiments on bubble migration under thermo-capillary force aboard the NASA Space shuttle in orbit [2]. The thermo-capillary migration velocity of a drop V in reduced gravity is

$$V=-\frac{2\left|\nabla T_{\infty }\right|{\sigma }_{T}R}{\mu \left(2+3\alpha \right)\left(2+\beta \right)} ,$$

(11)

where $\nabla T_{\infty }$ is the temperature gradient,${\sigma }_T$ is the rate of the change of the interfacial tension with temperature, $\alpha$ is the viscosity ratio of inside and outside the bubble, and $\beta$ is the thermal conductivity ratio of inside and outside the bubble. The size of the computational model is 0.12 m × 0.04 m (See Figure 2). Initially, the vapor bubble with a radius R = 0.005 m is positioned near the bottom wall. The gradient of temperature between two walls $\nabla T$ is set as 833.3 K/m. Using Equation (11), the migration velocity is calculated as V = 0.15625 m/s. Figure 2 shows the transient bubble position with zero gravity. As time t increases, it is found that the temperature contour tends to show linear distribution, and the bubble moves toward the hot (left boundary) wall and remains at the same height under no gravity. Our simulation result shows that the migration velocity is approximately 0.142 m/s which is smaller than the theoretical solution, but it falls in the same order which shows a qualitative agreement. However, Balasubramaniam’s experiment [2] also showed a smaller migration velocity compared with the theoretical results. According to YGB’s (Young, Goldstein, and Block) model [1], the deformation of the bubble requires large enough kinetic energy, and the bubble deformation is not seen with $V_{real}<V_{YGB}$, which is in good agreement.

4.2. Bubble Sizes from R = 0.01 m to 0.05 m

Since the capillary length is found at $R=3.23\times {10}^{-2} \mathrm{m}$, two groups of bubbles will be categorized for discussion. The bubble sizes of the first group are chosen from R = 0.01 m to 0.05 m, i.e., R = 0.01 m, 0.015 m, 0.025 m, 0.035 m, and 0.05 m. The temperature of the right wall T_R is ranged from 273 K to 373 K with a constant left wall temperature T_L = 373 K. The effects of the temperature gradient, bubble radius, and bubble shape will be discussed as follow.

4.2.1. Effect of Temperature Gradients

Figure 3a–f show the bubble migration under gravity for bubble sizes of R = 0.01 m, 0.015 m, 0.025 m, 0.035 m, and 0.05 m under various temperature gradients of 0 K, 20 K, 40 K, 60 K, 80 K, and 100 K, respectively. Since the temperature gradient is applied to the direction perpendicular to the gravitational force, we are only interested in the bubble movement in the x-direction. The x-axis represents the centroid bubble position, and the y-axis represents the time. The surface tension is more significant in a lower temperature than higher temperature, hence, a net force on the vapor bubble is induced to move the bubble. If there is no temperature difference $\Delta T=0 \mathrm{K}$, as shown in Figure 3a, the bubbles with different sizes almost stay near the center (x = 0.06 m). With a larger temperature gradient, it is much easier for the vapor bubble to reach the left wall in a shorter time. It is also noted that the bubble with a smaller size arrives at the left wall much faster than a larger bubble under the same temperature gradient. The bubble with a radius of 0.05 m only oscillates around its initial position, which is the center of the domain. We believe that the bubbles oscillate because of the Kelvin–Helmholtz instability. The Kelvin–Helmholtz instability is due to the perturbation of the velocity field, and the larger bubble will induce larger perturbation and hence larger oscillation.

4.2.2. Effect of Bubble sizes

Figure 4a–e show the centroid positions of the bubbles with different radii. When the temperature of the right wall decreases, the bubble tends to move toward the left wall, as shown in Figure 4a–d. With a radius R = 0.01 m or 0.015 m, the bubble tends to move faster to the left wall, except for no temperature gradient applied. The medium-sized bubble with a radius R = 0.025 m or 0.035 m tends to move gradually towards the left wall, and the biggest bubble, R = 0.05 m, prefers to stay near the center, as shown in Figure 4e. Therefore, the bubble size plays a significant role in the migration dynamics. Furthermore, a smaller bubble does not change its shape, so its migration velocity is faster. However, a larger bubble changes its shape from sphere to spheroid or spherical cap, reducing its flow velocity and moving slowly.

4.2.3. Bubble Shape Deformation and Rising Velocity

The vapor bubbles can be classified into three different shapes: sphere, spheroid, and spherical cap. The bubble is spherical when its inertial force is much lower than surface tension or viscous force. As bubble size increases, velocity rises since the larger bubble will induce larger buoyancy and hence higher oscillation. The shape of the bubble is transforming the drop from a sphere into oblate spheroid shapes. Due to the resistance in the liquid medium, when bubbles are large enough, they tend to have flat and often curved bases, breaking up the curve of the bubble shape. This is known as the spherical cap. Figure 5 shows the rising velocity and bubble shape for different bubble radii. The computational domain of Figure 5 is 0.12 m $\times$ 0.3 m to perform a better picture of the bubble rising. When the bubble size R is 0.01 m and 0.015 m, the bubble shape remains spherical. The bubble shows a spheroid shape with a radius R = 0.025 m. The bubble shape can be found in a spherical cap with a radius R = 0.035 m and 0.05 m. It is noted that the deformation of the bubble occurs near the capillary length (R < L_c $\approx 3.23\times {10}^{-2}$ m). Figure 6 shows both our simulation results (in red lines) and previous experiment work (in black lines) [24]. Since the fluid properties in the two works are different, the magnitudes of bubble size and velocity are different, as expected. However, our numerical work shows a similar trend to experiment one qualitatively.

4.3. Bubble Sizes from R = 0.0045–0.01 m

In this section, the bubble size is decreased to R = 0.0045 m to 0.01 m. The bubble sizes are R = 0.0045 m, 0.005 m, 0.006 m, 0.007 m, 0.008 m, and 0.01 m, respectively. The temperature of the right wall T_R is ranged from 273 K to 373 K, and the left wall remains the same as T_L = 373 K. The effects of the temperature gradient, bubble radius, and shape are discussed as follows.

4.3.1. Effect of Temperature Gradients

Figure 7a–f show the movement of the centroid position of the bubble in the x-direction with a radius from R = 0.0045–0.01 m under different temperature gradients, 0 K, 20 K, 40 K, 60 K, 80 K, and 100 K, respectively. It is found that the bubble stays at the center under the uniform temperature, as shown in Figure 7a. Once the temperature difference $\Delta T$ between the two walls increases, the vapor bubble starts to migrate slowly towards the hotter side, as shown in Figure 7b–d. Compared with the first group (R = 0.01 m to 0.05 m), it takes more time for smaller bubbles to reach the left wall. At a higher temperature gradient, $\Delta T=$ 80–100 K, it is much easier for the bubble to arrive at the left wall in a shorter time, as shown in Figure 7e–f.

4.3.2. Effect of Bubble Sizes

Figure 8 shows the transient bubble centroid position in the $\mathrm{x}$-axis with varying bubble radius from 0.0045 m to 0.01 m. In Figure 8a, when the bubble radius R is 0.0045 m, the centroid position of the bubble does not move for the cases of $\Delta T$= 0 and 20 K. With a larger $\Delta T$ > 40 K, the bubble moves faster. Figure 8b–f show that the vapor bubbles remain at the center in the x–direction when $\Delta T = 0 \mathrm{K}$. Once the temperature difference $\Delta T$ between two walls increases, the bubbles start to move toward the left wall. It is shown in Figure 8 that the bubble with a larger size can reach the left wall faster than a smaller one. It is because the bubble with a larger size will encounter a larger temperature difference and hence higher surface thermal capillary force. Therefore, it reaches the left wall faster. Compared with the first group of larger bubbles (R = 0.01 m to 0.05 m), the bubbles with smaller radii remain spherical in shape for all sizes of bubbles from 0.0045 m to 0.01 m. Figure 9 shows the transient state bubble with radius R = 0.006 m and ΔT = 40 K at different time instants in (a) velocity field, (b) temperature field, and (c) pressure field.

4.4. Force Analysis

In fluid mechanics, it is important to correlate parameters to meaningful dimensionless numbers. It is well known that the drag force is one of the four primary forces acting on an object moving through a fluid. In this way, the physics can be observed and discussed on a scaled-down model, e.g., an airfoil model in the wind tunnel, as long as the Re number is kept the same. In this work, we are goanna plot $C_{Th}$ against the Weber number (We) because this information might be helpful for future development of the correlation of $C_{Th}$ as a function of the We number.

The drag coefficient of $C_D$ is defined as

$$C_D=\frac{F_D}{\frac{1}{2}{\rho }_l{V}^{2}A}$$

(12)

where $F_D$ is the drag force, ${\rho }_l$ is the density of the liquid, $V$ is the velocity, and $A$ is the reference area. $C_D$ can be represented as a function of dimensionless parameters such as Reynolds number, Froude number, Mach number, and surface roughness. That is

$$C_D=f\left(\mathrm{shape}, \mathrm{Re}, \mathrm{Fr}, \mathrm{Ma},\dots \right).$$

(13)

In this work, we are interested in finding the relationship between the coefficient of lateral force ($C_{Th}$) and the Weber number since the thermo-capillary force is dominant in this system. Weber number is defined as:

$$W_e=\frac{2{\rho }_l{V}^{2}R}{\sigma },$$

(14)

where $R$ is the bubble radius, and $\sigma$ is the surface tension. Hence, to determine the lateral force acting on the vapor bubble, the following assumptions have been made: (1) there is no phase change between the liquid and vapor, and (2) the dominant forces applied to the bubble in the $x$-direction are the drag force and thermo-capillary force (lubrication force due to the wall is neglected).

$$\Sigma F_x=F_D+F_{La}=\frac{1}{2}{C}_D\rho V^{2}A+C_{th}\nabla \sigma A=m{a}_x.$$

(15)

This is the force balance equation in the $x$-direction, where the $C_D$ represents the drag coefficient, V is the bubble velocity, A is the projection area, and $C_{Th}$ is the lateral force coefficient due to the thermo-capillary force, and ∇σ is the surface tension gradient. For a two-dimensional problem, the bubble mass is $\rho \pi R^2$, and the projection area is 2 R. By substituting the above quantities into the force balance Equation (15), the coefficient of lateral force $C_{Th}$ is then derived as:

$$C_{Th}=\frac{F_{La}}{\nabla \sigma A}=\frac{m{a}_x-F_D}{\nabla \sigma A}=\frac{m{a}_x-\frac{1}{2}{C}_D\rho V^{2}A}{\nabla \sigma A}.$$

(16)

Figure 10 shows the relationship between the thermo-capillary coefficient (C_Th) and the Weber number (We). The bubble sizes range from $R=0.0045 \mathrm{m}-0.03 \mathrm{m}$, which are smaller than the Capillary length. The Weber numbers then fall between ${10}^{-4}-{10}^0$. It is found that as the bubble radius increases, the Weber number increases. Since a larger bubble size induces a higher drag force, the surface tension effect is not significant. In addition, the thermo-capillary coefficient is large due to a large thermo-capillary force.

For bubble sizes in the range of $0.0045 \mathrm{m}<R\le 0.007 \mathrm{m}$, Figure 11 shows the thermo-capillary coefficient ($C_{Th}$) vs. the Weber number (We) for R = 0.0045 m. It is noted that both the Weber number and thermo-capillary coefficient increase when ΔT increases.

For bubble sizes in the range of $0.008 \mathrm{m}\le R<$0.01 m, Figure 12 shows the thermo-capillary coefficient (C_Th) vs. the Weber number (We) for R = 0.008 m. It is found that the temperature difference $\Delta T$ does not affect the Weber number much. However, the thermo-capillary coefficient decreases when $\Delta T$ increases except for $\Delta T=20 \mathrm{K}$.

Also, for the bubble sizes in the range of $0.01 \mathrm{m}\le R\le 0.03 \mathrm{m}$, the temperature difference $\Delta T$ does not change the Weber number much (all around ~1). Figure 13 shows the thermo-capillary coefficient (C_Th) vs. the Weber number (We) for R = 0.03 m. When $\Delta T=20-60 \mathrm{K}$, both the Weber number and thermo-capillary coefficient increase when $\Delta T$ increases. For a larger $\Delta T=80 \mathrm{K} \& 100 \mathrm{K}$, both Weber number and thermo-capillary coefficient decrease when $\Delta T$ increases.

5. Conclusions

In this work, the thermo-capillary effect and related dynamics of vapor bubble migration under different temperature gradients, bubble sizes, and bubble shapes have been studied using VOF numerical modeling.

• The simulation based on the volume-of-fluid method and temperature-dependent surface tension has been carried out. The results have been compared with theoretical and experimental studies, which demonstrates qualitative agreement. The corresponding dynamics have been qualitatively captured.
• The bubble shape deforms with the sizes larger than the capillary length, and the bubble shapes remain spherical with the sizes smaller than the capillary length.
• For the first group of bubble sizes from R = 0.01 m to 0.05 m, the bubble with radius R = 0.05 m remains near the center of the domain even under the temperature gradient ΔT = 100 K. The bubble with size R = 0.01 m can reach the left wall in a shorter time.
• For the second group of bubble sizes from R = 0.0045 m to 0.01 m, the larger size bubble can arrive at the left wall faster compared with the smaller size bubble.
• Coefficient of thermo-capillary force (C_Th) is defined. It is found that as the bubble radius increases, the Weber number increases since a larger bubble size induces a higher drag force and the surface tension effect is not significant. Also, the thermo-capillary coefficient increases due to a large thermo-capillary force.

It is verified that this approach can be used to investigate the thermo-capillary effect and provide more physical insights into bubble migration dynamics. In addition, it can provide sound guidance for precise critical heat flux estimation and direction of accurate manipulation of bubble migration which is of great importance in various engineering applications. In the future, we will extend our work to study thermo-capillary bubble migration, including solid surface wettability [25].