- Email: [email protected]
A simulation of surface tension driven coalescence
A Simulation of Surface Tension Driven Coalescence YAEL HIRAM l AND AVINOAM NIR 2 Department o f Chemical Engineering, Technion-lsrael Institute o f Technology, Haifa 32000, Israel Received November 30, 1982; accepted February 7, 1983 Coalescence of particles in a continuous fluid occurs quite commonly, e.g., in engineering, physics, and biology. The motion is induced by the tendency of the interface to reduce its area in order to minimize the interracial energy. In this communication we present a numerical solution for the equation of motion describing the coalescence of two particles using an integral equation representation for the surface velocities. The results show a similar behavior for all particles to fluid viscosity ratios X. Processes such as polymer sintering (X ~ ~), biological particles fusion (X ~ 1) and bubble coalescence (X ~ 0) share a similar dynamics when simulated by surface tension driven flows. It is shown that Frenkel's constant power law for the evolution of the neck radius and its velocity of growth does not generally hold during the bulk of the coalescence process. The simulation predicts well the available experimental data describing the isothermal sintering of polymer particles. INTRODUCTION
Many physical processes involve coalescence of dispersed elements in a continuous medium to form bigger particles. Bubbles and droplets coalescence (1, 2), merger of membrane-coated subcellular vesicles in biological systems (3, 4), and sintering of amorphous polymer particles (5) are but a few examples of the phenomenon. The coalescence dynamics is thus an important step in a variety of processes such as multiphase chemical reactions, sedimentation, emulsification, liquid-liquid extraction, combustion, spray drying fermentation, biological mass transfer and uptake, emulsion polymerization, and polymer processing. These differ in the rheological properties of the fluids involved, the scale of the particles participating, the interfacial properties of the phases, and the time scale during which a complete coalescence is obtained. Nevertheless, in all these cases the driving force for coalescence is the excess in interfacial energy Present address: Department of Chemical Engineering, Imperial College, London SW7 2AZ, England. 2 To whom correspondence should be directed.
which exists once fusion between particles commences. Motion is induced by the tendency of the interface to reduce its area in order to minimize the interfacial energy. Real fluids and real interfaces have complex rheological properties. The interface may exhibit viscous, viscoelastic, or even elastic behavior with both spatial as well as temporal variations, depending on the surface chemical composition and the interaction with the bulk. Yet a simple flow model can serve to unify the approach to the various cases and study the important parameters involved. We consider here processes in which the flow can be assumed to be a motion of Newtonian fluids where the interface is characterized by interfacial tension, a. Furthermore, since the typical velocity is of 0(a/ (~ti + #0)), where ~i and ~0 are the viscosity of the fluid particles and the continuous phase, respectively, it is assumed that the typical Reynolds number based on this velocity and the particle dimension is small enough such that inertia effects can be neglected. Indeed this assumption is valid for the case of polymer sintering processes due to the high viscosities involved and for the case of bio-
462 0021-9797/83 $3.00 Copyright© 1983by AcademicPress,Inc. All rightsof reproductionin any formreserved.
Journalof Colloidand InterfaceScience, Vol,95, No. 2, October1983
SIMULATION
OF COALESCENCE
463
solution of the flow problem. Such an approach is not simple since one has to solve the equations of motion simultaneous with the evolution of the curved interface during the entire merger of the two droplets to the ultimate single spherical particle. In principle two intersecting spherical droplets can be described by a conformal mapping using bipolar coordinates. Unfortunately, this technique involves complex eigenfunctions and the general axisymmetric Stokes streamfunction is yet to be produced. Furthermore, such a mapping involves geometrical assumptions at the intersection contour similar to those introduced by Frenkel (6). Ross et al. (9) have attempted a computer simulation for the sintering of a continuous row of particles using a finite element method for the solution of the equation of motion. Although they have encountered accuracy problems in the initial stages, their finding is most important, where Ro is the radial dimension of the neck namely, Ro does not change with t through between the particles and t is the time. A is a single power law but rather the power varies a constant depending on the physical prop- monotonically during the neck growth. Roerties. Frenkel's result was not based on a senzweig (10) has attempted a finite element solution of the equations of motion, how- solution to the coalescence problem of two ever, it was widely used as a basis for studies particles. He suggested an overall power law involving sintering of spherical amorphous correlation similar to Frenkel's equation, polymer particles. Furthermore, the result even though his results indicated a monowas taken to apply even outside its originally tonic change of power over the limited region assumed region of validity, e.g., beyond the of contact neck radius which he considered. initial stage of coalescence. Two basic deviaIn this communication we solve the equations from Frenkel's (6) model were ob- tion of motion numerically using the integral served. Kuczynski et al. (7) found that the equation representation for the surface vedynamics of coalescence varies with temper- locities. This approach enables a solution for ature where both A and the power of t are the evolution of the interface during the conot constant. They have attributed this to the alescence or sintering process without the neglected non-Newtonian viscoelastic prop- necessity to evaluate the details of the flow erties of the polymeric material. Yet, their field within or out of the particles. We show evidence still supported the fact that the neck that the Frenkel model is generally not apsize grows as a power of time. Similar find- plicable to the entire process of neck growth ings were recently reported by Rosenzweig even in the case of Newtonian fluids with a and Narkis (8) who also obtained small os- constant surface tension. Furthermore, cillatory deviations from the logarithmic de- boundary layer considerations indicate that pendence (see Figs. 1 and 5 in this reference). in the initial transitional stages of coalesOf course, a direct way to check the model cence, deviations from Frenkel's model are and the extent of its validity is an explicit always possible. The general equations are
logical particle fusion in view of their small dimension. It is probably only roughly acceptable for atmospheric droplets merger and grossly violated in the case of bubbles coalescence. The sintering of polymer particles is an example which has been extensively studied experimentally and theoretically (6-9). An approximation based on a surface tension driven flow model for this process was proposed by Frenkel (6). He assumed that the interracial shape of two spherical droplets, during initial stages of coalescence, is described by the intersecting joint of two spheres. Approximating the energy dissipation from geometrical consideration and balancing it with similar estimate for the surface energy change Frenkel obtained the power law R o = A t 1/2 [1]
Journal of Colloid and Interface Science, Vol. 95, No. 2, October 1983
464
HIRAM AND NIR
given in the following section. We then describe the numerical procedure, and the resuits are given and discussed. INTEGRAL EQUATION FORMULATION
0/'/i
= 0
[5]
where n,- is a unit normal pointing outward and cr is a constant surface tension. The velocity at a point x on the surface S can be related to the velocities and forces at all points y at the interface via an integral equation (see Rallison and Acrivos (11), Ladyzhenskaya (12))
. [8]
KijkujnkdSy
1 -
2 fsr &~(V. n)nflSy.
[9]
The only parameter in Eq. [9], X, is the key which enables simulations of different physical processes as polymer sintering (k ---, oo), biological particles fusion (X ~ 0(1)), and bubble coalescence (X ---, 0). The approach to a solution of Eq. [9] is a numerical evaluation of the surface velocities with a quasistationary interfacial shape and a subsequent evolution of the surface. When the problem has an axis of symmetry it is possible to eliminate the angular dependence in Eq. [9]. The velocity and surface unit normal vector have only axial and radial components, uz, ur, nz, and nr, respectively, and [9] is reduced to the form Uz(Z)+ 2 11 -_X) ~ f t-t (Kzzuz(~) + Kzrur(~))d~
(1 + k)ui(x) + (1 - X)
2l f f l (Jzznz(() + J~nr(())V" n(~)d(
)< ~ s, giik(X -- y)uj(y)nk(y)dS,
1 fs, J0(x - y)A~rjk(y)nk(y)dSy 87r/~0
yl 3
bli(X ) At- 2 ~
On.._~j
_
Ix -
[31
where /x indicates a difference operation at the interface (inner minus outer), and a balance of stresses
1
(Xi -- Yi)(Xi -- Yi)
t~ii
Jij(x - y) - Ix ~ yl +
[21
within the droplets, and with a similar expression where tZo replaces ~i in the outer fluid. u, a, and p denote the velocity, stress, and pressure fields, respectively. The interfacial conditions are continuity of velocities, AU/= 0 [4]
Atrijn i = --tr OXi ni
[71
Equation [6] can be rendered nondimensional by scaling u with U = a/2a-~0(1 + X) and x by a. Then, using condition [5] for the difference in stress distribution we have
where
[Oui Ouj] alj = -P~ij + ~ i ~ x i + ~x,]
3 (X i -- yi)(xj - Yfl(Xk -- Yk) Ix - yl5 47r and
Consider two Newtonian droplets, each of radius a, coalescing in a continuous fluid. The fluids obey the Stokes equations Oail = 0 ;
Kijk(X -- y)
1 -Xf t u~(z) + 2 1 + k d-i (Krrur(~) -t- Krzuz(~))d~ [61
where X = ~i/#0. Here the kernels K and d derive from the fundamental solutions to the Stokes equations Journal of Colloid and InterfaceScience, Vol. 95, No. 2, October 1983
_
--
1 ff
"2 -t (J~rnr(~) + Jrznz(())V" n(~)d~ [101
where L = 2l is the length of the merging particles along the axis and
SIMULATION
V'n=
l, d z 2 l k
OF COALESCENCE
+ \-~z] ]
1 R
1+
\dz] ]
[111
with the surface given by r = R(z). The axisymmetric components of K and d are given by Youngren and Acrivos (13). NUMERICAL
SOLUTION
The integral Eq. [9] was solved using the following procedure. The initial shape was assumed to be two spherical droplets intersecting each other at a small joint neck. Surface points were distributed along the interface. The spacing between points was chosen such that a denser distribution was applied where derivatives and curvature were large. The equation was reduced to a set of simultaneous linear algebraic equations for the components of the velocities at the surface points. These velocities were used, in turn, to evaluate the normal displacement of the surface segments during a small interval of time to form a new interfacial shape. The procedure was applied repeatedly from the initial shape until all surface motion diminished. We have used an interactive language (APL) which is very suitable for operations on systems of simultaneous algebraic equations and is even more attractive for use in solutions of evolutionary processes since checks and changes can be made and intermediate surface shapes can be plotted as the solution progresses. The number of surface points, N, on one droplet was chosen as a compromise between the will to increase numerical accuracy and the tendency to reduce the size of matrices to be inverted. The sequence of calculation was as follows. First and second surface derivatives, the components of the unit normal, and the surface curvature at each point were evaluated using a straightforward second-order finitedifference scheme for unequal intervals. For
465
a given point x those were used to calculate the tensors K and d at the surface points y. It should be noted that d has a logarithmic singularity as y ~ x except when x is the point on the axis of symmetry where d is of 0(IX -- y[-l/2) The integral involving d was calculated using a trapezoidal procedure. The surface segment with the singularity was excluded while its contribution was evaluated analytically and added. Repeating this procedure for each point x on the two droplets, a set of 2 x ( 2 N - 1) simultaneous equations for the axial and radial velocities components at all 2N - 1 surface points was produced. Using the fore-and-aft symmetry reduced the rank of the problem to 2 × N. This set was readily inverted by the APL "domino" operation to yield the desired surface velocity. Next the surface evolution was evaluated by considering the motion of each surface point normal to the surface segment during a small time interval dt, i.e., the position increment was dx = (u. n)ndt. The fluid points move of course also tangentially, hence, after the increment dt a new set of fluid points is at hand. In order that the mathematical surface points remain approximately in an appropriate spatial separation after each step, especially in the narrow neck region, they were relocated using surface smoothing by spline functions (14). The magnitude of the time step was chosen small enough to ensure stability of the evolution scheme. At <~AZmi,/ Un max proved stable enough in all cases investigated as is also the case in hyperbolic evolutionary problems. Accuracy was measured by requiring that the volume of the droplets does not change. This is equivalent to satisfying the equation of continuity. In most cases investigated the total volume change during the entire process was a few percent (see Fig. 1). When X ~ 0, however, the changes in volume become large indicating, as was also pointed out by Rallison and Acrivos (11), that multiples of homogeneous eigensolutions of Eq. [9] appear in the scheme. Indeed the cumulative Journalof Colloidand InterfaceScience,Vol. 95, No. 2, October 1983
466
H I R A M A N D NIR
0.00
-
0.04
Error ( -AV -) V -0.08
-0.12
-0,16
_0] I
-0.
1
0
I
-0.6
I
-0.2
I
0.2
0.6
1.0
1-X_X I*X FIG. 1. Total relative change in volume V during the simulation process.
volume change for the case ), = 0.1 was as high as 22%, while for the case ), ~ oo it was an order of magnitude lower. RESULTS A N D DISCUSSION
The evolution of the interface shows a similar behavior for all viscosity ratios. Figure
'In,1 .
. . . . . . . .
2 depicts the dependence of the neck radius, R0, and its growth velocity., Ugo, on time for X ~ ~ and ), = 1. It shows three regions of evolution. An initial transitory period, an intermediate stage and a final decay of growth. For all X investigated, 0.1 < X < ~ , this qualitative dynamic similarity remains. It suggests that a wide spectrum of physical processes differs only slightly when they are driven by surface tension and retarded by viscous dissipation. Indeed, an indication that X is not a parameter with singular behavior is already given in Eq. [9] where the actual parameter appears to be (1 - ),)/(1 + X) and varies in the interval (-1, 1) when X changes from zero to infinity. A further evidence that the evolutionary characteristics of the process are similar at all ), is given in Fig. 3 which shows the inward velocity of the surface on the axis of symmetry, i.e., half the shrinkage speed of the coalescing doublet. At the initial steps the numerical scheme is inadequate to describe the evolutionary process. In the vicinity of the neck there exists a boundary layer where extremely high surface curvature and velocities are involved. This stage is short and its effect on the total 101
I
i
~ , i,,,,
I
Transitory
'"Fre n k ¢ l " ~
t i 1,1,,
t
t
' "Frenkel " "h..
stQ~ 10 o
i , ~1,1~
"
stage
-Transitory ' ~
stage
i
URo
~
5
10o
~
X=I i
10"1 10 -2
t
I
I
, , ,,,I
'
'
' ' ' ''
10 q
I(3' 10 0
I
10-3
I
I
I I I Ill
I
I
I
10.2
I I Illl
I
10-1
I
I
I I I I~
10 0
t a/U
t a/U
FIG. 2. Dimensionless evolution of the neck size, Ro, and its radial velocity, URo, as function of time. (a) X ~ 0% (b) X = 1 (a, initial droplet radius; U = a/2~rtz0(1 + X)).
Journal of Colloid andlnterface Science, Vol. 95, No. 2, October
1983
SIMULATION
0.0
dicate that, even at this early stage, deviations from Frenkel's 1/2 power law [1] may occur. In the inner region the variations in the axial, z, and radial, r, directions are of 0(e) where ~ is a small measure of the gap between the droplets (e.g., the gap before coalescence started or the neck size itself). Since the radial and axial components of the velocity are of different orders of magnitude the continuity equation reduces to
-0.5 u& U_I.(
Y -2£
i
O0
0.15
i
0.30
10 r ar (rUr) = 0. i
t a/U
467
OF COALESCENCE
0.45
0.60
FIG. 3. T h e d i m e n s i o n l e s s c h a n g e o f the v e l o c i t y o f the surface p o i n t o n t h e axis o f s y m m e t r y as f u n c t i o n w i t h t i m e a n d ), ( U = a/27r/~0(1 + X); a, i n i t i a l d r o p l e t radius).
dynamics is transitory. The inaccuracy of the numerical results in this region reflects the fact that for any refinement of the spacing of surface points there is always a short period during which the boundary layer is of a smaller size. Some useful information, however, can be Obtained from a boundary layer consideration. At this stage the surface can be separated into three distinct regions. An outer region where the droplets remain spherical having almost the original size, a narrow inner region including the connecting neck where the curvature and the normal velocity attain extremal values, and an intermediate region separating the two in which the surface gradients are high but the curvature passes from negative to positive values. The exact features of the inner region grossly depend on the relative position of the droplets when coalescence started. It is therefore of no general interest to obtain exact perturbation solution to this boundary layer. We shall only point out the characteristics of the inner region which are of interest in as much as they in-
[12]
Hence the flow outside the droplets, to a first approximation, is a source flow C1 Ur-- - r
[13]
where C1 is a constant source strength. At the neck we have of course URoR0 = C~ which apparently agrees with Frenkel's hypothesis. However, it will become clear shortly that C~ itself depends on R0 as well. From the r component of the Stokes equations it is clear that the pressure in the gap near the neck must be constant. It follows that the normal stress component outside the interface is dominated by viscous terms and is of the dimensionless form O'rr ~--- --2 --Cl /.2 "
[14]
Neglecting the contribution from the interior, and with a surface given by r = R(z), the normal stress balance is approximately of the form 0"rrnrn r =
-
-2
CI R2(I + R '2)
R" 1 + (1 + R'z) 3/2 R(1 + R'a) 1/2
[15] where R' denotes dR/dz. Equation [5] can be integrated once to become Journal of Colloid and Interface Science, Vol. 95, No. 2, October 1983
468
HIRAM AND NIR
(1
+ R'2) 1/2 = - C 1
+ C2R2
[16]
where C2 is a constant. Hence at the neck where z = R' = 0 we find C1 = C2R 2 - 1
[17]
where C2 must be determined by integrating Eq. [ 16] and matching with the outer regions. It is sufficient here to indicate that in view of Eq. [17] deviation from URoRo = constant are evident even at this early stage and a Frenkel 1/2 power law may be violated. Indeed, the boundary layer analysis applies only to a very limited period of time. After this transitory period it is expected that the interface shape will not depend on the specific features of the early stages, and will further evolve in a quasi-steady manner to complete coalescence. It is at this stage that application of a numerical scheme to solve Eq. [9] for the entire shape is appropriate. It is also at this stage that experimental results concerning sintering of polymeric particles are available for comparison (10). In the intermediate stage only a small change of surface area occurs, however, it is at this period where most of the neck growth is evident since dissipation is also small. Figure 2 shows that the neck evolution process does not obey Frenkel's constant power law for an extended period of time. This result was also suggested by Ross et al. (9) who applied a finite element technique to study the evolution of necks in a continuous array of coalescing spheres. In order to better illustrate the above conclusion we attempt to apply Frenkel's power m to a time interval within the intermediate region. If R o a t m then the decline in the velocity should obey URoOttm-1. F o r the limited interval of time it is possible to find m depending on the specific h. These are shown in Fig. 4. It is clear, however, from our results (Fig. 2) that the velocity can be only roughly correlated by a constant power m. This power depends strongly on the choice of time interval in which the correlation is considJournalof ColloidandInterfaceScience,Vol. 95, No. 2, October 1983
ered. Rosenzweig (10) obtained results for a finite element simulation of the sintering of two Newtonian fluid spheres (k --~ oo). He studied the period during which the neck varies over the interval 0.6 ~< R0 ~< 0.8 for which he obtained a power m = 0.3, markedly below the m = 0.5 Frenkel value as well as below our results for X ~ ~ . A careful differentiation of Rosenzweig's plot shows no existence of a constant power growth. In addition, the low value obtained for m merely reflects the fact that his results were obtained in a time interval almost beyond the pseudo-Frenkel stage and support our finding that m is decreasing with evolution. Beyond the intermediate stage comes the final period of coalescence where most of the changes in surface area occur accompanied by strong viscous dissipation and damping of the evolution. Variations in surface curvature are mild and the particle assumes its ultimate spherical shape following a slow first order exponential decay of the motion. Figure 5 depicts a semilogarithmic plot of the decay (1 - R o / R o nn~) and of the neck velocity for the entire process for two different X's. For all practical purposes it is clear that
0.6 O
0.5 Frenket's Power m
0.4
0.3
°:o
I
5
I
I
I
3o
oo
FIG. 4. T h e dependence of Frenkel's power, m, on the viscosity ratio, ~,.
469
SIMULATION OF COALESCENCE 101
101 k=30
3
~
a-finradius aL
k=0.3
3~'2a-radius final
10c
10 0
"
~ffl
0
~
~
1-, Ro
.
I
I
1
0.1
0.2
0.3
t a/U
161
0.4
c~.1
0
o'.2
0'.3
t a/U
o'.4
d.5
0.6
FIG. 5. Dimensionless evolution of(1 - Ro/Ro~n~) and uR0 as a function of time. (a) X = 30, (b) X = 0.3 (U = a / 2 r ~ ( 1 + X); a, initial droplet radius).
1 - ( R o / R o fin~a)decreases exponentially with time which suggests that a correlation of the type R 0 = R 0 fi.~a(1 - e - ' t ) would describe the development of the neck radius not worse than Frenkel's constant power law. A more refined examination is found by observing the behaviors of the neck velocity. The three stages are again clearly evident. The short initial transitory stage is followed by an in101
I
I
I
I
I
I
termediate pseudo-Frenkel stage which exists for only a small fraction of the time required for complete coalescence. During most of the time, however, the velocity shows a decay characteristic to processes involving first-order viscous dissipation. It is useful to compare the simulation resuits with the available experimental measurements of various dimensions. Such data l
I I
I
I
I
I
Cl- Initial I I
I
I
I
I
Ro -6-
1°°
I~ 10-1 100
I
I
I
I
I
I
I II
10 ~
1800C radius
I
t
t
~
102 t , (rnin)
FIG. 6. A comparison between experimental measurements and simulation results for neck size evolution of polystyrene doublets at 180°C (A, O, ×, +, represent four independent experimental runs).
Journal of Colloid and Interface Science, Vol. 9 5, No, 2, October 19 8 3
470
HIRAM AND NIR I
O
i
~-~..~
I
i
i
~. ~.~.
Ro'(gm)
i
I
I
OO
0
0
FRENKE
ESTIMATE
~'-~
~
o ^ ~
""-. "%%,~ %
20C
\\
~0
~\\ ~0 \\\\
100 PS 180-c
0
J 700
I
I
~ 800
r
900
I
I 1000
\\ ~
\
I
1100
L, (l.Lm)
I~G. 7. A comparison between experimental measurements and simulation results for the dependence of neck size on the total length of a doublet, L, (Polystyrene doublets at 180°C.)
are available for the sintering process of polymer particles. Figure 6 depicts the isothermal neck evolution of four different experiments of coalescence of two polystyrene spheres (10). The solid line is obtained by our simulation 0, ~ ~ ) with physical properties chosen to dimensionalize the time scale. A further comparison, independent of the choice of time scale, can be made by following the shrinkage of the longest dimension, L, of the doublet as the neck grows. This evolution is shown in Fig. 7. The experimental data are again from Rosenzweig (10) for polystyrene particles at 180°C. The dotted line is a consequence of Frenkel's m = 1/2 power law. The simulation prediction lies slightly lower than the corresponding measurements possibly due to deviations from pure Newtonian behaviour. The good agreement shown in Figs. 6 and 7 indicates that the polymer sintering process can indeed be simulated by balancing surface and viscous forces and supports the similar simulation to the other coalescence processes at low Reynolds numbers.
Journal of Colloid and Interface Science, Vol. 95, No. 2, October 1983
REFERENCES 1. Gal-Or, B., Klinzing, G. E., and Tavlarides, L. L., Ind. Eng. Chem. Fund. 61, 21 (1969). 2. Mason, B. J., "Clouds, Rain and Rainmaking." Cambridge Univ., Cambridge, 1975. 3. Poste, G., and Allison, A. C., Biochim. Biophys. Acta 300, 421 (1973). 4. Duncan, R., and Pratten, M. K., J. Theor. Biol. 66, 727 (1977). 5. Tadmor, Z., and Gogos, C. G., "Principles of Polymer Processing." Wiley, New York, 1979. 6. Frenkel, J., J. Phys. (USSR) 9, 385 (1945). 7. Kuczynski, G. C., Neuville, B., and Toner, H. P., J. Appl. Polym. Sci. 14, 2069 (1970). 8. Rosenzweig, N., and Narkis, M., Polym. Eng. Sci. 21, 1167 (1981). 9. Ross, J. W., Miller, W. A., and Weatherly, G. C., J. Appl. Phys. 52, 3884 (1981). 10. Rosenzweig, N., D.Sc. dissertation, Technion, Haifa, 1980. 11. Rallison, J. A., and Acrivos, A., J. FluidMech. 89, 191 (1978). 12. Ladyzhenskaya, O. A., "The Mathematical Theory of Viscous Incompressible Flow." Gordon Breach, New York, 1969. 13. Youngren, G. K., and Acrivos, A., J. Fluid Mech. 69, 377 (1975). 14. Reinsch, C. H., Numerische Math. 10, 177 (1967).
Many physical processes involve coalescence of dispersed elements in a continuous medium to form bigger particles. Bubbles and droplets coalescence (1, 2), merger of membrane-coated subcellular vesicles in biological systems (3, 4), and sintering of amorphous polymer particles (5) are but a few examples of the phenomenon. The coalescence dynamics is thus an important step in a variety of processes such as multiphase chemical reactions, sedimentation, emulsification, liquid-liquid extraction, combustion, spray drying fermentation, biological mass transfer and uptake, emulsion polymerization, and polymer processing. These differ in the rheological properties of the fluids involved, the scale of the particles participating, the interfacial properties of the phases, and the time scale during which a complete coalescence is obtained. Nevertheless, in all these cases the driving force for coalescence is the excess in interfacial energy Present address: Department of Chemical Engineering, Imperial College, London SW7 2AZ, England. 2 To whom correspondence should be directed.
which exists once fusion between particles commences. Motion is induced by the tendency of the interface to reduce its area in order to minimize the interfacial energy. Real fluids and real interfaces have complex rheological properties. The interface may exhibit viscous, viscoelastic, or even elastic behavior with both spatial as well as temporal variations, depending on the surface chemical composition and the interaction with the bulk. Yet a simple flow model can serve to unify the approach to the various cases and study the important parameters involved. We consider here processes in which the flow can be assumed to be a motion of Newtonian fluids where the interface is characterized by interfacial tension, a. Furthermore, since the typical velocity is of 0(a/ (~ti + #0)), where ~i and ~0 are the viscosity of the fluid particles and the continuous phase, respectively, it is assumed that the typical Reynolds number based on this velocity and the particle dimension is small enough such that inertia effects can be neglected. Indeed this assumption is valid for the case of polymer sintering processes due to the high viscosities involved and for the case of bio-
462 0021-9797/83 $3.00 Copyright© 1983by AcademicPress,Inc. All rightsof reproductionin any formreserved.
Journalof Colloidand InterfaceScience, Vol,95, No. 2, October1983
SIMULATION
OF COALESCENCE
463
solution of the flow problem. Such an approach is not simple since one has to solve the equations of motion simultaneous with the evolution of the curved interface during the entire merger of the two droplets to the ultimate single spherical particle. In principle two intersecting spherical droplets can be described by a conformal mapping using bipolar coordinates. Unfortunately, this technique involves complex eigenfunctions and the general axisymmetric Stokes streamfunction is yet to be produced. Furthermore, such a mapping involves geometrical assumptions at the intersection contour similar to those introduced by Frenkel (6). Ross et al. (9) have attempted a computer simulation for the sintering of a continuous row of particles using a finite element method for the solution of the equation of motion. Although they have encountered accuracy problems in the initial stages, their finding is most important, where Ro is the radial dimension of the neck namely, Ro does not change with t through between the particles and t is the time. A is a single power law but rather the power varies a constant depending on the physical prop- monotonically during the neck growth. Roerties. Frenkel's result was not based on a senzweig (10) has attempted a finite element solution of the equations of motion, how- solution to the coalescence problem of two ever, it was widely used as a basis for studies particles. He suggested an overall power law involving sintering of spherical amorphous correlation similar to Frenkel's equation, polymer particles. Furthermore, the result even though his results indicated a monowas taken to apply even outside its originally tonic change of power over the limited region assumed region of validity, e.g., beyond the of contact neck radius which he considered. initial stage of coalescence. Two basic deviaIn this communication we solve the equations from Frenkel's (6) model were ob- tion of motion numerically using the integral served. Kuczynski et al. (7) found that the equation representation for the surface vedynamics of coalescence varies with temper- locities. This approach enables a solution for ature where both A and the power of t are the evolution of the interface during the conot constant. They have attributed this to the alescence or sintering process without the neglected non-Newtonian viscoelastic prop- necessity to evaluate the details of the flow erties of the polymeric material. Yet, their field within or out of the particles. We show evidence still supported the fact that the neck that the Frenkel model is generally not apsize grows as a power of time. Similar find- plicable to the entire process of neck growth ings were recently reported by Rosenzweig even in the case of Newtonian fluids with a and Narkis (8) who also obtained small os- constant surface tension. Furthermore, cillatory deviations from the logarithmic de- boundary layer considerations indicate that pendence (see Figs. 1 and 5 in this reference). in the initial transitional stages of coalesOf course, a direct way to check the model cence, deviations from Frenkel's model are and the extent of its validity is an explicit always possible. The general equations are
logical particle fusion in view of their small dimension. It is probably only roughly acceptable for atmospheric droplets merger and grossly violated in the case of bubbles coalescence. The sintering of polymer particles is an example which has been extensively studied experimentally and theoretically (6-9). An approximation based on a surface tension driven flow model for this process was proposed by Frenkel (6). He assumed that the interracial shape of two spherical droplets, during initial stages of coalescence, is described by the intersecting joint of two spheres. Approximating the energy dissipation from geometrical consideration and balancing it with similar estimate for the surface energy change Frenkel obtained the power law R o = A t 1/2 [1]
Journal of Colloid and Interface Science, Vol. 95, No. 2, October 1983
464
HIRAM AND NIR
given in the following section. We then describe the numerical procedure, and the resuits are given and discussed. INTEGRAL EQUATION FORMULATION
0/'/i
= 0
[5]
where n,- is a unit normal pointing outward and cr is a constant surface tension. The velocity at a point x on the surface S can be related to the velocities and forces at all points y at the interface via an integral equation (see Rallison and Acrivos (11), Ladyzhenskaya (12))
. [8]
KijkujnkdSy
1 -
2 fsr &~(V. n)nflSy.
[9]
The only parameter in Eq. [9], X, is the key which enables simulations of different physical processes as polymer sintering (k ---, oo), biological particles fusion (X ~ 0(1)), and bubble coalescence (X ---, 0). The approach to a solution of Eq. [9] is a numerical evaluation of the surface velocities with a quasistationary interfacial shape and a subsequent evolution of the surface. When the problem has an axis of symmetry it is possible to eliminate the angular dependence in Eq. [9]. The velocity and surface unit normal vector have only axial and radial components, uz, ur, nz, and nr, respectively, and [9] is reduced to the form Uz(Z)+ 2 11 -_X) ~ f t-t (Kzzuz(~) + Kzrur(~))d~
(1 + k)ui(x) + (1 - X)
2l f f l (Jzznz(() + J~nr(())V" n(~)d(
)< ~ s, giik(X -- y)uj(y)nk(y)dS,
1 fs, J0(x - y)A~rjk(y)nk(y)dSy 87r/~0
yl 3
bli(X ) At- 2 ~
On.._~j
_
Ix -
[31
where /x indicates a difference operation at the interface (inner minus outer), and a balance of stresses
1
(Xi -- Yi)(Xi -- Yi)
t~ii
Jij(x - y) - Ix ~ yl +
[21
within the droplets, and with a similar expression where tZo replaces ~i in the outer fluid. u, a, and p denote the velocity, stress, and pressure fields, respectively. The interfacial conditions are continuity of velocities, AU/= 0 [4]
Atrijn i = --tr OXi ni
[71
Equation [6] can be rendered nondimensional by scaling u with U = a/2a-~0(1 + X) and x by a. Then, using condition [5] for the difference in stress distribution we have
where
[Oui Ouj] alj = -P~ij + ~ i ~ x i + ~x,]
3 (X i -- yi)(xj - Yfl(Xk -- Yk) Ix - yl5 47r and
Consider two Newtonian droplets, each of radius a, coalescing in a continuous fluid. The fluids obey the Stokes equations Oail = 0 ;
Kijk(X -- y)
1 -Xf t u~(z) + 2 1 + k d-i (Krrur(~) -t- Krzuz(~))d~ [61
where X = ~i/#0. Here the kernels K and d derive from the fundamental solutions to the Stokes equations Journal of Colloid and InterfaceScience, Vol. 95, No. 2, October 1983
_
--
1 ff
"2 -t (J~rnr(~) + Jrznz(())V" n(~)d~ [101
where L = 2l is the length of the merging particles along the axis and
SIMULATION
V'n=
l, d z 2 l k
OF COALESCENCE
+ \-~z] ]
1 R
1+
\dz] ]
[111
with the surface given by r = R(z). The axisymmetric components of K and d are given by Youngren and Acrivos (13). NUMERICAL
SOLUTION
The integral Eq. [9] was solved using the following procedure. The initial shape was assumed to be two spherical droplets intersecting each other at a small joint neck. Surface points were distributed along the interface. The spacing between points was chosen such that a denser distribution was applied where derivatives and curvature were large. The equation was reduced to a set of simultaneous linear algebraic equations for the components of the velocities at the surface points. These velocities were used, in turn, to evaluate the normal displacement of the surface segments during a small interval of time to form a new interfacial shape. The procedure was applied repeatedly from the initial shape until all surface motion diminished. We have used an interactive language (APL) which is very suitable for operations on systems of simultaneous algebraic equations and is even more attractive for use in solutions of evolutionary processes since checks and changes can be made and intermediate surface shapes can be plotted as the solution progresses. The number of surface points, N, on one droplet was chosen as a compromise between the will to increase numerical accuracy and the tendency to reduce the size of matrices to be inverted. The sequence of calculation was as follows. First and second surface derivatives, the components of the unit normal, and the surface curvature at each point were evaluated using a straightforward second-order finitedifference scheme for unequal intervals. For
465
a given point x those were used to calculate the tensors K and d at the surface points y. It should be noted that d has a logarithmic singularity as y ~ x except when x is the point on the axis of symmetry where d is of 0(IX -- y[-l/2) The integral involving d was calculated using a trapezoidal procedure. The surface segment with the singularity was excluded while its contribution was evaluated analytically and added. Repeating this procedure for each point x on the two droplets, a set of 2 x ( 2 N - 1) simultaneous equations for the axial and radial velocities components at all 2N - 1 surface points was produced. Using the fore-and-aft symmetry reduced the rank of the problem to 2 × N. This set was readily inverted by the APL "domino" operation to yield the desired surface velocity. Next the surface evolution was evaluated by considering the motion of each surface point normal to the surface segment during a small time interval dt, i.e., the position increment was dx = (u. n)ndt. The fluid points move of course also tangentially, hence, after the increment dt a new set of fluid points is at hand. In order that the mathematical surface points remain approximately in an appropriate spatial separation after each step, especially in the narrow neck region, they were relocated using surface smoothing by spline functions (14). The magnitude of the time step was chosen small enough to ensure stability of the evolution scheme. At <~AZmi,/ Un max proved stable enough in all cases investigated as is also the case in hyperbolic evolutionary problems. Accuracy was measured by requiring that the volume of the droplets does not change. This is equivalent to satisfying the equation of continuity. In most cases investigated the total volume change during the entire process was a few percent (see Fig. 1). When X ~ 0, however, the changes in volume become large indicating, as was also pointed out by Rallison and Acrivos (11), that multiples of homogeneous eigensolutions of Eq. [9] appear in the scheme. Indeed the cumulative Journalof Colloidand InterfaceScience,Vol. 95, No. 2, October 1983
466
H I R A M A N D NIR
0.00
-
0.04
Error ( -AV -) V -0.08
-0.12
-0,16
_0] I
-0.
1
0
I
-0.6
I
-0.2
I
0.2
0.6
1.0
1-X_X I*X FIG. 1. Total relative change in volume V during the simulation process.
volume change for the case ), = 0.1 was as high as 22%, while for the case ), ~ oo it was an order of magnitude lower. RESULTS A N D DISCUSSION
The evolution of the interface shows a similar behavior for all viscosity ratios. Figure
'In,1 .
. . . . . . . .
2 depicts the dependence of the neck radius, R0, and its growth velocity., Ugo, on time for X ~ ~ and ), = 1. It shows three regions of evolution. An initial transitory period, an intermediate stage and a final decay of growth. For all X investigated, 0.1 < X < ~ , this qualitative dynamic similarity remains. It suggests that a wide spectrum of physical processes differs only slightly when they are driven by surface tension and retarded by viscous dissipation. Indeed, an indication that X is not a parameter with singular behavior is already given in Eq. [9] where the actual parameter appears to be (1 - ),)/(1 + X) and varies in the interval (-1, 1) when X changes from zero to infinity. A further evidence that the evolutionary characteristics of the process are similar at all ), is given in Fig. 3 which shows the inward velocity of the surface on the axis of symmetry, i.e., half the shrinkage speed of the coalescing doublet. At the initial steps the numerical scheme is inadequate to describe the evolutionary process. In the vicinity of the neck there exists a boundary layer where extremely high surface curvature and velocities are involved. This stage is short and its effect on the total 101
I
i
~ , i,,,,
I
Transitory
'"Fre n k ¢ l " ~
t i 1,1,,
t
t
' "Frenkel " "h..
stQ~ 10 o
i , ~1,1~
"
stage
-Transitory ' ~
stage
i
URo
~
5
10o
~
X=I i
10"1 10 -2
t
I
I
, , ,,,I
'
'
' ' ' ''
10 q
I(3' 10 0
I
10-3
I
I
I I I Ill
I
I
I
10.2
I I Illl
I
10-1
I
I
I I I I~
10 0
t a/U
t a/U
FIG. 2. Dimensionless evolution of the neck size, Ro, and its radial velocity, URo, as function of time. (a) X ~ 0% (b) X = 1 (a, initial droplet radius; U = a/2~rtz0(1 + X)).
Journal of Colloid andlnterface Science, Vol. 95, No. 2, October
1983
SIMULATION
0.0
dicate that, even at this early stage, deviations from Frenkel's 1/2 power law [1] may occur. In the inner region the variations in the axial, z, and radial, r, directions are of 0(e) where ~ is a small measure of the gap between the droplets (e.g., the gap before coalescence started or the neck size itself). Since the radial and axial components of the velocity are of different orders of magnitude the continuity equation reduces to
-0.5 u& U_I.(
Y -2£
i
O0
0.15
i
0.30
10 r ar (rUr) = 0. i
t a/U
467
OF COALESCENCE
0.45
0.60
FIG. 3. T h e d i m e n s i o n l e s s c h a n g e o f the v e l o c i t y o f the surface p o i n t o n t h e axis o f s y m m e t r y as f u n c t i o n w i t h t i m e a n d ), ( U = a/27r/~0(1 + X); a, i n i t i a l d r o p l e t radius).
dynamics is transitory. The inaccuracy of the numerical results in this region reflects the fact that for any refinement of the spacing of surface points there is always a short period during which the boundary layer is of a smaller size. Some useful information, however, can be Obtained from a boundary layer consideration. At this stage the surface can be separated into three distinct regions. An outer region where the droplets remain spherical having almost the original size, a narrow inner region including the connecting neck where the curvature and the normal velocity attain extremal values, and an intermediate region separating the two in which the surface gradients are high but the curvature passes from negative to positive values. The exact features of the inner region grossly depend on the relative position of the droplets when coalescence started. It is therefore of no general interest to obtain exact perturbation solution to this boundary layer. We shall only point out the characteristics of the inner region which are of interest in as much as they in-
[12]
Hence the flow outside the droplets, to a first approximation, is a source flow C1 Ur-- - r
[13]
where C1 is a constant source strength. At the neck we have of course URoR0 = C~ which apparently agrees with Frenkel's hypothesis. However, it will become clear shortly that C~ itself depends on R0 as well. From the r component of the Stokes equations it is clear that the pressure in the gap near the neck must be constant. It follows that the normal stress component outside the interface is dominated by viscous terms and is of the dimensionless form O'rr ~--- --2 --Cl /.2 "
[14]
Neglecting the contribution from the interior, and with a surface given by r = R(z), the normal stress balance is approximately of the form 0"rrnrn r =
-
-2
CI R2(I + R '2)
R" 1 + (1 + R'z) 3/2 R(1 + R'a) 1/2
[15] where R' denotes dR/dz. Equation [5] can be integrated once to become Journal of Colloid and Interface Science, Vol. 95, No. 2, October 1983
468
HIRAM AND NIR
(1
+ R'2) 1/2 = - C 1
+ C2R2
[16]
where C2 is a constant. Hence at the neck where z = R' = 0 we find C1 = C2R 2 - 1
[17]
where C2 must be determined by integrating Eq. [ 16] and matching with the outer regions. It is sufficient here to indicate that in view of Eq. [17] deviation from URoRo = constant are evident even at this early stage and a Frenkel 1/2 power law may be violated. Indeed, the boundary layer analysis applies only to a very limited period of time. After this transitory period it is expected that the interface shape will not depend on the specific features of the early stages, and will further evolve in a quasi-steady manner to complete coalescence. It is at this stage that application of a numerical scheme to solve Eq. [9] for the entire shape is appropriate. It is also at this stage that experimental results concerning sintering of polymeric particles are available for comparison (10). In the intermediate stage only a small change of surface area occurs, however, it is at this period where most of the neck growth is evident since dissipation is also small. Figure 2 shows that the neck evolution process does not obey Frenkel's constant power law for an extended period of time. This result was also suggested by Ross et al. (9) who applied a finite element technique to study the evolution of necks in a continuous array of coalescing spheres. In order to better illustrate the above conclusion we attempt to apply Frenkel's power m to a time interval within the intermediate region. If R o a t m then the decline in the velocity should obey URoOttm-1. F o r the limited interval of time it is possible to find m depending on the specific h. These are shown in Fig. 4. It is clear, however, from our results (Fig. 2) that the velocity can be only roughly correlated by a constant power m. This power depends strongly on the choice of time interval in which the correlation is considJournalof ColloidandInterfaceScience,Vol. 95, No. 2, October 1983
ered. Rosenzweig (10) obtained results for a finite element simulation of the sintering of two Newtonian fluid spheres (k --~ oo). He studied the period during which the neck varies over the interval 0.6 ~< R0 ~< 0.8 for which he obtained a power m = 0.3, markedly below the m = 0.5 Frenkel value as well as below our results for X ~ ~ . A careful differentiation of Rosenzweig's plot shows no existence of a constant power growth. In addition, the low value obtained for m merely reflects the fact that his results were obtained in a time interval almost beyond the pseudo-Frenkel stage and support our finding that m is decreasing with evolution. Beyond the intermediate stage comes the final period of coalescence where most of the changes in surface area occur accompanied by strong viscous dissipation and damping of the evolution. Variations in surface curvature are mild and the particle assumes its ultimate spherical shape following a slow first order exponential decay of the motion. Figure 5 depicts a semilogarithmic plot of the decay (1 - R o / R o nn~) and of the neck velocity for the entire process for two different X's. For all practical purposes it is clear that
0.6 O
0.5 Frenket's Power m
0.4
0.3
°:o
I
5
I
I
I
3o
oo
FIG. 4. T h e dependence of Frenkel's power, m, on the viscosity ratio, ~,.
469
SIMULATION OF COALESCENCE 101
101 k=30
3
~
a-finradius aL
k=0.3
3~'2a-radius final
10c
10 0
"
~ffl
0
~
~
1-, Ro
.
I
I
1
0.1
0.2
0.3
t a/U
161
0.4
c~.1
0
o'.2
0'.3
t a/U
o'.4
d.5
0.6
FIG. 5. Dimensionless evolution of(1 - Ro/Ro~n~) and uR0 as a function of time. (a) X = 30, (b) X = 0.3 (U = a / 2 r ~ ( 1 + X); a, initial droplet radius).
1 - ( R o / R o fin~a)decreases exponentially with time which suggests that a correlation of the type R 0 = R 0 fi.~a(1 - e - ' t ) would describe the development of the neck radius not worse than Frenkel's constant power law. A more refined examination is found by observing the behaviors of the neck velocity. The three stages are again clearly evident. The short initial transitory stage is followed by an in101
I
I
I
I
I
I
termediate pseudo-Frenkel stage which exists for only a small fraction of the time required for complete coalescence. During most of the time, however, the velocity shows a decay characteristic to processes involving first-order viscous dissipation. It is useful to compare the simulation resuits with the available experimental measurements of various dimensions. Such data l
I I
I
I
I
I
Cl- Initial I I
I
I
I
I
Ro -6-
1°°
I~ 10-1 100
I
I
I
I
I
I
I II
10 ~
1800C radius
I
t
t
~
102 t , (rnin)
FIG. 6. A comparison between experimental measurements and simulation results for neck size evolution of polystyrene doublets at 180°C (A, O, ×, +, represent four independent experimental runs).
Journal of Colloid and Interface Science, Vol. 9 5, No, 2, October 19 8 3
470
HIRAM AND NIR I
O
i
~-~..~
I
i
i
~. ~.~.
Ro'(gm)
i
I
I
OO
0
0
FRENKE
ESTIMATE
~'-~
~
o ^ ~
""-. "%%,~ %
20C
\\
~0
~\\ ~0 \\\\
100 PS 180-c
0
J 700
I
I
~ 800
r
900
I
I 1000
\\ ~
\
I
1100
L, (l.Lm)
I~G. 7. A comparison between experimental measurements and simulation results for the dependence of neck size on the total length of a doublet, L, (Polystyrene doublets at 180°C.)
are available for the sintering process of polymer particles. Figure 6 depicts the isothermal neck evolution of four different experiments of coalescence of two polystyrene spheres (10). The solid line is obtained by our simulation 0, ~ ~ ) with physical properties chosen to dimensionalize the time scale. A further comparison, independent of the choice of time scale, can be made by following the shrinkage of the longest dimension, L, of the doublet as the neck grows. This evolution is shown in Fig. 7. The experimental data are again from Rosenzweig (10) for polystyrene particles at 180°C. The dotted line is a consequence of Frenkel's m = 1/2 power law. The simulation prediction lies slightly lower than the corresponding measurements possibly due to deviations from pure Newtonian behaviour. The good agreement shown in Figs. 6 and 7 indicates that the polymer sintering process can indeed be simulated by balancing surface and viscous forces and supports the similar simulation to the other coalescence processes at low Reynolds numbers.
Journal of Colloid and Interface Science, Vol. 95, No. 2, October 1983
REFERENCES 1. Gal-Or, B., Klinzing, G. E., and Tavlarides, L. L., Ind. Eng. Chem. Fund. 61, 21 (1969). 2. Mason, B. J., "Clouds, Rain and Rainmaking." Cambridge Univ., Cambridge, 1975. 3. Poste, G., and Allison, A. C., Biochim. Biophys. Acta 300, 421 (1973). 4. Duncan, R., and Pratten, M. K., J. Theor. Biol. 66, 727 (1977). 5. Tadmor, Z., and Gogos, C. G., "Principles of Polymer Processing." Wiley, New York, 1979. 6. Frenkel, J., J. Phys. (USSR) 9, 385 (1945). 7. Kuczynski, G. C., Neuville, B., and Toner, H. P., J. Appl. Polym. Sci. 14, 2069 (1970). 8. Rosenzweig, N., and Narkis, M., Polym. Eng. Sci. 21, 1167 (1981). 9. Ross, J. W., Miller, W. A., and Weatherly, G. C., J. Appl. Phys. 52, 3884 (1981). 10. Rosenzweig, N., D.Sc. dissertation, Technion, Haifa, 1980. 11. Rallison, J. A., and Acrivos, A., J. FluidMech. 89, 191 (1978). 12. Ladyzhenskaya, O. A., "The Mathematical Theory of Viscous Incompressible Flow." Gordon Breach, New York, 1969. 13. Youngren, G. K., and Acrivos, A., J. Fluid Mech. 69, 377 (1975). 14. Reinsch, C. H., Numerische Math. 10, 177 (1967).