- Email: [email protected]
The stokes motion of a gas bubble due to interfacial tension gradients at low to moderate Marangoni numbers
The Stokes Motion of a Gas Bubble Due to Interfacial Tension Gradients at Low to Moderate Marangoni Numbers N . S H A N K A R 1 AND R. S H A N K A R
SUBRAMANIAN
Department of Chemical Engineering, Clarkson University, Potsdam, New York 13676 Received May 29, 1987; accepted July 21, 1987 The steady migration velocity of a gas bubble placed in a large liquid body possessing a uniform temperature gradient in the undisturbed state is calculated under conditions of negligible inertia. Proper account is taken of convective energy transport effects, and results are obtained by solving the energy equation using a numerical scheme. The range of utility of an expansion for the bubble velocity valid in the limit of small Marangoni numbers is extended by the use of series improvement techniques. The improved series is good for Ma ~< 20 whereas the original series is useful for only Ma ~ 1. Also, an excellent empirical fit of the numerical results for the bubble velocity is obtained in the limit of large values of the Marangoni number. This result is good for Ma -~-20. Streamline plots are used to demonstrate the existence of a recirculafion region which moves closer to the bubble as the Marangoni number is increased. Isotherms are included to display the presence of an extended thermal wake behind the bubble whose extent increases as the Marangoni number is increased. © 1988AcademicPress,Inc.
Thus, under such conditions, the interface can become the principal driver of fluid motion. A gas bubble will migrate in a fluid if subAn example wherein interfacial tension grajected to a gradient of interfacial tension. The dients can be used to remove gas bubbles from motion is a consequence of the discontinuity a glass melt processed in space is illustrated in the tangential stress across the bubble surby Subramanian (1). face. The traction exerted on the fluid adjacent The fact that gas bubbles will move due to to the interface will cause it to move in the interfacial tension gradients was first demondirection of the pole with the greatest interstrated experimentally by Young et al. (2). facial tension, and by reaction, the bubble will Subsequent experimental work as well as themove in the opposite direction. ory has been reviewed earlier (3, 4), and will Interfacial tension gradients can be mainnot be discussed here. Our attention, in the rained at a bubble surface by arranging to have present article, will be devoted to the theorettemperature or concentration gradients in the ical problem of the steady motion of a gas continuous phase fluid. As a mechanism for bubble in a large liquid body possessing a unidriving the motion of bubbles, when compared form temperature gradient in the undisturbed to a body force field such as gravity, interfacial state. tension gradients become significant at small Young et al. (2) considered the above probbubbles sizes. Furthermore, in "space prolem under conditions wherein convective cessing," a term which describes working with transport effects were completely negligible. In materials and processes in the near free fall their treatment, the general case of a fluid environment aboard orbiting spacecraft, gravdroplet was analyzed. For a gas bubble, one itational effects are relatively unimportant. can make the simplification that the viscosity and the thermal conductivity of the gas are 1 Present address: Martin Marietta Laboratories, 1450 South Rolling Road, Baltimore, M D 21227. negligible compared to the same properties in 1. INTRODUCTION
512 0021-9797/88 $3.00 Copyright© 1988by AcademicPress,Inc. All rightsof reproductionin any formreserved.
Journalof ColloidandInterfaceScience.Vol. 123,No. 2, June 1988
MOTION
the surrounding liquid. In this special case, in the absence of gravitational effects, the result given by Young et al. reduces to u =
a l d a / d T [ IVTo~ [ 2~
[1]
Here, u is the steady migration velocity of the bubble and VT~ is the undisturbed constant temperature gradient in the fluid, d a / d T is the gradient of interfacial tension with temperature, assumed constant, a is the radius of the bubble, and # is the viscosity of the continuous phase, also assumed constant. One direction in which the assumptions of Young et al. can be relaxed is to include the contribution from convective transport. In applications to highly viscous fluids, such as glass melts, the Reynolds numbers will be quite negligible, and inertial effects can be safely ignored. However, the Prandtl numbers in such systems can be very large (10 4 and larger). Thus, the Peclet number for heat transfer can be substantial, and convective energy transport cannot always be neglected. In this problem, the Peclet number when based on an appropriate velocity scale, is commonly designated the Marangoni number. The above problem was considered by Subramanian (3) for small values of the Marangoni number, Ma. He obtained a solution of the governing field equations by the method of matched asymptotic expansions, and gave a result for the bubble speed good to 0(Ma2), U=
1
2
301 --Ma 14,400
2+ . - . .
[2]
In Eq. [2], U is the scaled bubble velocity (U = U/Vo). The reference velocity is Vo = (al da/ dTIl~TZoo])/~. Thus, at zero Marangoni number the result reduces correctly to that predicted by Young et al. Subsequently, Merritt (5) found the next higher term in this expansion, U-
1
2
301 - Ma 2 14,400 13 + 1 - ~ Ma3 + " " ""
513
OF A GAS BUBBLE
[3]
It is interesting to note that, to this order, there appear to be no logarithmic terms in the asymptotic sequence. This is in contrast to the behavior of a similar heat transfer problem associated with a rigid sphere in creeping motion through a fluid which was studied by Acrivos and Taylor (6). In principle, one can carry the calculation to higher orders. However, this process is a formidable exercise in algebra. Furthermore, it is useful to be able to check the utility of the above expressions by comparison with an independently obtained solution of the governing equations. This would permit us to determine the range of Marangoni numbers over which Eqs. [2] or [3] are useful. For the reasons cited above, we have solved the governing equations here numerically using the method of finite differences, and have calculated the bubble velocity up to Ma = 200. Furthermore, we have used series improvement techniques to discover useful information concealed in the coefficients in Eq. [3] and to develop an analytical result which appears to be good up to Ma = 25. In contrast, Eq. [3] diverges from the numerical results beyond Ma = 1. In Section 2, the theoretical model equations are stated, and the method of solution is discussed briefly. In Section 3, results from the calculations are displayed and discussed. We conclude with a few remarks in Section 4.
2. A N A L Y S I S
Consider a gas bubble present in an unbounded liquid body which, in the undisturbed state, possesses a uniform temperature gradient, VToo. Gravitational effects are neglected, and only the motion of the bubble due to the interfacial tension variation resulting from the temperature gradient on the interface is considered. It is assumed that, over the spatial domain of interest, the physical properties of the liquid are constant with the exception ofinterfacial tension. This quantity is assumed to decrease linearly with increasing Journal of Colloid and Interface Science, Vol. 123, No. 2, June 1988
5 14
S H A N K A R A N D S U B R A MA N IA N
temperature. Furthermore, it is assumed that the bubble size does not change appreciably during the migration process. Under the above assumptions, it is possible for the interfacial tension gradient on the bubble surface to achieve a steady representation, resulting in steady migration of the bubble in the direction ofVT~. This is the problem we shall analyze below. As stated earlier, inertial effects are ignored, and the motion of the fluid is assumed to be incompressible and Newtonian. The viscosity and the thermal conductivity of the gas phase are assumed negligible compared to the corresponding properties in the liquid. Thus, only the transport problem in the liquid needs to be treated. In general, the bubble will be deformed from the spherical shape even when inertial effects are ignored. This is because the interfacial tension is not a constant over the bubble surface. However, provided the variation in interfacial tension over the bubble surface is small when compared to the value of the interfacial tension, the capillary number will be small, and the bubble boundary may be assumed spherical with little error. It is convenient to use a reference frame traveling with the bubble. As recognized by Bratukhin (7), an observer riding with the bubble will perceive an unsteady temperature field at all times. However, the important entry which determines the achievement of a steady migration velocity, subject to the assumptions made earlier, is the temperature gradient field. The gradient field, after an initial transient, will achieve an asymptotic steady state. To obtain this, as done by Bratukhin (7), we pose a problem for the difference between the temperature at any point and the temperature at the bubble's equatorial plane at infinity. This difference field will reach a nontrivial steady state after an initial transient, and its gradient is identical to that of the actual temperature field. The following governing equations may be written for the scaled velocity and temperature fields, Journal of ColloM and Interface Science, Vol. 123,No. 2, June 1988
VZv = Vp
[4]
Ma[U + v. X7T] = V2T.
[51
In the above, distances are scaled by the bubble radius, a, velocities by v0, defined earlier, and the dynamic pressure by IVZo~l Ida~ dTI, The scaled temperature is the difference between the actual temperature and that at infinity on the equatorial plane of the bubble, divided by a lVTo~I. The Marangoni number :is defined by Ma = (aVo)/a where o~is the thermal diffusivity of the liquid. The boundary conditions on the fields may be written as follows. Far away from the bubble, the velocity and temperature fields attain their free stream values, v --~ -Uiz ]
[6] as
r--~ ~ .
T --~ z
[7]
Here, iz is-a unit vector in the direction of the temperature gradient VToo in the undisturbed fluid, and z is-distance measured from the equatorial plane in this direction, r is the position vector with the origin at the bubble center. At the bubble surface, Sb, the vanishing of the normal velocity, the balance of tangential stress, and the vanishing of the normal energy flux yield the following boundary conditions. v.nb = 0
]
rib" ~r- tb = --Vs T. t b l XTT.n~ = 0
[8] on
Sb.
[9] [10]
In the above, ,c is the scaled hydrodynamic stress tensor in the liquid. Vs is the surface gradient operator (8), nb is a unit normal vector to the surface, and tb is an arbitrarily oriented unit vector lying in a tangent plane on the bubble surface. Finally, setting the hydrodynamic force on the bubble to zero completes the problem statement. A formal solution for the velocity field may be written immediately from a known general solution of the Stokes equation in spherical polar coordinates by suitable specialization.
M O T I O N OF A GAS BUBBLE
The solution, given by Subramanian (3), contains an infinite set of constants which are coefficients in the expansion of the temperature field on the bubble surface in a series of Legendre polynomials. It is the appearance of these constants in the velocity field which makes the energy equation nonlinear. We solved the above problem for the temperature field using the method of finite differences. In doing so, a new temperature field, T' = T - z, was defined to facilitate the application of the boundary condition at infinity. Also, in recognition of the existence of the sharpest gradients near the bubble surface, the radial coordinate r, measured from the bubble center, was transformed via ~ = log(r). This transformation was employed by Dennis et al. (9), who studied steady heat transfer in forced convection past a solid sphere, and found to be satisfactory. In writing the finite difference equations, central differences were written except for the convective transport terms which were handled using an upwind differencing scheme, as discussed by Roache (10) and Gosman et al. (11). As observed by Dennis et al. (9), the accuracy of the solution can be improved by further replacing the infinity condition by an "improved" one which can be used at a sufficiently large value of the radial coordinate (see also Ref. (12)). By following their procedure which involves the recognition that at large distances from the bubble the velocity field may be closely approximated by the free stream (the Oseen approximation), followed by the solution of the appropriate Oseen equation, the ratio of the temperature fields at adjacent grid points can be obtained. After sufficient testing of this condition against the earlier one, the improved condition was employed in all the calculations at ~ = 3. This corresponds to a location approximately 20 bubble radii from the bubble center. It was verified that using the condition at ~ = 2 gave bubble velocities which were different by no more than 0.5%. In some runs at low values of the Marangoni number, the infinity condition was applied at further radial locations such as 100 bubble ra-
515
dii to confirm that there was no change in the results. An 81 X 81 grid in the two spatial coordinates was employed in all the calculations, and the results checked for accuracy by using a 51 X 51 grid at selected values of the Marangoni number. The latter produced results within 1% of the former over the entire range of values of the Marangoni number investigated here. Typically, on an IBM 4341 computer, it was found that it took approximately 30 rain of CPU time to perform the entire set of computations for a single value of the Marangoni number. Values of the Marangoni number ranging from 0 to 200 were covered. Further details concerning the numerical solution procedure may be found in Shankar (13). 3. RESULTS AND DISCUSSION
Streamlines Typical streamlines from the numerical calculations are displayed in a meridian section in Figs. 1a-1 c for an increasing progression of values of the Marangoni number. The streamlines for Ma = 200 were plotted, but were found to be virtually identical to those for Ma = 100, and hence are not included. Figure 1a corresponds to Ma = 0 and was obtained from the analytical solution given by Young et al. All the plots are given as seen by an observer in a laboratory reference frame. The bubble is moving upward in all the figures. The velocity field is sensitive to the value of the Marangoni number through its dependence on the surface temperature field. The single most striking feature in Figs. 1a - l c is the appearance of an unattached recirculation region in the wake of the bubble. There is perfect fore-aft symmetry in the plot for Ma = 0 since the temperature gradient field on the bubble surface in this case is perfectly symmetric about the equator of the bubble. Such symmetry is lost in the presence of convective energy transport. A streamline on which the scaled Stokes stream function ff = 0 appears in the wake region, and separates the recirJournal of ColloM and Interface Science, Vol. t23, No. 2, June 1988
516
SHANKAR
AND SUBRAMANIAN
I
2-o.ot 1171-o.I I 3 -o.oz II8 o.ool I -oo4 191 o ooz5 I
I -0.005 5-0.08 2 -0.02 6 "0.1 3 -0.04 7 -0.125
4-0.06
-0.06
]18-0.15
0 0.004
I -0005 6-o.o8 [ z -o.oz Ilrl-°.l / 3-0.03
8 0.004 I
4-o.o4 -o.o6
91 0.006 / IOl o.ooe /
C
FIG. 1. (a) Streamlines, M a = 0. (b) Streamlines, M a = 5. (c) Streamlines, M a = 100.
culating region from the rest of the flow field. This may be regarded as a separatrix. The recirculating "eddy" in Figs. lb and lc must close at infinity. The intensity of flow in this region appears to increase with increasing values of the Marangoni number. Also, the recirculation pattern moves closer to the rear Journal of Colloid and Interface Science, Vol. 123, No, 2, June 1988
stagnation point as the Marangoni number is increased. It is interesting to determine whether the above features can be predicted from the small Ma expansion. Subramanian's result for the stream function in a laboratory reference frame, good to 0(MaZ), may be written as
MOTION
OF A GAS BUBBLE
517
The above equation has one real root and a pair of complex conjugate roots for Ma ~<4.76. At Ma = 0.25, the real root is 34.2. 3 This compares favorably with the finite difference result, r0 = 34.6. At Ma = 1, the above equation yields r0 = 8.2 whereas the numerical where solution gives r0 = 9.3. Thus, for small values I2 = - 1 ÷ ~ Ma 2 [12a] of Ma, Eq. (14) can be used to obtain a good estimate of r0. The physical origin of the recirculating eddy 7 13 = 1 ~ Ma [ 12b] can be identified. When convective energy transport is neglected, the resulting tempera29 ture field on the bubble surface is a pure Le14 - Ma 2 [ 12c] 11200 gendre mode, namely Pl(s). The surface gradient of this field appears as the tangential Here, r is the radial coordinate in a spherical stress (in scaled form) and drives the flow in polar coordinate system centered in the bubthe liquid. The basic flow due to this temperble. s = cos 0 where 0 is the polar angle measured from the positive z-direction defined ature gradient is illustrated in Fig. 1a. Crucially important is the fact that the magnitude of the earlier. It may be seen that Eq. [11] correctly pre- velocity in the fluid decays away from the dicts the existence of a dividing streamline or bubble as 1/r 3. When convective transport of energy is inseparatrix for nonzero values of the Marangoni number. Setting ff = 0 leads to the axis of cluded, higher Legendre modes of the surface symmetry as one solution (s = _+1), and the temperature field are excited. For n >1 2, each equation for the separatrix as the other solu- pure mode Pn(s) leads to a flow which decays tion: as r-n as r ~ ~ . Thus, it is clear that flow at large distances from the bubble must be dominated by the P2 mode of the surface temper4r ~I3 -ls ature field. This flow is toward the bubble in the equatorial region, and away from it in the +-~14 (5s 2 - 1 ) = 0 . [13] polar regions. Thus, in the region of the north pole (forward region) of the bubble, this flow When streamlines are plotted from Eq. [ 11 ], is in the same general direction as the basic for representative values of the Marangoni flow due to the P1 mode, and the two simply number, the qualitative features are similar to reinforce each other. In the rear, the two opthose presented in Figs. lb and lc. Further- pose each other. Thus, flow must be away from more, Eq. [ 13] may be used to predict the value the bubble at large distances into the wake and of ro where the separatrix intersects the axis of toward the bubble in its vicinity. This, coupled symmetry. Substituting the coordinates of this with the inward flow at the equator induced point (r0, - 1 ) in Eq. [13], the following cubic by the P2 mode, leads to the observed separated equation is obtained: eddy which moves closer to the bubble as the P2 mode gains strength. As a matter of interest, / 289 r03 + ~7--~ Ma - ~Ma)r~ the P3 mode of the temperature field induces flows which reinforce those due to the P1 mode 261 in the polar regions, and which decay also - r o - 1--9~ Ma = 0. [141 as 1/r 3.
s2,{4 +
3
1
1)s
),5s21,}
Lll
Journal of Colloid and Interface Science, Vol. 123, No. 2, June 1988
518
SHANKAR AND SUBRAMANIAN
Isotherms Isotherms in a meridian plane are plotted for a series of values of the Marangoni number in Figs. 2a-2d. Figure 2d is drawn for the same value of Ma (=200) as Fig. 2c, but over a greater spatial domain to show the extent of the thermal wake region. The plots clearly establish the appearance of a long thermal wake whose extent increases with increasing values of the Marangoni number. The fore-aft symmetry present in the isotherms for the case of pure conduction is destroyed in the presence of convective transport. It may be observed from the plots that the overall temperature variation over the bubble surface is reduced by the inclusion of convec-
tive energy transport. Thus, one might expect that the magnitude of the scaled migration velocity would decrease with increase in Ma. This is indeed the case as will be seen shortly. The condition of no energy flux into the bubble requires that the isotherms meet the bubble surface in a normal direction. At first, some of the isotherms in Fig. 2c appear to violate this requirement. However, this is just an artifact due to the limitation on spatial resolution in this drawing. At large values of the Marangoni number, convective energy transport is dominant in determining the temperature field almost everywhere in the liquid. However, there must exist a thermal boundary layer in the vicinity of the bubble surface where conduction and convection must balance each
Z.Q
2._66
1.0
-I.0
I
T'-2.6 T'-3.0 z~
0
T=-7 FIG. 2. (a) Isotherms, M a = 0, A T = 0.4. (b) Isotherms, M a = 5, A T = 0.4. (c) Isotherms, M a = 200, A T = 0.4. (d) Isotherms, M a = 200, A T = 1. A larger region is shown.
Journal of Colloid and Interface Science, Vo|. 123, No. 2, June 1988
519
MOTION OF A GAS BUBBLE
other. In this boundary layer, the temperature gradient field is adjusted rapidly so that the no flux condition at the bubble surface is met. In making the finite difference calculations, we were careful to check that the spatial resolution is adequate to account for the rapid variation in the gradient in this boundary layer.
Bubble Velocity The principal quest of this work is the calculation of the magnitude of the steady migration velocity of the bubble as a function of the Marangoni number. A comparison of the numerical results for this quantity against the predictions of the expansions for small Ma given in Eqs. [2] and [3] is displayed in Fig. 3a, and over a larger range of Marangoni number values in Fig. 3b. It is seen from Fig. 3a that the small Ma expansions have only limited utility, being good only up to Ma ~< 1.
a
0.7
I
I
f
06 "
~
0.5
-
......
O(M@)
-'-" -"-
PADE" E I / 2 3 PADE' E2./I~
°.°ll
: I
/
F I N I T E DIFFERENCE _ _ _ O(Mo 2) =) "
The figures also contain other approximations which are discussed below. Up to 0(Ma3), the series for the scaled bubble velocity given in Eq. [3] is a power series in Ma. However, while inclusion of the 0(Ma 3) term improves the performance of the series over that when truncation is performed at the previous order, divergence appears to occur quickly as the value of the Marangoni number is increased. This led us to explore the possibility of the presence of a singularity in the complex Ma plane which might be restricting the convergence of a possible power series for the bubble velocity (14-16). We recognized at the outset that such series improvement methods might not yield much information here for two reasons. First, only the first four coefficients are available, and one of them is zero. Thus, a Domb-Sykes plot does not prove useful. Second, in other analogous problems, the asymptotic sequence contains terms involving the logarithm of the pertur-
."
....
J D,I
...'""
I
PADE" E 2 / l ' l
0.4
,,,!
lid ~
I
"-~. ~ . ~ . ~\
J
"~",~. ,
0.3
m
m o w .J
I\ 0.2 t
\
; I
\
0.2
I
W -J
0"41
'~
I
,, FINITE DIFFERENCE --O(Mo 2) ....... O(Ma ~) . . . . PADE" C I/2"1
/"
o
t-~
i
I
I
I
2
MARANGONI
J
3
NUMBER,
1 ~ [ 4
Mo
%1
I
i
I
I
I
I
20
40
60
80
MARANGONI
NUMBER,
I00
Ma
FIG. 3a. Comparison of the scaled bubble velocity, U, as a function of the Marangoni number, Ma, calculated from the finite difference solution, the [1/2] Pad~ approximant, the [2/1] Pad~ approximant, the 0(Ma2) result, and the 0(Ma 3) result. (b) Comparison of the scaled bubble velocity, U, as a function of the Marangoni number, Ma, calculated from the finite difference solution, the [1/2] Pad6 approximant, the [2/1 ] Pad6 approximant, the 0(Ma2) result, and the 0(Ma3) result. A larger range of Ma values is shown. Journal of Colloid and Interface Science, Vol. 123, No. 2, June 1988
520
SHANKAR AND SUBRAMANIAN
bation parameter. If that were to be the case here, series improvement via the usual techniques may not be possible. Thus, we were pleasantly surprised by the results reported below. Pad6 approximants are rational fraction approximants to a power series. The [M/N] Pad6 approximant is defined via f(~) =
PM(E) QN(e)
[15] "
denominator of the [2/1] approximant (Eq. [17]) is located at Ma = -2.32. Those of the [1/2] approximant are at Ma = -3.5 and -6.83. When a singularity is located on the negative real axis, at Ma = -~0, an Euler transformation (Ma = Ma/(Ma + ~o))may be used to map the singularity to the point at infinity. Thus, we attempted this with each of the above guessed locations of a singularity. We found that the results from the [1/2] approximant's singularities performed very poorly. However, using the guess from the [2/1] approximant, considerable success was achieved. The modified expansion for U is
Here f ( 0 is the known power series in e of degree M + N. PM(0 and QN(Oare polynomials in e of degree M and N, respectively. Given a power series to a certain order, various Pad6 approximants can be constructed in a U = 0.5 - 0.1125 Ma + 2.32] straightforward fashion. The properties and use of Pad6 approximants are detailed by - 0.1123(MaMa ,/3 Baker (17). According to Van Dyke (14), "The + 2.32] " [18] approximants with M and N nearly equal posThe bubble velocity calculated from Eq. sess remarkable, though somewhat mysterious, [ 18] is compared with the numerical solution properties of analytic continuation." Hence, these approximants may be used to find the in Fig. 4. The improvement, when contrasted approximate location of singularities in a per- with the performance of the earlier approxiturbation solution by calculating the zeroes of mations, is dramatic. For Ma ~<25, the above the denominator, Qi(e). All the possible Pad6 approximation is indistinguishable from the approximants were constructed from the numerical result. Beyond that point, the nupower series in Eq. [3] and compared with the merical result continues to decrease with furnumerical solution for the bubble velocity. ther increase in the Marangoni number while The best performers were the [1/2] and the Eq. [18] gives a prediction which levels off at [2/1 ] approximants. These are given below and approximately Ma ~ 50. For the benefit of the user, it is good to fit shown in Figs. 3a and 3b: the numerical results for U versus Ma for large Udl/2] values of Ma to a simple equation. In addition, such a fit can lend insight into the analytical 1 + (130/301)Ma [161 structure of an asymptotic expansion for the 2 + (260/301)Ma + (301/3600)Ma 2 fields in the limit of large values of Ma. Thus, we explored this topic with some success and (1/2) + (65/301)Ma report on the results below. (301/14400)Ma 2 Up[2~1] = [17] Examining the decay of U with Ma at large 1 + (130/301)Ma Ma values, we first attempted to fit the nuThe figures show that while both the Pad6 merical results in the asymptotic limit of large approximants do better than the original series Ma with a power law relation of the form U from which they were obtained, the improve- = a(Ma) -b. For this purpose, all the data from ment is not substantial. The next step is to the finite difference results for 50 ~
( M_a /2
-
Journal of Colloid and Interface Science, Vol. 123,No. 2, June 1988
521
MOTION OF A GAS BUBBLE 0.5
I
I
200. We report below the result obtained by fitting the numerical results for Ma i> 75:
I
--
FINITE D I F F E R E N C E
---
EULER
TRANSFORMATION
1.59
1.59
>..~ 0.4
U = log(Ma) + 1.84 "
Results calculated from Eq. [ 18] are plotted in Fig. 4. The values are within 0.6% of the finite difference results for Ma ~ 25 and are indistinguishable from the latter on the figure in this range. Even at Ma = 10, the deviation is only 4%. The excellent nature of the above fit cannot be overlooked in future attempts to construct asymptotic expansions in the limit M a ---~ ~ .
0 m
la.I 0 3 > lad .J eel :~
0.2-
a bJ
,.I t
OJ
-
f,t)
00
[19]
I
50
i
I00
MARANGONINUMBER,
4. CONCLUDING REMARKS
I
150
200
Mo
FIG. 4. Comparison of the scaled bubble velocity, U, as a function of the Marangoni number, Ma, calculated from the finite difference solution, the improved series from the Euler transformation, Eq. [ 18], and the fit for large values ofMa, Eq. [19].
data point at a time from the low Marangoni number end, we found that the constants a and b changed monotonically. In particular, b was found to decrease as more points from the low end were omitted. This suggested that the decay of the bubble velocity with increasing Marangoni number is slower than predicted by a power law. Hence, other possibilities were explored. Prompted by a suggestion from A. Sangani (personal communication), we tried the form U = a/(log Ma + b). The constants a and b were calculated as before by including data for 50 ~< Ma ~< 200. The calculation was repeated by successively omitting one data point from the low Ma end each time. The constants in this case were virtually unaffected by this change in the set of data used, especially when several points were included. The bubble velocities calculated from the three fits, obtained from data for Ma >~ 50, 75, and 100, respectively, were within 1% of each other over the entire range of Ma values from 25 to
It should be noted that while scaled bubble velocities show a decrease with increasing Marangoni number, actual velocities will increase with increasing ~TTo~over the range of Ma values explored. We suggest that Eq. [ 18] be used to calculate the bubble velocity for 0 ~< Ma ~< 20, and that Eq. [19] be used for 20 ~< Ma ~< 200. Unfortunately, experimental data under zero gravity conditions consistent with the assumptions made in this work are not available for comparison. Thus, confirmation of the predictions made here will have to await future experiments on this subject. ACKNOWLEDGMENTS We express our gratitude to Professor A. Sangani of Syracuse University for pointing out the possibility of the logarithmic fit at high values of the Marangoni number. This research was supported by the Microgravity Sciences Division of the National Aeronautics and Space Administration through NASA Contract NAS8-32944 from the Marshall Space Flight Center to Clarkson University, and by the National Science Foundation through NSF Grant CBT-8315048. REFERENCES 1. Subramanian, R. S., Perspect. Comput. 4(2/3), 4 (1984). 2. Young, N. O., Goldstein, J. S., and Block, M. J., J. FluidMech. 6, 350 (1959). 3. Subramanian, R. S., AIChE J. 27, 646 (1981). Journal of Colloid and Interface Science. Vol. 123,No. 2, June 1988
522
SHANKAR AND SUBRAMANIAN
4. Subramartian, R. S., in "Advances in Space Research" (Y. Malmejec, Ed.), Vol. 3, No. 5, p. 145. Pergamon, Great Britain, 1983. 5. Merritt, R. M., Term Project Report, Chemical Engineering Analysis Course (unpublished), Department of Chemical Engineering, Clarkson University, 1984. 6. Acrivos, A., and Taylor, T. D., Phys. Fluids 5(4), 387 (1962). 7. Bratukhin, Y. K., lzv Akad. Nauk SSSR, Mekh. Zhidk. Gaza No. 5, 156 (1975); original in Russian, NASA Technical Translation NASA TT 17093, June 1976. 8. Scriven, L. E., Chem. Eng. Sci. 12, 98 (1960). 9. Dennis, S. C. R., Walker, J. D. A., and Hudson, J. D., J. Fluid Mech. 60, 273 (t973). 10. Roache, P. J., "Computational Fluid Dynamics." Hermosa Publishers, Albuquerque, NM, 1972. 11. Gosman, A. D., Pun, W. M., Runchal, A. K., Spalding,
Journal of Colloid and Interface Science, Vo[. 123, No. 2, June 1988
12. 13.
14. 15. 16.
17.
D. B., and Wolfshtein, M., "Heat and Mass Transfer in Recirculating Flows." Academic Press, Orlando/London, 1969. Brabston, D. C., and Keller, H. B., J. FluidMech. 69, 179 (1975). Shankar, N., "Motion of Bubbles Due to Thermocapillary Effects." Ph.D. Thesis in Chemical Engineering, Clarkson University, 1984. Van Dyke, M., Q J. Mech. Appl. Math. XXVII, 423 (1974). Van Dyke, M., "Perturbation Methods in Fluid Mechanics." The Parabolic Press, Stanford, CA, 1975. Van Dyke, M., "Extension, Analysis and Improvement of Perturbation Series," p. 449. Tenth Symposium, Naval Hydrodynamics, Washington, DC, 1976. Baker, G. A., Jr., in "Advances in Theoretical Physics" (K. A. Breuckner, Ed.), Vol. 1. Academic Press, New York, 1965.
SUBRAMANIAN
Department of Chemical Engineering, Clarkson University, Potsdam, New York 13676 Received May 29, 1987; accepted July 21, 1987 The steady migration velocity of a gas bubble placed in a large liquid body possessing a uniform temperature gradient in the undisturbed state is calculated under conditions of negligible inertia. Proper account is taken of convective energy transport effects, and results are obtained by solving the energy equation using a numerical scheme. The range of utility of an expansion for the bubble velocity valid in the limit of small Marangoni numbers is extended by the use of series improvement techniques. The improved series is good for Ma ~< 20 whereas the original series is useful for only Ma ~ 1. Also, an excellent empirical fit of the numerical results for the bubble velocity is obtained in the limit of large values of the Marangoni number. This result is good for Ma -~-20. Streamline plots are used to demonstrate the existence of a recirculafion region which moves closer to the bubble as the Marangoni number is increased. Isotherms are included to display the presence of an extended thermal wake behind the bubble whose extent increases as the Marangoni number is increased. © 1988AcademicPress,Inc.
Thus, under such conditions, the interface can become the principal driver of fluid motion. A gas bubble will migrate in a fluid if subAn example wherein interfacial tension grajected to a gradient of interfacial tension. The dients can be used to remove gas bubbles from motion is a consequence of the discontinuity a glass melt processed in space is illustrated in the tangential stress across the bubble surby Subramanian (1). face. The traction exerted on the fluid adjacent The fact that gas bubbles will move due to to the interface will cause it to move in the interfacial tension gradients was first demondirection of the pole with the greatest interstrated experimentally by Young et al. (2). facial tension, and by reaction, the bubble will Subsequent experimental work as well as themove in the opposite direction. ory has been reviewed earlier (3, 4), and will Interfacial tension gradients can be mainnot be discussed here. Our attention, in the rained at a bubble surface by arranging to have present article, will be devoted to the theorettemperature or concentration gradients in the ical problem of the steady motion of a gas continuous phase fluid. As a mechanism for bubble in a large liquid body possessing a unidriving the motion of bubbles, when compared form temperature gradient in the undisturbed to a body force field such as gravity, interfacial state. tension gradients become significant at small Young et al. (2) considered the above probbubbles sizes. Furthermore, in "space prolem under conditions wherein convective cessing," a term which describes working with transport effects were completely negligible. In materials and processes in the near free fall their treatment, the general case of a fluid environment aboard orbiting spacecraft, gravdroplet was analyzed. For a gas bubble, one itational effects are relatively unimportant. can make the simplification that the viscosity and the thermal conductivity of the gas are 1 Present address: Martin Marietta Laboratories, 1450 South Rolling Road, Baltimore, M D 21227. negligible compared to the same properties in 1. INTRODUCTION
512 0021-9797/88 $3.00 Copyright© 1988by AcademicPress,Inc. All rightsof reproductionin any formreserved.
Journalof ColloidandInterfaceScience.Vol. 123,No. 2, June 1988
MOTION
the surrounding liquid. In this special case, in the absence of gravitational effects, the result given by Young et al. reduces to u =
a l d a / d T [ IVTo~ [ 2~
[1]
Here, u is the steady migration velocity of the bubble and VT~ is the undisturbed constant temperature gradient in the fluid, d a / d T is the gradient of interfacial tension with temperature, assumed constant, a is the radius of the bubble, and # is the viscosity of the continuous phase, also assumed constant. One direction in which the assumptions of Young et al. can be relaxed is to include the contribution from convective transport. In applications to highly viscous fluids, such as glass melts, the Reynolds numbers will be quite negligible, and inertial effects can be safely ignored. However, the Prandtl numbers in such systems can be very large (10 4 and larger). Thus, the Peclet number for heat transfer can be substantial, and convective energy transport cannot always be neglected. In this problem, the Peclet number when based on an appropriate velocity scale, is commonly designated the Marangoni number. The above problem was considered by Subramanian (3) for small values of the Marangoni number, Ma. He obtained a solution of the governing field equations by the method of matched asymptotic expansions, and gave a result for the bubble speed good to 0(Ma2), U=
1
2
301 --Ma 14,400
2+ . - . .
[2]
In Eq. [2], U is the scaled bubble velocity (U = U/Vo). The reference velocity is Vo = (al da/ dTIl~TZoo])/~. Thus, at zero Marangoni number the result reduces correctly to that predicted by Young et al. Subsequently, Merritt (5) found the next higher term in this expansion, U-
1
2
301 - Ma 2 14,400 13 + 1 - ~ Ma3 + " " ""
513
OF A GAS BUBBLE
[3]
It is interesting to note that, to this order, there appear to be no logarithmic terms in the asymptotic sequence. This is in contrast to the behavior of a similar heat transfer problem associated with a rigid sphere in creeping motion through a fluid which was studied by Acrivos and Taylor (6). In principle, one can carry the calculation to higher orders. However, this process is a formidable exercise in algebra. Furthermore, it is useful to be able to check the utility of the above expressions by comparison with an independently obtained solution of the governing equations. This would permit us to determine the range of Marangoni numbers over which Eqs. [2] or [3] are useful. For the reasons cited above, we have solved the governing equations here numerically using the method of finite differences, and have calculated the bubble velocity up to Ma = 200. Furthermore, we have used series improvement techniques to discover useful information concealed in the coefficients in Eq. [3] and to develop an analytical result which appears to be good up to Ma = 25. In contrast, Eq. [3] diverges from the numerical results beyond Ma = 1. In Section 2, the theoretical model equations are stated, and the method of solution is discussed briefly. In Section 3, results from the calculations are displayed and discussed. We conclude with a few remarks in Section 4.
2. A N A L Y S I S
Consider a gas bubble present in an unbounded liquid body which, in the undisturbed state, possesses a uniform temperature gradient, VToo. Gravitational effects are neglected, and only the motion of the bubble due to the interfacial tension variation resulting from the temperature gradient on the interface is considered. It is assumed that, over the spatial domain of interest, the physical properties of the liquid are constant with the exception ofinterfacial tension. This quantity is assumed to decrease linearly with increasing Journal of Colloid and Interface Science, Vol. 123, No. 2, June 1988
5 14
S H A N K A R A N D S U B R A MA N IA N
temperature. Furthermore, it is assumed that the bubble size does not change appreciably during the migration process. Under the above assumptions, it is possible for the interfacial tension gradient on the bubble surface to achieve a steady representation, resulting in steady migration of the bubble in the direction ofVT~. This is the problem we shall analyze below. As stated earlier, inertial effects are ignored, and the motion of the fluid is assumed to be incompressible and Newtonian. The viscosity and the thermal conductivity of the gas phase are assumed negligible compared to the corresponding properties in the liquid. Thus, only the transport problem in the liquid needs to be treated. In general, the bubble will be deformed from the spherical shape even when inertial effects are ignored. This is because the interfacial tension is not a constant over the bubble surface. However, provided the variation in interfacial tension over the bubble surface is small when compared to the value of the interfacial tension, the capillary number will be small, and the bubble boundary may be assumed spherical with little error. It is convenient to use a reference frame traveling with the bubble. As recognized by Bratukhin (7), an observer riding with the bubble will perceive an unsteady temperature field at all times. However, the important entry which determines the achievement of a steady migration velocity, subject to the assumptions made earlier, is the temperature gradient field. The gradient field, after an initial transient, will achieve an asymptotic steady state. To obtain this, as done by Bratukhin (7), we pose a problem for the difference between the temperature at any point and the temperature at the bubble's equatorial plane at infinity. This difference field will reach a nontrivial steady state after an initial transient, and its gradient is identical to that of the actual temperature field. The following governing equations may be written for the scaled velocity and temperature fields, Journal of ColloM and Interface Science, Vol. 123,No. 2, June 1988
VZv = Vp
[4]
Ma[U + v. X7T] = V2T.
[51
In the above, distances are scaled by the bubble radius, a, velocities by v0, defined earlier, and the dynamic pressure by IVZo~l Ida~ dTI, The scaled temperature is the difference between the actual temperature and that at infinity on the equatorial plane of the bubble, divided by a lVTo~I. The Marangoni number :is defined by Ma = (aVo)/a where o~is the thermal diffusivity of the liquid. The boundary conditions on the fields may be written as follows. Far away from the bubble, the velocity and temperature fields attain their free stream values, v --~ -Uiz ]
[6] as
r--~ ~ .
T --~ z
[7]
Here, iz is-a unit vector in the direction of the temperature gradient VToo in the undisturbed fluid, and z is-distance measured from the equatorial plane in this direction, r is the position vector with the origin at the bubble center. At the bubble surface, Sb, the vanishing of the normal velocity, the balance of tangential stress, and the vanishing of the normal energy flux yield the following boundary conditions. v.nb = 0
]
rib" ~r- tb = --Vs T. t b l XTT.n~ = 0
[8] on
Sb.
[9] [10]
In the above, ,c is the scaled hydrodynamic stress tensor in the liquid. Vs is the surface gradient operator (8), nb is a unit normal vector to the surface, and tb is an arbitrarily oriented unit vector lying in a tangent plane on the bubble surface. Finally, setting the hydrodynamic force on the bubble to zero completes the problem statement. A formal solution for the velocity field may be written immediately from a known general solution of the Stokes equation in spherical polar coordinates by suitable specialization.
M O T I O N OF A GAS BUBBLE
The solution, given by Subramanian (3), contains an infinite set of constants which are coefficients in the expansion of the temperature field on the bubble surface in a series of Legendre polynomials. It is the appearance of these constants in the velocity field which makes the energy equation nonlinear. We solved the above problem for the temperature field using the method of finite differences. In doing so, a new temperature field, T' = T - z, was defined to facilitate the application of the boundary condition at infinity. Also, in recognition of the existence of the sharpest gradients near the bubble surface, the radial coordinate r, measured from the bubble center, was transformed via ~ = log(r). This transformation was employed by Dennis et al. (9), who studied steady heat transfer in forced convection past a solid sphere, and found to be satisfactory. In writing the finite difference equations, central differences were written except for the convective transport terms which were handled using an upwind differencing scheme, as discussed by Roache (10) and Gosman et al. (11). As observed by Dennis et al. (9), the accuracy of the solution can be improved by further replacing the infinity condition by an "improved" one which can be used at a sufficiently large value of the radial coordinate (see also Ref. (12)). By following their procedure which involves the recognition that at large distances from the bubble the velocity field may be closely approximated by the free stream (the Oseen approximation), followed by the solution of the appropriate Oseen equation, the ratio of the temperature fields at adjacent grid points can be obtained. After sufficient testing of this condition against the earlier one, the improved condition was employed in all the calculations at ~ = 3. This corresponds to a location approximately 20 bubble radii from the bubble center. It was verified that using the condition at ~ = 2 gave bubble velocities which were different by no more than 0.5%. In some runs at low values of the Marangoni number, the infinity condition was applied at further radial locations such as 100 bubble ra-
515
dii to confirm that there was no change in the results. An 81 X 81 grid in the two spatial coordinates was employed in all the calculations, and the results checked for accuracy by using a 51 X 51 grid at selected values of the Marangoni number. The latter produced results within 1% of the former over the entire range of values of the Marangoni number investigated here. Typically, on an IBM 4341 computer, it was found that it took approximately 30 rain of CPU time to perform the entire set of computations for a single value of the Marangoni number. Values of the Marangoni number ranging from 0 to 200 were covered. Further details concerning the numerical solution procedure may be found in Shankar (13). 3. RESULTS AND DISCUSSION
Streamlines Typical streamlines from the numerical calculations are displayed in a meridian section in Figs. 1a-1 c for an increasing progression of values of the Marangoni number. The streamlines for Ma = 200 were plotted, but were found to be virtually identical to those for Ma = 100, and hence are not included. Figure 1a corresponds to Ma = 0 and was obtained from the analytical solution given by Young et al. All the plots are given as seen by an observer in a laboratory reference frame. The bubble is moving upward in all the figures. The velocity field is sensitive to the value of the Marangoni number through its dependence on the surface temperature field. The single most striking feature in Figs. 1a - l c is the appearance of an unattached recirculation region in the wake of the bubble. There is perfect fore-aft symmetry in the plot for Ma = 0 since the temperature gradient field on the bubble surface in this case is perfectly symmetric about the equator of the bubble. Such symmetry is lost in the presence of convective energy transport. A streamline on which the scaled Stokes stream function ff = 0 appears in the wake region, and separates the recirJournal of ColloM and Interface Science, Vol. t23, No. 2, June 1988
516
SHANKAR
AND SUBRAMANIAN
I
2-o.ot 1171-o.I I 3 -o.oz II8 o.ool I -oo4 191 o ooz5 I
I -0.005 5-0.08 2 -0.02 6 "0.1 3 -0.04 7 -0.125
4-0.06
-0.06
]18-0.15
0 0.004
I -0005 6-o.o8 [ z -o.oz Ilrl-°.l / 3-0.03
8 0.004 I
4-o.o4 -o.o6
91 0.006 / IOl o.ooe /
C
FIG. 1. (a) Streamlines, M a = 0. (b) Streamlines, M a = 5. (c) Streamlines, M a = 100.
culating region from the rest of the flow field. This may be regarded as a separatrix. The recirculating "eddy" in Figs. lb and lc must close at infinity. The intensity of flow in this region appears to increase with increasing values of the Marangoni number. Also, the recirculation pattern moves closer to the rear Journal of Colloid and Interface Science, Vol. 123, No, 2, June 1988
stagnation point as the Marangoni number is increased. It is interesting to determine whether the above features can be predicted from the small Ma expansion. Subramanian's result for the stream function in a laboratory reference frame, good to 0(MaZ), may be written as
MOTION
OF A GAS BUBBLE
517
The above equation has one real root and a pair of complex conjugate roots for Ma ~<4.76. At Ma = 0.25, the real root is 34.2. 3 This compares favorably with the finite difference result, r0 = 34.6. At Ma = 1, the above equation yields r0 = 8.2 whereas the numerical where solution gives r0 = 9.3. Thus, for small values I2 = - 1 ÷ ~ Ma 2 [12a] of Ma, Eq. (14) can be used to obtain a good estimate of r0. The physical origin of the recirculating eddy 7 13 = 1 ~ Ma [ 12b] can be identified. When convective energy transport is neglected, the resulting tempera29 ture field on the bubble surface is a pure Le14 - Ma 2 [ 12c] 11200 gendre mode, namely Pl(s). The surface gradient of this field appears as the tangential Here, r is the radial coordinate in a spherical stress (in scaled form) and drives the flow in polar coordinate system centered in the bubthe liquid. The basic flow due to this temperble. s = cos 0 where 0 is the polar angle measured from the positive z-direction defined ature gradient is illustrated in Fig. 1a. Crucially important is the fact that the magnitude of the earlier. It may be seen that Eq. [11] correctly pre- velocity in the fluid decays away from the dicts the existence of a dividing streamline or bubble as 1/r 3. When convective transport of energy is inseparatrix for nonzero values of the Marangoni number. Setting ff = 0 leads to the axis of cluded, higher Legendre modes of the surface symmetry as one solution (s = _+1), and the temperature field are excited. For n >1 2, each equation for the separatrix as the other solu- pure mode Pn(s) leads to a flow which decays tion: as r-n as r ~ ~ . Thus, it is clear that flow at large distances from the bubble must be dominated by the P2 mode of the surface temper4r ~I3 -ls ature field. This flow is toward the bubble in the equatorial region, and away from it in the +-~14 (5s 2 - 1 ) = 0 . [13] polar regions. Thus, in the region of the north pole (forward region) of the bubble, this flow When streamlines are plotted from Eq. [ 11 ], is in the same general direction as the basic for representative values of the Marangoni flow due to the P1 mode, and the two simply number, the qualitative features are similar to reinforce each other. In the rear, the two opthose presented in Figs. lb and lc. Further- pose each other. Thus, flow must be away from more, Eq. [ 13] may be used to predict the value the bubble at large distances into the wake and of ro where the separatrix intersects the axis of toward the bubble in its vicinity. This, coupled symmetry. Substituting the coordinates of this with the inward flow at the equator induced point (r0, - 1 ) in Eq. [13], the following cubic by the P2 mode, leads to the observed separated equation is obtained: eddy which moves closer to the bubble as the P2 mode gains strength. As a matter of interest, / 289 r03 + ~7--~ Ma - ~Ma)r~ the P3 mode of the temperature field induces flows which reinforce those due to the P1 mode 261 in the polar regions, and which decay also - r o - 1--9~ Ma = 0. [141 as 1/r 3.
s2,{4 +
3
1
1)s
),5s21,}
Lll
Journal of Colloid and Interface Science, Vol. 123, No. 2, June 1988
518
SHANKAR AND SUBRAMANIAN
Isotherms Isotherms in a meridian plane are plotted for a series of values of the Marangoni number in Figs. 2a-2d. Figure 2d is drawn for the same value of Ma (=200) as Fig. 2c, but over a greater spatial domain to show the extent of the thermal wake region. The plots clearly establish the appearance of a long thermal wake whose extent increases with increasing values of the Marangoni number. The fore-aft symmetry present in the isotherms for the case of pure conduction is destroyed in the presence of convective transport. It may be observed from the plots that the overall temperature variation over the bubble surface is reduced by the inclusion of convec-
tive energy transport. Thus, one might expect that the magnitude of the scaled migration velocity would decrease with increase in Ma. This is indeed the case as will be seen shortly. The condition of no energy flux into the bubble requires that the isotherms meet the bubble surface in a normal direction. At first, some of the isotherms in Fig. 2c appear to violate this requirement. However, this is just an artifact due to the limitation on spatial resolution in this drawing. At large values of the Marangoni number, convective energy transport is dominant in determining the temperature field almost everywhere in the liquid. However, there must exist a thermal boundary layer in the vicinity of the bubble surface where conduction and convection must balance each
Z.Q
2._66
1.0
-I.0
I
T'-2.6 T'-3.0 z~
0
T=-7 FIG. 2. (a) Isotherms, M a = 0, A T = 0.4. (b) Isotherms, M a = 5, A T = 0.4. (c) Isotherms, M a = 200, A T = 0.4. (d) Isotherms, M a = 200, A T = 1. A larger region is shown.
Journal of Colloid and Interface Science, Vo|. 123, No. 2, June 1988
519
MOTION OF A GAS BUBBLE
other. In this boundary layer, the temperature gradient field is adjusted rapidly so that the no flux condition at the bubble surface is met. In making the finite difference calculations, we were careful to check that the spatial resolution is adequate to account for the rapid variation in the gradient in this boundary layer.
Bubble Velocity The principal quest of this work is the calculation of the magnitude of the steady migration velocity of the bubble as a function of the Marangoni number. A comparison of the numerical results for this quantity against the predictions of the expansions for small Ma given in Eqs. [2] and [3] is displayed in Fig. 3a, and over a larger range of Marangoni number values in Fig. 3b. It is seen from Fig. 3a that the small Ma expansions have only limited utility, being good only up to Ma ~< 1.
a
0.7
I
I
f
06 "
~
0.5
-
......
O(M@)
-'-" -"-
PADE" E I / 2 3 PADE' E2./I~
°.°ll
: I
/
F I N I T E DIFFERENCE _ _ _ O(Mo 2) =) "
The figures also contain other approximations which are discussed below. Up to 0(Ma3), the series for the scaled bubble velocity given in Eq. [3] is a power series in Ma. However, while inclusion of the 0(Ma 3) term improves the performance of the series over that when truncation is performed at the previous order, divergence appears to occur quickly as the value of the Marangoni number is increased. This led us to explore the possibility of the presence of a singularity in the complex Ma plane which might be restricting the convergence of a possible power series for the bubble velocity (14-16). We recognized at the outset that such series improvement methods might not yield much information here for two reasons. First, only the first four coefficients are available, and one of them is zero. Thus, a Domb-Sykes plot does not prove useful. Second, in other analogous problems, the asymptotic sequence contains terms involving the logarithm of the pertur-
."
....
J D,I
...'""
I
PADE" E 2 / l ' l
0.4
,,,!
lid ~
I
"-~. ~ . ~ . ~\
J
"~",~. ,
0.3
m
m o w .J
I\ 0.2 t
\
; I
\
0.2
I
W -J
0"41
'~
I
,, FINITE DIFFERENCE --O(Mo 2) ....... O(Ma ~) . . . . PADE" C I/2"1
/"
o
t-~
i
I
I
I
2
MARANGONI
J
3
NUMBER,
1 ~ [ 4
Mo
%1
I
i
I
I
I
I
20
40
60
80
MARANGONI
NUMBER,
I00
Ma
FIG. 3a. Comparison of the scaled bubble velocity, U, as a function of the Marangoni number, Ma, calculated from the finite difference solution, the [1/2] Pad~ approximant, the [2/1] Pad~ approximant, the 0(Ma2) result, and the 0(Ma 3) result. (b) Comparison of the scaled bubble velocity, U, as a function of the Marangoni number, Ma, calculated from the finite difference solution, the [1/2] Pad6 approximant, the [2/1 ] Pad6 approximant, the 0(Ma2) result, and the 0(Ma3) result. A larger range of Ma values is shown. Journal of Colloid and Interface Science, Vol. 123, No. 2, June 1988
520
SHANKAR AND SUBRAMANIAN
bation parameter. If that were to be the case here, series improvement via the usual techniques may not be possible. Thus, we were pleasantly surprised by the results reported below. Pad6 approximants are rational fraction approximants to a power series. The [M/N] Pad6 approximant is defined via f(~) =
PM(E) QN(e)
[15] "
denominator of the [2/1] approximant (Eq. [17]) is located at Ma = -2.32. Those of the [1/2] approximant are at Ma = -3.5 and -6.83. When a singularity is located on the negative real axis, at Ma = -~0, an Euler transformation (Ma = Ma/(Ma + ~o))may be used to map the singularity to the point at infinity. Thus, we attempted this with each of the above guessed locations of a singularity. We found that the results from the [1/2] approximant's singularities performed very poorly. However, using the guess from the [2/1] approximant, considerable success was achieved. The modified expansion for U is
Here f ( 0 is the known power series in e of degree M + N. PM(0 and QN(Oare polynomials in e of degree M and N, respectively. Given a power series to a certain order, various Pad6 approximants can be constructed in a U = 0.5 - 0.1125 Ma + 2.32] straightforward fashion. The properties and use of Pad6 approximants are detailed by - 0.1123(MaMa ,/3 Baker (17). According to Van Dyke (14), "The + 2.32] " [18] approximants with M and N nearly equal posThe bubble velocity calculated from Eq. sess remarkable, though somewhat mysterious, [ 18] is compared with the numerical solution properties of analytic continuation." Hence, these approximants may be used to find the in Fig. 4. The improvement, when contrasted approximate location of singularities in a per- with the performance of the earlier approxiturbation solution by calculating the zeroes of mations, is dramatic. For Ma ~<25, the above the denominator, Qi(e). All the possible Pad6 approximation is indistinguishable from the approximants were constructed from the numerical result. Beyond that point, the nupower series in Eq. [3] and compared with the merical result continues to decrease with furnumerical solution for the bubble velocity. ther increase in the Marangoni number while The best performers were the [1/2] and the Eq. [18] gives a prediction which levels off at [2/1 ] approximants. These are given below and approximately Ma ~ 50. For the benefit of the user, it is good to fit shown in Figs. 3a and 3b: the numerical results for U versus Ma for large Udl/2] values of Ma to a simple equation. In addition, such a fit can lend insight into the analytical 1 + (130/301)Ma [161 structure of an asymptotic expansion for the 2 + (260/301)Ma + (301/3600)Ma 2 fields in the limit of large values of Ma. Thus, we explored this topic with some success and (1/2) + (65/301)Ma report on the results below. (301/14400)Ma 2 Up[2~1] = [17] Examining the decay of U with Ma at large 1 + (130/301)Ma Ma values, we first attempted to fit the nuThe figures show that while both the Pad6 merical results in the asymptotic limit of large approximants do better than the original series Ma with a power law relation of the form U from which they were obtained, the improve- = a(Ma) -b. For this purpose, all the data from ment is not substantial. The next step is to the finite difference results for 50 ~
( M_a /2
-
Journal of Colloid and Interface Science, Vol. 123,No. 2, June 1988
521
MOTION OF A GAS BUBBLE 0.5
I
I
200. We report below the result obtained by fitting the numerical results for Ma i> 75:
I
--
FINITE D I F F E R E N C E
---
EULER
TRANSFORMATION
1.59
1.59
>..~ 0.4
U = log(Ma) + 1.84 "
Results calculated from Eq. [ 18] are plotted in Fig. 4. The values are within 0.6% of the finite difference results for Ma ~ 25 and are indistinguishable from the latter on the figure in this range. Even at Ma = 10, the deviation is only 4%. The excellent nature of the above fit cannot be overlooked in future attempts to construct asymptotic expansions in the limit M a ---~ ~ .
0 m
la.I 0 3 > lad .J eel :~
0.2-
a bJ
,.I t
OJ
-
f,t)
00
[19]
I
50
i
I00
MARANGONINUMBER,
4. CONCLUDING REMARKS
I
150
200
Mo
FIG. 4. Comparison of the scaled bubble velocity, U, as a function of the Marangoni number, Ma, calculated from the finite difference solution, the improved series from the Euler transformation, Eq. [ 18], and the fit for large values ofMa, Eq. [19].
data point at a time from the low Marangoni number end, we found that the constants a and b changed monotonically. In particular, b was found to decrease as more points from the low end were omitted. This suggested that the decay of the bubble velocity with increasing Marangoni number is slower than predicted by a power law. Hence, other possibilities were explored. Prompted by a suggestion from A. Sangani (personal communication), we tried the form U = a/(log Ma + b). The constants a and b were calculated as before by including data for 50 ~< Ma ~< 200. The calculation was repeated by successively omitting one data point from the low Ma end each time. The constants in this case were virtually unaffected by this change in the set of data used, especially when several points were included. The bubble velocities calculated from the three fits, obtained from data for Ma >~ 50, 75, and 100, respectively, were within 1% of each other over the entire range of Ma values from 25 to
It should be noted that while scaled bubble velocities show a decrease with increasing Marangoni number, actual velocities will increase with increasing ~TTo~over the range of Ma values explored. We suggest that Eq. [ 18] be used to calculate the bubble velocity for 0 ~< Ma ~< 20, and that Eq. [19] be used for 20 ~< Ma ~< 200. Unfortunately, experimental data under zero gravity conditions consistent with the assumptions made in this work are not available for comparison. Thus, confirmation of the predictions made here will have to await future experiments on this subject. ACKNOWLEDGMENTS We express our gratitude to Professor A. Sangani of Syracuse University for pointing out the possibility of the logarithmic fit at high values of the Marangoni number. This research was supported by the Microgravity Sciences Division of the National Aeronautics and Space Administration through NASA Contract NAS8-32944 from the Marshall Space Flight Center to Clarkson University, and by the National Science Foundation through NSF Grant CBT-8315048. REFERENCES 1. Subramanian, R. S., Perspect. Comput. 4(2/3), 4 (1984). 2. Young, N. O., Goldstein, J. S., and Block, M. J., J. FluidMech. 6, 350 (1959). 3. Subramanian, R. S., AIChE J. 27, 646 (1981). Journal of Colloid and Interface Science. Vol. 123,No. 2, June 1988
522
SHANKAR AND SUBRAMANIAN
4. Subramartian, R. S., in "Advances in Space Research" (Y. Malmejec, Ed.), Vol. 3, No. 5, p. 145. Pergamon, Great Britain, 1983. 5. Merritt, R. M., Term Project Report, Chemical Engineering Analysis Course (unpublished), Department of Chemical Engineering, Clarkson University, 1984. 6. Acrivos, A., and Taylor, T. D., Phys. Fluids 5(4), 387 (1962). 7. Bratukhin, Y. K., lzv Akad. Nauk SSSR, Mekh. Zhidk. Gaza No. 5, 156 (1975); original in Russian, NASA Technical Translation NASA TT 17093, June 1976. 8. Scriven, L. E., Chem. Eng. Sci. 12, 98 (1960). 9. Dennis, S. C. R., Walker, J. D. A., and Hudson, J. D., J. Fluid Mech. 60, 273 (t973). 10. Roache, P. J., "Computational Fluid Dynamics." Hermosa Publishers, Albuquerque, NM, 1972. 11. Gosman, A. D., Pun, W. M., Runchal, A. K., Spalding,
Journal of Colloid and Interface Science, Vo[. 123, No. 2, June 1988
12. 13.
14. 15. 16.
17.
D. B., and Wolfshtein, M., "Heat and Mass Transfer in Recirculating Flows." Academic Press, Orlando/London, 1969. Brabston, D. C., and Keller, H. B., J. FluidMech. 69, 179 (1975). Shankar, N., "Motion of Bubbles Due to Thermocapillary Effects." Ph.D. Thesis in Chemical Engineering, Clarkson University, 1984. Van Dyke, M., Q J. Mech. Appl. Math. XXVII, 423 (1974). Van Dyke, M., "Perturbation Methods in Fluid Mechanics." The Parabolic Press, Stanford, CA, 1975. Van Dyke, M., "Extension, Analysis and Improvement of Perturbation Series," p. 449. Tenth Symposium, Naval Hydrodynamics, Washington, DC, 1976. Baker, G. A., Jr., in "Advances in Theoretical Physics" (K. A. Breuckner, Ed.), Vol. 1. Academic Press, New York, 1965.