Transient behavior of the thermocapillary migration of drops under the influence of deformation
aa r X i v : . [ phy s i c s . f l u - dyn ] D ec Transient behavior of the thermocapillary migration of dropsunder the influence of deformation
L. Chang, Z. Yin, ∗ and W. Hu National Microgravity Laboratory, Institute of Mechanics,Chinese Academy of Sciences, Beijing 100190, P.R.China (Dated: July 4, 2018)
Abstract
The transient thermocapillary migration of drops with nontrivial deformation is studied. Thefinite difference method is employed to solve the incompressible Navier-Stokes equations coupledwith the energy equation; the front-tracking method is adopted to track the moving deformabledrop interface. In the hot region, deformations of drops increase with the decrease of interfacialtensions. In order to indicate the temperature impact on the interfacial tension, a local capillarynumber ( Ca l ) is introduced. It is found that, when the drop density is smaller/larger than that ofthe bulk fluid, the drop velocity decreases/increases with the increase of the drop deformation. ∗ Electronic address: [email protected] . INTRODUCTION Thermocapillary migration, which is induced by the variance of the interfacial tensionwhen a temperature gradient is imposed on the bulk fluid, is an important phenomenon indrop/bubble transportation. With the development of aeronautics and astronautics, it isnecessary to study the thermocapillary motion because of its practical roles in the materialprocessing, and the management of heat and fluids in space [1, 2]. In most studies on thethermocapillary motion, it is assumed that the interfacial tension linearly depends on thetemperature: σ ( T ) = σ + σ T ( T − T ) , (1)where the constant σ T is the interfacial tension coefficient of temperature, and σ the inter-facial tension at the reference temperature T .Considering the balance between the capillary force and viscosity on the bubble or drop,the reference velocity can be defined as U = | σ T ||∇ T ∞ | R /µ , in which |∇ T ∞ | is the temperature gradient imposed on bulk fluid (the subscript 1 repre-sented), µ the viscosity, and R the radius of the spherical drop/bubble. σ (10 − N/m ) σ T (10 − N/ ( m · K ))Gas bubble/Silicone oil[3] 17.3 -0.061Fc-75 drop/Silicone oil[3] 3.47 -0.036Monotectic alloy[4] 5.7 -0.23TABLE I: Typical interfacial tensions and their temperature coefficients. The capillary number is defined as Ca = µ U/σ = | σ T ||∇ T ∞ | R /σ . In most cases, the Ca numbers in thermocapillary migration are small (0.01 or less, e.g.,the materials in the first two lines of Table 1), and the deformations of drops are alsovery small. Therefore, previous studies usually assumed drops to be nondeformable [5].However, in practice, there are many situations where the interface tension is very small,2r where the deformation introduced by capillary effect is very large. For the case of thepreparation of monotectic alloys in micro-gravity conditions [6], the Ca number of the Al-based alloy ( Al . Bi . ) Sn is 0.12 at the temperature of 1200k if R = 300 µm and |∇ T ∞ | = 100 K/cm [7] (see the third line in Table 1). Such a big Ca number can causeobvious deformation on the drop. Similar situations happen in the chemical flooding methodsapplied in enhanced oil recovery [8–12]. Finally, experiments in space also observed someslight deformation of fairly large bubbles [3]. Therefore, it is necessary to investigate thedeformed drops with relatively large Ca ’s or small interfacial tensions in the thermocapillarystudies.When inertia and thermal convections are neglected, the velocity of the spherical drop inthermocapillary motion is presented in the pioneering work of Young et al. [13]: V Y GB = 2
AU, where A = 1(2 + 3 α )(2 + λ ) , (2)and α = µ /µ and λ = k /k are the ratios of the kinematic viscosity and thermal conduc-tivity between droplet (the subscript 2 represented) and background liquid, respectively.There are a lot of investigations studying the thermocapillary motion, which are reviewedin [14, 15]. Most of the researches sofar focus on the convection effects of inertia and energy,and only a few of them are about the drop/bubble deformation. When the thermal convec-tion is very small, the exact solution of the momentum equation was provided in [16], wheresmall inertial deformations of drops were also calculated. Under the same assumption, usingthe Lorentz reciprocal theorem, Haj-Hariri et al. calculated the small inertial deformationsas well as the consequent correction on the temperature field and the migration velocity[17]. Later, the work is extended to moderate parameters in a steady-state numerical inves-tigation [18]. It is concluded that even a small deformation will retard the motion of thedrop, and the steady-state migration can still be reached even for large capillary numbersup to 0.5. However, the actual continuous decrease of the interfacial tension when the dropmigrates to the hot region was not considered in their study. While taking this decrease intoconsideration, bubbles with large capillary numbers can not reach steady migrating states[19]. 3 r T ? q R q ( ) wallWallwall r z FIG. 1: The sketch of thermocapillary motion in the axisymmetric model.
According to our knowledge, transient migration behaviors due to drop deformation havenever been reported before, which is the main topic in this work. This paper is arrangedas follows: the governing equations and numerical methods are introduced in Section 2, thevalidation tests are presented in Section 3, and the results and discussions are in Section 4.
II. GOVERNING EQUATIONS AND NUMERICAL METHODS
In our axisymmetric simulation, the initially-spherical droplet with a radius of R issurrounded by the bulk liquid, which is imposed on a constant temperature gradient alongthe z axis (Fig. 1). The calculated domain is Ω(( r, z )) ( r ∈ [0 , r ], z ∈ [0 , z ]), and thecenter of the droplet is at r = 0. The governing equations for the thermocapillary migrationof droplets are: ∇ · u = 0 , (3) ∂ ( ρ u ) ∂t + ∇ · ( ρ uu ) = −∇ p + ∇ · ( µ ( ∇ u + ∇ T u )) + F σ , (4) ρC p ( ∂T∂t + u · ∇ T ) = ∇ · ( k ∇ T ) , (5)in which u = ( v r , v z ) is the velocity, ρ the density, p the pressure , and C p the specific heat. F σ is the body force term due to the interfacial tension, and its axisymmetrical formula ispresented in [20, 21]: F σ = Z B δ ( x − x f )( σκ n + ∂σ∂s τ ) ds. (6)4ere, x = ( r, z ) is the space vector, x f = ( r f , z f ) the position of the cell f on the interface B , and s the natural coordinate along the interface. n = ( n r , n z ) and τ denote the normaland tangential unit vectors of the interface, respectively. κ = κ + κ is the sum of twoprincipal curvatures of the interface, κ is the in-plane curvature and κ is given as [22]: κ = − n r r f . (7)We defined nondimensional variables as follows:¯ u = u /U, ¯ x = x /R , ¯ t = t/ ( R U ) , ¯ p = p/ ( ρ U ) , ¯ T = T / ( |∇ T ∞ | R ) , ¯ ρ = ρ/ρ , ξ = ρ /ρ , ¯ µ = µ/µ , (8)¯ k = k/k , ¯ C p = C p /C p , γ = C p /C p , ¯ F σ = F σ R / ( ρ U ) , M a = ρ C p U R /k . In total, the problem is governed by seven nondimensional parameters:
Re, M a, Ca, α, λ, γ and ξ . The nondimensional equations in the axisymmetrical model are: ∂ ¯ v r ∂ ¯ r + ¯ v r ¯ r + ∂ ¯ v z ∂ ¯ z = 0 , (9) ∂ ( ¯ ρ ¯ v r ) ∂ ¯ t + 1¯ r ∂ (¯ r ¯ ρ ¯ v r ¯ v r ) ∂ ¯ r + ∂ ( ¯ ρ ¯ v z ¯ v r ) ∂ ¯ z = − ∂ ¯ p∂ ¯ r + 1 Re " r ∂∂ ¯ r (¯ r ¯ µ ( ∂ ¯ v r ∂ ¯ r )) + ∂∂ ¯ z (¯ µ ( ∂ ¯ v r ∂ ¯ z + ∂ ¯ v z ∂ ¯ r )) − µ ¯ v r ¯ r + ¯ F σr , (10) ∂ ( ¯ ρ ¯ v z ) ∂ ¯ t + 1¯ r ∂ (¯ r ¯ ρ ¯ v r ¯ v z ) ∂ ¯ r + ∂ ( ¯ ρ ¯ v z ¯ v z ) ∂ ¯ z = − ∂ ¯ p∂ ¯ z + 1 Re " r ∂∂ ¯ r (¯ r ¯ µ ( ∂ ¯ v z ∂ ¯ r + ∂ ¯ v r ∂ ¯ z )) + 2 ∂∂ ¯ z (¯ µ ( ∂ ¯ v z ∂ ¯ z )) + ¯ F σz , (11)¯ ρ ¯ C p ∂ ¯ T∂ ¯ t + ¯ v r ∂ ¯ T∂ ¯ r + ¯ v z ∂ ¯ T∂ ¯ z ! = 1 M a " r ∂∂ ¯ r (¯ r ¯ k ∂ ¯ T∂ ¯ r ) + ∂∂ ¯ z (¯ k ∂ ¯ T∂ ¯ z ) . (12)Using Eq. 1, the nondimensional interfacial tension term ¯ F σ can be written as:¯ F σ = Z B δ (¯ x − ¯ x f )(( 1 ReCa − Re ( ¯ T − ¯ T )) κ n − Re ∂ ¯ T∂ ¯ s τ ) d ¯ s. (13)Eqs. 9-12 are valid for both phases of the drop and the bulk fluid. The discontinu-ities of the physical parameters across the interface are handled with an indicator function(see [23] for more details). The governing equations are discretized by the second-order5enter-difference method on fixed, staggered, Cartesian grids. The projection method [24]is adopted to solve the Navier-Stokes equations.On the symmetric axis, the symmetrical condition is employed:¯ v r | ¯ r =0 = 0 , ∂ ¯ v z ∂ ¯ r | ¯ r =0 = 0 , ∂ ¯ T∂ ¯ r | ¯ r =0 = 0 , (14)and all other boundaries for velocities adopt rigid boundary conditions, and those for tem-perature use constant temperature boundary conditions. Following the tradition in thisfield, the initial conditions are:¯ v r | ¯ t =0 = ¯ v z | ¯ t =0 = 0 , ¯ T | ¯ t =0 = ¯ z. (15)In the following, symbols without bars will be used to denote nondimensional values.In addition, when the momentum equation is solved at each new time step, the frontelement f is supposed to have the same velocity ( u f ) as that of the surrounding fluid. Inorder to update the front with a conserved drop volume, a velocity correction ∆ u n is imposedon u f in the normal direction ( n f ) of the element:∆ u n = − P f ∆ s f u f · n f P f ∆ s f . (16)With this correction, the drop volume is almost conserved throughout our simulations witha maximum volume loss of less than 0.1% of the initial volume. By comparison, without thevelocity correction, there will be about 2% loss at the simulating time t = 20.The computing domain in this paper is ( r, z ) ∈ [0 , × [0 , × × − for most cases, and 1 × − for the ξ = 0 .
01 case insubsection 4.2.
III. THE VALIDATION OF THE COMPUTING PROGRAM
In tradition, the scaled deviation of the drop profile from the sphere is defined as f ( θ ) = R ( θ ) /R − , where, the polar angle θ is measured from the front stagnation point and R ( θ ) denotes thedistance between the interface and the mass center of the drop (Fig. 1).6 f ( θ ) FIG. 2: Deviation of the drop profile from the sphere for
M a = 1, Ca = 0 . α = λ = ξ = 0 . γ = 0 .
25, and Re as indicated. Dash line: the theoretical prediction in [17]; solid line: the presentwork. θ f ( θ ) FIG. 3: Deviations of drop profiles from spheres in [17], [18], and this work, where Re = 25, M a = 1, Ca = 0 . α = λ = ξ = 0 .
5, and γ = 0 . Before this work, there was a theoretical result of small inertial deformations of dropswith small Ca and Re numbers [17]: f ( θ ) = 38 A ReCa ( ξ − θ − . (17)The earlier numerical simulation [18], however, used a different formula for the interfacialtension: F σ = Z B δ ( x − x f )( 1 ReCa κ n − Re ∂T∂s τ ) ds. (18)Compared with Eq. 13, the effect of capillary force in Eq. 18 was assumed constant on thenormal direction of the interface, or, there is no decrease in the interfacial tension when the7roplet moves to the hotter region [18]. There are two distinguished differences betweenthese two assumptions:1. As it will be shown later in this paper, Eq. 13 leads to a bigger drop deformation thanEq. 18.2. With Eq. 13, the simulation must be stopped before the drop moves a distance of(1 /Ca −
1) in the z direction. This is because the interfacial tension of the drop willbe smaller than zero beyond that location, and not only the simulation will collapse,but also it is not physically possible.With Eq. 18, on the other hand, the simulation can be extended to the infinity.In this section, Eq. 18 is adopted in our codes to have some comparisons with theprevious researches. It seems that our simulations have very good agreements with theanalytical work (Fig. 2), and that a better similarity is achieved when the Reynolds numberis smaller. Compared with the previous numerical work [18], our numerical simulations arecloser to the theoretical analysis (Fig. 3).In the rest of this paper, the interfacial tension is calculated by Eq. 13, and simulationsare stopped before the capillary force on the drop becomes negative. IV. RESULTS AND DISCUSSIONSA. Transient thermocapillary migrations of droplets with large capillary numbers
In this subsection, we simulate the thermocapillary process with the same set of param-eters in [18]: Re = M a = 50, Ca = 0 . α = λ = ξ = 0 .
5, and γ = 0 .
25. The time evolutionof migration velocity is shown in Fig. 4 (the solid line). When the capillary number is suffi-ciently small ( Ca = 0 .
05, the dashed line in Fig. 4), the drop will reach a steady migratingstate after the initial accelerating process. However, the velocity with Ca = 0 . t >
60, and can not reach a steady state. This result is differentfrom that of [18], in which the steady velocity state was assumed in advance.The deviation of the drop profile from the sphere for Ca = 0 . V FIG. 4: Time evolution of drop velocity for Re = M a = 50, α = λ = ξ = 0 .
5, and γ = 0 . Ca = 0 . Ca = 0 .
05 (dashed line). θ f ( θ ) FIG. 5: Deviation of the drop profile from thesphere at indicated times for Re = M a = 50, Ca = 0 . α = λ = ξ = 0 .
5, and γ = 0 . the drop becomes flattened, and the drop looks like a cap (Fig. 6(b)). The deformationat t = 60 in the current study is about twice bigger than that in the simulation of [18],and even bigger difference at the late stage. There are two reasons for the rapid decline invelocity:1. Because the drop becomes flatter, the resistance on the drop is increased.2. Because the distance between the front and rear stagnation points at the late stage isshorter than that in the beginning, their temperature difference also decreases (Fig.7). As a result, the thermocapillary driving force on the drop becomes smaller.We also studied the influence of different Re numbers, M a numbers, and Ca numbers.Some typical results are shown in Figs. 8. With other parameters fixed, the migratingvelocity tails off earlier for the larger Reynolds number.With Re = M a = 50, Ca = 0 .
1, and ξ = 0 .
5, we studied the influence of α , λ and γ .These three parameters play trivial roles in the transient motions of deformable drops. Withsmall α , λ or γ , drops migrate faster (see [5, 25]) and need shorter time to reach the hotregions with near-zero interfacial tensions. As a result, the migration velocities tail off atearlier times.To sum up, the influence of all other parameters except ξ depends on when the dropmigrates to the warm region with very low interfacial tension, in other words, the speed of9 . . r z . . . r z FIG. 6: Isotherms around the drop for Re = M a = 50, Ca = 0 . α = λ = ξ = 0 .
5, and γ = 0 . t = 40; b) t = 70. θ T ( θ )- T ( ) FIG. 7: The temperature distribution on the drop interface for Re = M a = 50, Ca = 0 . α = λ = ξ = 0 .
5, and γ = 0 .
25 when t = 40 and t = 70. V (a) Re = 1 t V (b) Re = 50 FIG. 8: Time evolutions of the drop velocity for different Ca numbers. M a = 100, α = λ = ξ = 0 . γ = 0 .
25. (a) Re = 1; (b) Re = 50. the drop. The influence of various parameters on drop speeds has been discussed in detailin our earlier work [5, 26]. B. Influence of density ratio ξ In the last subsection, densities of drops are always smaller than those of bulk fluids, andwe will adopt denser drops ( ξ = 2) in the following. Several Ca numbers are studied, allother parameters are fixed to make the discussion simpler ( Re = M a = 50, α = λ = 0 . γ = 0 . Ca = 0 . z direction (Fig.10(b)), and the temperature difference between the front and rear stagnation pointsof the drop becomes larger (Fig. 11). As a result, the thermocapillary driving forceon the drop becomes larger.2. Because the drop becomes slender, the resistance on the drop decreases.To further discuss drop deformations, it is essential to introduce the local capillary number Ca l to indicate the local magnitude of interfacial tension, which is defined on the moving11 V FIG. 9: Time evolutions of drop velocities for different capillary numbers. Re = 50, M a = 50, α = λ = 0 . γ = 0 .
25, and ξ = 2. . . . r z (a) t= 40 . . r z (b) t = 80 FIG. 10: Isotherms around the drop at t = 40 and t = 80. Re = 50, M a = 50, Ca = 0 . α = λ = 0 . γ = 0 .
25, and ξ = 2. T ( θ )- T ( ) FIG. 11: Temperature distribution on the drop interface at t = 40 and t = 80. Re = 50, M a = 50, Ca = 0 . α = λ = 0 . γ = 0 .
25, and ξ = 2. t V ξ = 0.01 ξ = 0.5 ξ = 1 ξ = 2 ξ = 4 (a) θ f ( θ ) ξ = 0.01 ξ = 0.5 ξ = 1 ξ = 2 ξ = 4 (b) FIG. 12: (a) Time evolutions of drop velocities with different ξ ’s, the circles on the velocity curvesindicate the locations at Ca l = 0 .
5; (b) deviations of drop profiles from spheres when Ca l = 0 . Re = M a = 50, Ca = 0 . α = λ = 0 .
5, and γ = 0 . drop center ( z c ): Ca l = Ca − Ca ( T c − T ) = Ca − Ca ( z c − z ) , (19)where, T c represents the initially temperature at z c . Time evolutions of drop velocities fordifferent ξ ’s are plotted in Fig. 12(a). To investigate the ξ influence on deformation, wefixed Ca l at the value of 0 .
5, when the drop starts to have an obviously different velocityfrom that of the non-deformable drop. In the current case ( Ca = 0 . Ca l = 0 . ξ = 4 ξ = 0.01 ξ = 0.5 ξ = 1 ξ = 2t = 70Ca l = 0.5 FIG. 13: Time evolutions of drop profiles for different ξ ’s, where Re = M a = 50, Ca = 0 . α = λ = 0 .
5, and γ = 0 .
25. From the bottom to the top, the profiles in each column arecorresponding to the time at t = 0 , , , , , that the drop migrates a distance of 8 (represented by circles on the curves of Fig. 12(a)). Itis found that the drop profile deforms to oblate/slender when ξ is smaller/larger than unit.When the absolute value of ξ − Re & M a numbers here. Fig. 13 presents shape evolutions of drops with different ξ ’s.It should be noticed that only slight deformations on large bubbles ( ξ ∼ − ) have beenobserved in previous space experiments, and no measurable deformations for those heavydrops ( ξ = 1 .
98) [3]. The reason for the difference between experiments and simulationsmight be:1. Bubbles/Drops in experiments do not migrate to areas with high temperatures, sotheir interface tensions are still large enough to hold perfect spherical shapes;2. The assumption for the σ ( T )’s linear dependence on temperature is not correct when14 + 1/2 ρ V p + 1/2 ρ V + 2 σ H FIG. 14: Sketch of normal stress balance on the stagnation point of the drop. σ ( T ) → ξ ’s and drop-deforming orientations, a simplifiedphysical explanation is introduced in the following. Only the front stagnation point isstudied.As shown in Fig. 14, the total pressures are p + ρ V inside the drop and p + ρ V outside the drop; here V is the local fluid velocity, and p and p denote the static pressuresat the two sides of the stagnation point, respectively. The extra pressure caused by theinterfacial tension is 2 σH (H is the average curvature). When the interfacial tension is largeenough, the normal stress balance on the interface is dominated by the interfacial tensionand the droplet keeps spherical. On the other hand, when the interfacial tension becomessmall at the hot region (or, Ca l > . if ξ > , the dynamic pressure inside the drop ( ρ V ) is larger than that outside ( ρ V ).Larger extra pressure (2 σH ) is needed to balance the normal stress on the interface,so the drop is elongated in the axis direction to have a larger H . if ξ < , the dynamic pressure inside is lower than that outside. Smaller extra pressure isneeded to balance the normal stress, so the drop is compacted in the axis direction tohave a smaller H . V. CONCLUSIONS
In this work, we found that the influence of deformations on drops is much more com-plicated than that on bubbles. When the drop density is smaller than the bulk fluid, the15igration velocity will decrease with the increase of the deformation. When the drop densityis larger than the bulk fluid, the migration velocity will increase with the increase of thedeformation. With the assumption adopted in this paper, we believe that in the hot regionthe drop will always start to deform, and there is no way to have any constant migratingvelocity. Hence, keeping Ca and Ca l small throughout simulations is the only possible wayfor the drop to keep spherical and reach the steady migrating state.All the discussions above are based on the assumption of the linear temperature depen-dence of the interfacial tension. More complicated temperature dependence will be usedin the future work, and variations of physical properties with the temperature will also beconsidered. Acknowledgments
This project is supported by the Knowledge Innovation Program of the Chinese Academyof Sciences, Grant No. KJCX2-YW-L08. [1] Ratke L, Korekt G. Solidification of Al-Pb-base alloys in low gravity. Z Metallkde, 2000, 91:919-927[2] Ostrach S. Low gravity fluid flows. Annu Rev Fluid Mech, 1982, 14:13-345[3] Hadland P H, Balasubramaniam R, Wozniak G, et al. Thermocapillary migration of bubblesand drops at moderate to large Marangoni number and moderate Reynolds number in reducedgravity. Exp Fluids, 1999, 26:240-248[4] Hoyer W, Kaban I. Experimental and calculated liquid-liquid interface tension in dimixingmetal alloys. Rare Met, 2006, 25(5):452-456[5] Yin Z H, Gao P, Hu W R, et al. Thermocapillary migration of nondeformable drops. PhysFluids, 2008, 20: 082101[6] Carlberg T, Fredriksson H. The influence of microgravity on the solidification of Zn-Bi immis-cible alloys. Metal Trans A, 1980, 11(10):1665-1676[7] Li H L, Zhao J Z. Convective effect on the microstructure evolution during a liquid-liquiddecomposition. Appl Phys Lett, 2008, 92(24):241902
8] Gurgel A, Moura M C P A, Dantas T N C, et al. A review on chemical flooding methodsapplied in enhanced oil recovery. Braz J Pet Gas, 2008, 2(2):83-95[9] Xia H F, Liu R Q, Ju Y, et al. Mechanism of betaine surfactant solution on residual oil afterwater flooding in ultralow interfacial tension. J Daqing Pet Inst, 2006, 30(6):24-27(in Chinese)[10] Lyford P A, Pratt H R C, Shallcross D C, et al. The marangoni effect and enhanced oilrecovery Part 1. Porous media studies. Can J Chem Eng, 1998, 76(2):167-174[11] Lyford P A, Pratt H R C, Shallcross D C, et al. The marangoni effect and enhanced oil recoveryPart 2. Interfacial tension and drop instability. Can J Chem Eng, 1998, 76(2): 175-182[12] Morrow N R, Mason G. Recovery of oil by spontaneous imbibition. Curr Opin Coll Int Sci,2001, 6(4):321-337[13] Young N O, Goldstein J S, Block M J. The motion of bubbles in a vertical temperaturegradient. J Fluid Mech, 1959, 11: 350-356[14] Subramanian R S, and Balasubramaniam R. The Motion of Bubbles and Drops in ReducedGravity. Cambridge:Cambridge University Press, 2001[15] Subramanian R S, Balasubramaniam R, Wozniak G. Fluid Mechanics of Bubbles and Drops.London: Taylor & Francis, 2002:149-177[16] Balasubramaniam R, Chai A. Thermocapillary migration of droplets: an exact solution forsmall marangoni numbers. J Colloid Interface Sci, 1987, 119:531-538[17] Haj-Hariri H, Nadim A, Borhan A. Effect of inertia on the thermocapillary velocity of a drop.J Colloid Interface Sci, 1990, 140(1):277-286[18] Haj-Hariri H, Shi Q, Borhan A. Thermocapillary motion of deformable drops at finite Reynoldsand Marangoni numbers. Phys Fluids, 1997, 9(4):845-855[19] Welch S W J. Transient thermocapillary migration of deformable bubbles. J Colloid InterfaceSci, 1998, 208: 500-508[20] Unverdi S O, Tryggvason G. A front-tracking method for viscous, incompressible, multi-fluidflows. J Comput Phys, 1992, 100:25-37[21] Tryggvason G, Bunner B, Esmaeeli A, et al. A front-tracking method for the computations ofmultiphase flow. J Comput Phys, 2001, 169:708-759[22] Kreyszig E. Differential Geometry. New York: Dover Publications, 1991[23] Hua J S, Lou J. Numerical simulation of bubble rising in viscous liquid. J Comput Phys, 2007,222:769-795
24] Chorin A J. A numerical method for solving incompressible viscous flow problems. J ComputPhys, 1967, 2:12-26[25] Nas S. Computational Investigation of Thermocapillary Migration of Bubbles and Drops inZero Gravity. Doctor Dissertation. Michigan: The University of Michigan, 1995[26] Gao P. Numerical Investigation of the Drop Thermocapillary Migration. Doctor Disserta-tion(in Chinese). Beijing: Graduate University of Chinese Academy of Science, 2007[27] Moldover M R. Interfacial tension of fluids near critical points and two-scale-factor universality.Phys Rev A, 1985, 31(2):1022-103324] Chorin A J. A numerical method for solving incompressible viscous flow problems. J ComputPhys, 1967, 2:12-26[25] Nas S. Computational Investigation of Thermocapillary Migration of Bubbles and Drops inZero Gravity. Doctor Dissertation. Michigan: The University of Michigan, 1995[26] Gao P. Numerical Investigation of the Drop Thermocapillary Migration. Doctor Disserta-tion(in Chinese). Beijing: Graduate University of Chinese Academy of Science, 2007[27] Moldover M R. Interfacial tension of fluids near critical points and two-scale-factor universality.Phys Rev A, 1985, 31(2):1022-1033