Greek Letters t = distance to surface r = gamma function 4 = void volume fraction Z = surface
Pearson, K.,Ed., "Tables of the Incomplete Gamma Function," McGraw-Hill, New York, N.Y., 1953. Reck, R. A., Prager, S. J., J. Cbem. Phys., 42, 3027 (1965). Reck, R . A.. Reck. G. P.. J. Cbern. Phys., 49, 701 (1968a). Reck, R. A., Reck. G. P..J. Chem. Pbys., 49, 3618 (1968b). Somorjai, G. A., Joyner, R. W., Lang, E., J. Cafal., 27, 405 (1972). Weissberg, H. L., Prager, S . , Pbys. Fluids, 5, 1390 (1962).
Literature Cited Received for reuiew June 26, 1974 Accepted May 2 2 , 1975
Ablow. C M., Motz, H.. Wise, H., J. Cbem. Pbys., 43, 10 (1965). Hegedus, L. L., Prepr.. Div. Pet. Cbem., Am. Cbem. SOC., 18 (3), 487 (1973).
A Theoretical and Experimental Investigation of the Effect of Internal Circulation on the Drag of Spherical Droplets Falling at Terminal Velocity in Liquid Media Ahmed H. Abdel-Alim' and A. E. Hamielec' Chemical EngineeringDepartment, McMaster University, Hamilton, Ontario, Canada
Numerical solutions of the Navier-Stokes equations have been obtained for steady flow around circulating liquid spheres in liquid media in the Reynolds number range 1-50. The results are given in the form of surface vorticity and pressure distributions and drag coefficients. The computed total drag coefficients are in excellent agreement with experimental values for a number of liquid-liquid systems.
Introduction A gas bubble or liquid drop released from rest in a viscous fluid under the influence of gravity will accelerate initially. After a finite time the drop will acquire a constant velocity and a particular shape. At this point the force of gravity will be balanced by buoyancy and the drag force of the exterior fluid. At the fluid-fluid interface there will be an equilibrium of normal forces. The forces acting inward will be due to the dynamic stress and static head of the exterior fluid and interfacial tension. Those acting outward will be due to the dynamic stress and static head of the interior fluid. If the drop is spherical, all the forces will lie on a radial line and the interfacial tension force will be the same on all parts of the surface. For this case the interfacial tension force will influence the flows of neither interior nor exterior fluids. Friction on the surface is transmitted across the fluidfluid interface setting the interior fluid in motion. Circulation of the interior fluid decreases the boundary layer thickness of the exterior flow and reduces the angle of flow separation with a concomitant reduction in drag and increase in terminal velocity. The circulatory motion in bubbles and drops has been observed visually by several investigators (Savic, 1953; Garner and Hammerton, 1954; Garner and Skelland, 1955; Johnson and Braida, 1957). These observations showed that below a critical drop size, peculiar to each system, motion in the interior did not exist. This is often attributed to the accumulation of surfactants on the interface, which influences smaller droplets to a greater degree. Research Department, Imperial Oil Ltd., Sarnia, Ontario, Canada. 308
Ind. Eng. Chem., Fundam., Vol. 14, No. 4, 1975
The literature contains a voluminous amount of terminal velocity data for drops and bubbles moving through viscous, incompressible media under the influence of gravity. An excellent review article is that by Haberman and Morton (1953). In a recent study (LeClair et al., 1972), the complete Navier-Stokes equations for both interior and exterior flows were solved numerically for spherical water drops falling a t terminal velocity in air. Surface velocities and pressure distributions and drag coefficients were calculated. There was excellent agreement between calculated and measured surface velocities and drag coefficients. Such a gas-liquid system is characterized by very small viscosity and density ratios. No calculations have so far been reported for density and viscosity ratios of the order unity. The purpose of the present investigation was therefore to compute velocities and drag forces for liquid-liquid systems by numerically solving the Navier-Stokes equations for a range of viscosity ratios and density ratios of the order unity and to compare computed drag coefficients with available experimental values for a range of intermediate Reynolds numbers.
Theory Details of the theory involved in the calculation of flow fields around circulating spheres are given by LeClair (1970). A brief summary is given in this article. Using the dimensionless variables rr = r / a , W = q/U,a2, (' = {a/U,, where Jr is the stream function and ( t h e vorticity, and dropping the superscript primes, the NavierStokes equations of motion for viscous, incompressible, axisymetric flow can be written in spherical coordinates and in nondimensional form as (LeClair, 1970)
Table I. Numerical Parameters Used in Computation .I1 AY
Re A
'Y 'YIAV
30
Number of angular mesh points Number of radial mesh points Reynolds number Radial step size
65
50 25
55
10
40 50
0.05
Viscosity Ratio Number of radial mesh points inside
E2* = {r sin 0
60 5
1
System I 0.0995. 0.301, 0.554
System I1 0.266, 0.708, 1.4
41
(1) and the surface velocity distribution using the relationship
(
a -
ar r sin 0
)] sin 0 = E2({r sin 0)
In the case of a circulating fluid sphere these equations are applied to the interior as well as the exterior fluid with the appropriate parameters. An exponential step size in the radial direction was used in the exterior flow field with the substitution r = ez, and the use of equal intervals in z . For the interior flow field a constant step size in the r-direction was used. With these conditions the Navier-Stokes equations of motion can be written for the exterior flow field as
e* sin 8 = e Z Z E 2( roeZ sin 0)
(3)
where k ( a , 8 ) = (p(a,0) - pa)/(1/2 p u * _ ) . T h e form and friction drag coefficients were computed using the relationships
CDP = CDF = Re
J"k ( a , 0 )sin 28 dd
1"[(E o
(10)
u ) sin 0 - ar
]
2 (A)sin 0 = E* ({ir sin 0) ar r sin 0 hr sin 0 = E2*i
(5)
46)
T h e following boundary conditions are to be satisfied. (1) Far from the drop at r = r,, z = z,, undisturbed parallel flow is assumed to occur \ko = (1/2)ez- sin2 8
=0
(2) Along the axes of symmetry 0 = 0, P *i
CDT
= CDP
+ CDF
(12)
Details of the numerical analysis regarding appropriate step sizes, etc., are given elsewhere (LeClair, 1970).
and for the interior flow field as
So = {o =
(9)
(11) where u and u are radial and tangential velocity components; the total drag coefficient is given by
(23I
(0
1
ui(r = l , 0 ) = - sin 6 iaa:i)r=i
(2)
=
{i
=0
(3) There is no material transfer across the fluid-fluid interface. r=l,z=O;\ko=\ki=O (4) There is continuity of tangential stress a t the fluidfluid interface.
TO(r,B)= 7i(r,0) from which, a t r = 1
T h e frontal stagnation pressure was computed using the relationship
The surface pressure distribution was computed using the relationship
Experimental Section Experimental techniques and data on drag coefficients for liquid-liquid systems may be found elsewhere (Hamielec, 1961). The data of interest in this study include total drag coefficient for two liquid-liquid systems. These are: System I: cyclohexanol (continuous phase)-water (dispersed phase); Re range, 1-10; viscosity ratio range (x = , L L ~ / ~ o0.1-0.6. ), System 11: n-butyl lactate (continuous phase)-water (dispersed phase); Re range, 5-50; viscosity ratio range, 0.25-1.5. For both systems the density ratio was almost unity. Different viscosity ratios were obtained by adding small quantities of carboxymethyl cellulose to the dispersed phase. This did not cause a significant change in density or in interfacial tension. Results and Discussions Details of the numerical parameters employed are summarized in Table I. The accuracy of the solutions was examined by varying both the radial and angular step sizes by factors of 2. This test indicated that the solutions are probably accurate to better than 1%for both local and integrated values. Solutions were assumed to have converged when values of vorticity and stream function a t all mesh points changed by less than 0.001 per iteration. Figure 1 shows the variation of the surface velocity with polar angle for different Reynolds numbers a t a viscosity ratio of 0.708. The variation of surface velocity with polar angle for different viscosity ratios at a Reyonlds number of 50 is shown in Figure 2. The curve for zero viscosity ratio (gas bubble) was calculated previously (Hamielec et al., 1967). Figure 3 is the angular variation of surface pressure at different Reynolds numbers at viscosity ratio of 0.708. Figure 4 shows the angular surface pressure distribution a t Re = 50 for different viscosity ratios. The curve for zero viscosity ratio was calculatInd. Eng. Chem., Fundam., Vol. 14, No. 4, 1975
309
z
I
I
X:0.708
08-
10
,Re=50 Re = 25
s
Re
:
5 06
04
iiE!5Yh=o
w
go2
ao
a2
i
POLAR ANGLE
w
M
POLAR ANGLE
00
B
Figure 1. Variation of surface velocity with polar angle for different Reynolds numbers at a viscosity ratio of 0.708.
i? - 0 2 -04
30
0
60
120
90 POLAR ANGLE
I80
150
Figure 4. Variation of surface pressure with polar angle for different viscosity ratios at a Reynolds number of 50.
POLAR ANGLE
Figure 2. Variation of surface angular velocity with polar angle for different viscosity ratios at a Reynolds number of 50.
0
I
NUMERICAL SOLUTION
OZC
01
02
04
03
05
06
X
Figure 5. Variation of surface velocity at 0 = */2 with viscosity ratio at Re = 1.0. Continuous line is the Hadamard-Rybczynski prediction.
60-
n - BUTYL LACTATE -WATER 50
-N W R I C A L X
CDT
0
SOLUTION
X:14 X.0708
30 -
20
‘(h-xz
10-
00 0
*------A-
I
1
100
20 0
I
1
30 0
40 0
500
Re -06
0
1 30
I
I
60 90 POLAR ANGLE
1 120
I 150
Figure 3. Variation of surface pressure with polar angle for different Reynolds numbers at a viscosity ratio of 0.708.
ed previously (Hamielec e t al., 1967). In Figure 5 the surface velocity at 0 = 7~/2variation with viscosity ratio at Re = 1.0 is compared with the prediction of the HadamardRybczynski theory, which is a limiting solution (Re 0) r = 1,8=5)
=f(l+,) 1
-
where x is the viscosity ratio. T h e disagreement is probably a result of the fact that the limiting solution is in error at a Reynolds number as large as unity. The computed drag coefficients for different Reynolds numbers and viscosity ratios are compared with the experi310
Ind. Eng. Chem., Fundam.. Vol. 14, No. 4, 1975
Figure 6. Total drag coefficient vs. Reynolds number at different viscosity ratios for the system n-butyl lactate-water. mentally measured coefficients (Hamielec, 1961) in Tables I1 and 111. Table I1 for the cyclohexanol-water system also shows the drag coefficients for Re = 1.0 calculated using the Hadamard-Rybczynski theory. cDT =
8
(-)3 x + 2
Re x + l
where CDT is the total drag coefficient. The other limiting case (high Re) was treated by Harper and Moore (1968). They calculated t h e drag coefficient assuming the validity of the boundary-layer theory 48 CDT = - (1 1 . 5 ~ ) Re
+
Table 11. Predicted and Experimental Drag Coefficients for the System Cyclohexanol-Water Viscosity ratio
Re
C DF
CDP
GT
Etptl CDT
0.0995 0.301 0.554 0.0995 0.301 0.554 0.0995 0.301 0.554
1.o 1.o 1.o 5 .O 5 .O 5 .O 10.0 10.0 10.0
11.59 12.87 13.29 2.78 3.18 3.30 1.65 1.90 2.04
5 -91 6.33 6.72 1.61 1.63 1.82 0.85 0.99 1.07
17.50 19.20 20.01 4.39 4.81 5.12 2.50 2.89 3.11
18 19 20 4.6 5 .O 5.2 2.8 3.1 3.3
Hadamard-Rybczynski 8 3x+2 c,, = Re x + l 16.7 17.9 18.9
Table 111. Predicted and Experimental Drag Coefficients for the System n-Butyl Lactate-Water Harper and Moore Viscosity ratio
Exptl
Re
CD F
' D P
0.266 0.708 1.40 0.266 0.708 1.40 0.266
5 .O 5 .O 5.0 25.0 25.0 25.0
0.708
50.0 50.0
3.05 3.24 3.26 0.88 1.05 1.12 0.57 0.73 0.78
1.49 1.66 1.76 0.43 0.55 0.57 0.27 0.38 0.45
1.40
50.0
This expression is only a first approximation. Second-order terms may be included to give more accurate estimates. If the Reynolds number is not high enough, however, these terms lead to erroneous predictions. In Table I11 the data for Re = 50 are compared with Harper-Moore first approximation. The poor agreement is attributed to the fact that a Reynolds number of 50 is not high enough. In fact when second-order terms were included the theory predicted negative drag coefficients clearly indicating its inadequacy for the Reynolds number range studied here. Variation of total drag coefficient with Reynolds number is shown in Figures 6 and 7 for the n-butyl lactate-water system and the cyclohexanol-water system, respectively. In Figure 7 1/Re is used instead of Re since the variation of C[)T with Re in this range is rather large. Experimental points are also shown for different viscosity ratios. The data in Figures 6 and 7 indicate that the variation of CDT with viscosity ratio is not great. However, there is a slight but clearly noticeable increase in the drag coefficient with viscosity ratio. I t was pointed out earlier (LeClair et al., 1972) that the drag coefficients for rigid and circulating spheres do not greatly differ. Agreement between experimental drag coefficients and calculated ones is excellent. This suggests that numerical techniques may be used accurately to predict velocity distributions for viscous flow around circulating spherical droplets moving a t terminal velocity in a liquid medium. These profiles could then be used to predict rates of forced convection heat or mass transfer. To present the numerical results of this study in a more usable form they were fitted in an empirical equation that gives the variation of the drag coefficient with Reynolds number and viscosity ratio 26.5
[ (1.3 +
CDT = -
(1.3
x)'
- 0.5
+ x)(2 + x )
1
4.54 4.90 5 -02 1.31 1.60 1.69 0.84 1.ll
1.23
CDT
CD, = 48 (1 + 1 . 5 ~ ) Re
4.5 4.9 5 .O 1.45 1.75 1.8 0.9 1.2 1.25
1.34 1.98 2.98
25 CYCLOHEXANOL- WATER
I
-NUMERICAL
SOLUTION
5\
Figure 7. Total drag coefficient vs. Reynolds number a t different viscosity ratios for t h e system cyclohexanol-water.
For gas bubbles ( x
-
O), the above correlation reduces to
12.1
CDT = Re0.74 This information should be of some interest in the design of liquid-liquid contactors for such chemical engineering operations as liquid-liquid extraction and direct contact heat transfer. Ind. Eng. Chem., Fundam., Vol. 14, No. 4, 1975
311
Nomenclature
T
a = sphereradius Ct)p = form drag coefficient CDF = friction drag coefficient Cryr = total drag coefficient = CDP CDF E 2 = differential operator h = dimensionless pressure p = pressure r = dimensionless radial coordinate Re = Reynoldsnumber u = radial velocity (dimensionless) u,uo = tangential velocity (dimensionless) x = viscosity ratio = wi/wo z = defined asln r
+
Greek Symbols \Tr = dimensionless stream function = dimensionless vorticity 0 = angular spherical coordinate (angle from forward stagnation point) = viscosity p = density
