- Email: [email protected]
Robust Digital Image Analysis of Pendant Drop Shapes
JOURNAL OF COLLOID AND INTERFACE SCIENCE ARTICLE NO.
177, 658–665 (1996)
0080
Robust Digital Image Analysis of Pendant Drop Shapes DAVID B. THIESSEN, DOMINIC J. CHIONE, CLINTON B. MCCREARY,
AND
WILLIAM B. KRANTZ 1
Department of Chemical Engineering, University of Colorado, Boulder, Colorado 80309-0424 Received May 10, 1994; accepted June 23, 1995
The advent of digital image analysis permits the use of several hundred data points defining the profile of a pendant drop in order to determine the surface tension from a fit to a solution to the Young–Laplace equation. This paper discusses several improvements to both the experimental apparatus for carrying out pendant drop measurements and to the methodology for analyzing the drop profiles in order to determine the surface tension. Several standard minimization algorithms were evaluated for fitting the theoretical and experimental drop profiles. These algorithms are very reliable, although there is a tradeoff between robustness and speed of computation. The most robust algorithm was found to converge irrespective of the ‘‘roughness’’ of the experimental drop profile or the accuracy of the initial estimates of the fitting parameters. Whereas the standard deviation of the surface-tension data provides a measure of the precision, the objective function is shown to be an indicator of systematic error in the data. Surface-tension measurements are reported for two standard liquids, methanol and toluene. A standard deviation of 0.016 dyn/cm is obtained for methanol using the full drop profile; this is significantly below the value of 0.036 dyn/cm obtained using the selected plane method, thus indicating that this improved pendant drop technique can give very precise measurements of surface tension. q 1996 Academic Press, Inc.
1. INTRODUCTION
The analysis of the shape of a pendant liquid drop is the basis for a very accurate method of measuring surface tension (1). Moreover, it offers the advantages of requiring only very small samples and being amenable to measuring dynamic interfacial tensions. This method has been applied to common liquids under ambient conditions as well as to high-temperature molten metals (2) and high-pressure liquid/vapor systems near the critical point (3). In recent years a number of groups have developed techniques to automate the pendant drop method using digital image analysis (4– 11). Rapid analysis of surface tension by computer has been applied to tracking the dynamically changing surface tension in a nonequilibrium system undergoing mass transfer (7). Computer analysis of surface tension by the pendant drop 1
To whom correspondence should be addressed.
0021-9797/96 $12.00 Copyright q 1996 by Academic Press, Inc. All rights of reproduction in any form reserved.
/ m4585$3915
12-22-95 21:57:10
method involves first locating the edge of the drop from a digital image which yields a set of points spaced along the drop profile. Next, the Young–Laplace equation for an axisymmetric interface is solved for an initial estimate of the parameters on which it depends. The parameters which govern the theoretical drop profile are adjusted and the equation is solved repeatedly until a best fit to the experimental points is obtained. The best fit is determined by minimizing an objective function which is typically a summation of the squared distances between the experimental points and the theoretical drop profile. Several problems are still encountered when using the pendant drop method to determine surface tension. The apparatus which is commercially available has some shortcomings related to operation at pressures different from ambient, temperature control, and stability of the light source used for illuminating the drop. Some methods used for determining the experimental drop profile yield ‘‘rough’’ shapes which cause problems in fitting the Young–Laplace equation. Moreover, some choices employed for the objective function used in fitting the theoretical and experimental drop profiles are not optimal. In addition, the Newton–Raphson numerical technique used to determine the fitting parameters is prone to nonconvergence for ‘‘rough’’ drop profiles or when the initial estimates of the parameters are not sufficiently accurate. Finally, although the standard deviation provides a measure of the precision of the data, no attempts have been made to assess any systematic error which might be present in the fit between the theoretical and experimental drop profiles. This paper addresses all of the aforementioned problems in an attempt to improve the precision of the pendant drop method for measuring surface tension. Moreover, this presents the first implementation of software to analyze pendant drop data using the Apple/Macintosh operating system. This paper is organized as follows. First, a brief review of prior studies will be given along with a critique of these studies which elaborates on some of the above mentioned problems. This is followed by a summary of the experimental apparatus which incorporates improvements which have been made to the commercially available instrument. The experimental procedure and data analysis then are summa-
658
coida
AP: Colloid
DIGITAL IMAGE ANALYSIS OF PENDANT DROPS
rized; the latter incorporates several improvements which have been made in the edge-detection method, calibration, and fitting algorithm; in particular, several standard minimization algorithms for determining the fitting parameters (which include the surface tension) are evaluated. The most robust algorithm was found to converge irrespective of the roughness of the drop profile or the accuracy of the initial guesses for the fitting parameters. 2. PRIOR STUDIES
659
to the fitting parameters. The method also prevented large oscillations in the parameters as the minimum was approached. Lin et al. (7) used a system having a resolution of 480 1 512 pixels with 256 levels of gray. These investigators fully utilized this high resolution in their edge detection routine and employed a local threshold and interpolation technique to provide smoother drop profiles; the local threshold was chosen to be halfway between the local maximum and minimum light intensity near the drop interface. However, they still encountered convergence problems owing to their use of a Newton–Raphson method to determine the four parameters involved in fitting the Young–Laplace equation to the experimental drop profile. The recent work of Hansen and Rodsrud (10) employed a system with a resolution of 256 1 256 pixels with 64 levels of gray. Their edge detection routine was similar to that used by Lin et al. (7). However, the objective function in their minimization routine used the square of the horizontal rather than the minimum distance between the theoretical and experimental drop profiles; as such, this method incurred more error near the apex of the drop where the horizontal distance can be quite significant. Nonetheless, their standard deviations were quite small, 0.01 to 0.06 dyn/ cm although some of their surface tensions differed significantly from currently accepted values; this could indicate some type of systematic error. Hansen (11) has improved on the method of Hansen and Rodsrud by utilizing the sum of squared normal distances as an objective function. This brief review indicates that image digitization can significantly increase the accuracy of the pendant drop technique by employing many points to determine the experimental drop shape from which the surface tension can be extracted. Efficient routines for determining the fitting parameters are based on least-squares minimization employing an objective function based on the square of the minimum distance between the experimental and predicted drop profile. Determination of the fitting parameters which include the surface tension via Newton–Raphson methods is complicated by convergence problems owing to ‘‘rough’’ drop profiles or inaccurate estimates of the initial parameter values. Whereas the standard deviation provides a measure of the scatter in the data, systematic errors can result in anomalously high or low values of the surface tension determined by the pendant drop method; clearly it would be of value to develop some measure of any such systematic errors. Although personal computers are now used to implement the image digitization and profile fitting, no software has been developed for the Apple/Macintosh operating system which offers a very powerful user-friendly graphical interface.
Early applications of the pendant drop technique by Andreas et al. (12) employed either the inflection plane or the selected plane method; these methods correlated the surface tension in terms of a drop-shape factor. Subsequent studies extended the tables for the drop-shape factor (13–17). These methods employed relatively few data points taken from the drop shape and hence were less accurate than other methods for determining the surface tension. The advent of image digitization permitted more accurate fitting of the drop shape using mainframe computers (18), thereby allowing determination of the surface tension by fitting the Young–Laplace equation to the full drop profile; however, the digitization was done on enlarged photographs of the drop. The first use of direct digitization of the drop without using any photographs was done by Girault et al. (19). Rapid advances in image digitization permitted resolution as high as 1066 1 575 pixels (20). Rotenberg et al. (4) improved the method of fitting the Young–Laplace equation to the data by employing rigorous least-squares minimization of an objective function equal to the sum of the squares of the minimum distance between the predicted and experimental drop profile. This minimization involves four parameters, the coordinates and radius of curvature of the apex of the drop and the Bond number which contains the unknown surface tension. These parameters are determined by a Newton-Raphson method employing exact derivatives of the objective function and incremental loading. Anastasiadis et al. (5) carried out image digitization and fitting of the Young–Laplace equation for the drop shape using a personal computer. They used a global thresholding method together with a smoothing algorithm for edge detection. They utilized a robust shape comparison method for fitting the Young–Laplace equation to the drop profile. This method is less sensitive to large errors in individual edge points. Optimization was reduced to a single shape parameter and carried out by exhaustive search. The standard deviation of their data, 0.3 dyn/cm, was larger than that typically obtained for the generally less accurate method of selected plane (0.1 dyn/cm). Jennings and Pallas (6) made a significant improvement in 3. EXPERIMENTAL APPARATUS the efficiency of the fitting method by utilizing a modified Newton–Raphson method which did not require evaluation of A commercial pendant drop apparatus (Rame-Hart, Inc.) the second derivatives of the objective function with respect with several modifications was used in this work. The appa-
/ m4585$3915
12-22-95 21:57:10
coida
AP: Colloid
660
THIESSEN ET AL.
ratus as purchased consisted of an optical rail on which was mounted a three-axis stage supporting an environmental chamber and heated syringe assembly, a three-axis camera support stage, and a tungsten filament lamp with a variac power supply. Temperature control of the environmental chamber could be accomplished by circulation of water at a controlled temperature through a jacket or by resistive heating with feedback control. Temperature control of the syringe could be accomplished only by resistive heating regulated by a separate controller. This resulted in temperature differences between the environmental chamber and the syringe. Consequently the syringe housing was modified by replacing the resistive heater with a water jacket and circulating water from the same temperature-controlled bath through both the chamber jacket as well as this syringe jacket. This resulted in a negligible temperature difference between the pendant drop and the chamber walls and maintained the temperature constant to within {0.17C. The syringe housing also was modified to accommodate a Hamilton (Model 1001) gas-tight syringe. A gas-tight syringe is necessary when working with reduced chamber pressures and to eliminate cooling effects owing to evaporation of volatile liquids. The final modification of the apparatus involved replacing the lamp variac with a regulated power supply which reduced voltage drift from 0.2 V to less than 0.05 V. The pendant-drop apparatus rests on a bench-top vibration-isolation table (Ehrenreich Photo-Optical Industries Model 78221), which consists of massive steel plates floating on three air bladders. Water hoses and electrical cords which are attached to the apparatus are supported from springs suspended from the lab ceiling to reduce vibrations transmitted to the apparatus. A Panasonic CCD camera with a Nikon lens (micro NIKKOR 55 mm) and a 27.5-mm extension tube are used for imaging the drop. A digital image of a pendant drop is acquired using an image-capture card and associated software (Video Image 1200, Scion Corp.) on a Macintosh IIfx personal computer. 4. EXPERIMENTAL PROCEDURE
First, the air table and the optical stand must be adjusted to ensure that a pendant drop will hang exactly perpendicular with respect to the camera and the light source. Analysis via the technique developed here then requires that a column of pixels in the image be parallel to the gravity vector. The apparatus is aligned using a plumb bob consisting of a weight hanging from a fine copper wire attached to the apparatus. The weight is submerged in a beaker of water to damp oscillations. The plumb bob wire is imaged with the CCD camera. The level of the apparatus is adjusted until the edge of the wire in the image aligns with a column of pixels. This leveling procedure must be carried out periodically because
/ m4585$3915
12-22-95 21:57:10
coida
of changes in the volume of the air bladders that support the vibration-isolation table. Before surface-tension measurements are made, a calibration image of a 3.0 { 0.0025 mm precision ball bearing (Industrial Tectonics) is captured. This calibration yields the dimensions of a pixel. The horizontal dimension of a pixel, lx , is not necessarily the same as its vertical dimension, lz , owing to an aspect-ratio error which could be caused by the camera and/or the image-capture system. This calibration involves suspending the ball bearing at the same location as the pendant drop; the camera adjustments and light-source settings are not changed after the calibration image is made. Surface-tension measurements were made on HPLC-grade methanol and toluene without further purification. The syringe was cleaned and dried and then rinsed several times with the liquid to be measured. A pool is made on the bottom of the environmental chamber with liquid from the final rinse. The syringe then is filled and placed in the syringe housing. Water from a constant temperature bath then is circulated through the chamber and syringe housing jackets; the liquid and air in the chamber are allowed to equilibrate for 1–2 h. The circulating bath is turned off for 60 s before a pendant drop image is captured to allow vibrations to damp out. The room lights are then turned off and the image is captured by the computer directly from the CCD camera; this yields considerably improved image quality over capturing a videotaped image. After a drop has been imaged, it is allowed to fall to the bottom of the chamber and a new drop is formed. The bath is turned on for 60 s and then the image capture sequence is repeated. Typically ten different drops are imaged in a single run. The digital images of the calibration sphere and pendant drops are stored as graphics files on the computer hard drive until they are analyzed. 5. DATA ANALYSIS
A typical experimental run yields one graphics file containing an image of the calibration sphere and a number of graphics files containing pendant drop images. The images are in 560 1 420 pixel format with 8 bits per pixel gray scale resolution (256 shades of gray). Image-analysis software was developed on the Macintosh IIfx computer in the C programming language. Picture files stored in the PICT, TIFF, or Encapsulated PostScript (EPSF) formats can be imported by the program. The precision sphere image is analyzed to obtain calibration parameters which are used in the subsequent analysis of the pendant drop images for surface-tension calculation. A subroutine was written to open the graphics file containing the image data, read the size of the image, read the pixel intensities, and store these values in a large two-dimensional array. The edge of an object defined by a digital image is found by searching along rows or columns of pixels looking for an abrupt intensity change. The sides of the sphere or drop
AP: Colloid
661
DIGITAL IMAGE ANALYSIS OF PENDANT DROPS
are found by searching along rows while the bottom is found by searching along columns. An intensity difference between adjacent pixels greater than some specified value indicates the vicinity of the interface. Once the approximate position of the interface has been located, a local threshold intensity is determined, which is taken to define the ‘‘true’’ edge. The local threshold is defined to be the average of 20 pixels in the same row or column, 10 on each side of the interface, starting 5 pixels away from the approximate edge position. The fractional coordinate defining the edge of the object along the search direction is then found by linear interpolation between adjacent pixels which bound the threshold intensity. Successive edge coordinates defining both sides of the axisymmetric object are then stored in an array. This technique locates the interface within less than one pixel with floating point/integer profile coordinates rather than the integer/integer profile coordinates used in earlier less accurate techniques. A subroutine called CALIBRATE analyzes the image of a precision sphere with known radius, Rs , in order to determine the dimensions of a pixel. The horizontal and vertical pixel dimensions, lx and lz , are determined by fitting a circle to the edge coordinates (xi , zi ). An objective function defined as M
Fc Å ∑ {[l 2x (xi 0 x0 ) 2 / l 2z (zi 0 z0 ) 2 ] 1 / 2 0 Rs } 2
[1]
i Å1
is minimized with respect to the four parameters lx , lz , x0 , and z0 , where the latter two are the pixel coordinates of the center of the circle being fit to the data. The parameters are initialized from the analytical solution for a circle which passes through three widely spaced edge points. A minimization subroutine then is called to minimize the objective function given by Eq. [1]; these subroutines will be discussed following a description of how the surface tension is determined. To analyze for the surface tension of a pendant drop, the edge of the drop is located as described above. In order to provide good initial estimates of the parameters, the surface tension is first estimated using the method of selected plane. The equatorial diameter, de , of the pendant drop is found by searching for a maximum diameter starting from the top of the drop. The selected plane diameter, ds , corresponding to the diameter of the drop at a location equal to one equatorial diameter above the drop apex, is then determined. These measurements are corrected for the aspect ratio error determined from the calibration using the precision sphere. The ratio S å ds /de is correlated to a shape factor, 1/H, using Misak’s (16) equation. The surface tension, g, is determined from the equation gÅ
/ m4585$3915
Drgd 2e , H
12-22-95 21:57:10
[2]
coida
where Dr is the density difference between the drop and its surrounding fluid and g is the gravitational acceleration. This initial estimate of the surface tension then is used in a procedure to determine the surface tension by fitting a solution to the Young–Laplace equation to the experimentally determined profile coordinates; approximately 1000 points were used to determine the drop profile. The Young– Laplace equation relates the curvature of an interface to the pressure difference across that interface. In a coordinate system whose origin is at the drop apex, the Young–Laplace equation for an axisymmetric pendant drop can be recast as three coupled first-order differential equations [Rotenberg et al. (4)], df sin f Å 2 0 bz 0 ds x dx Å cos f ds dz Å sin f, ds
[3]
where b å DrgR 20 / g is the Bond number; x and z are the horizontal and vertical coordinates, respectively, nondimensionalized with respect to R0 , the radius of curvature at the drop apex; s is the arc length along the drop surface measured from the drop apex nondimensionalized with respect to R0 ; and f is the angle between the surface tangent and the horizontal. The differential equations described by [3] must satisfy the conditions that f Å x Å z Å 0 at s Å 0. An initial estimate of b is obtained from a correlation between b and the diameter ratio S [Hansen and Rodsrud (10)], the latter being obtained from the method of selected plane. The differential equations given by [3] are integrated by expanding each of the variables f, x, and z in a four-term Taylor series starting at s Å 0 and propagating forward in s; the values of the first, second, and third derivatives in this Taylor series are obtained from Eq. [3], evaluating them at the point of interest. The experimental edge coordinates in terms of pixels, (xi , zi ), are converted to dimensionless form as ji Å
lx (xi 0 x0 ) , R0
zi Å
lz (zi 0 z0 ) , R0
[4]
where lx and lz are the horizontal and vertical pixel dimensions, respectively, determined from the calibration; x0 and z0 now denote the pixel coordinates of the drop apex; and R0 is the radius of curvature at the drop apex. The four parameters b, x0 , z0 , and R0 are varied in order to minimize the difference between the theoretical and experimental drop shapes. An objective function of the form used by Rotenberg et al. (4) is minimized,
AP: Colloid
662
THIESSEN ET AL. M
Fd (au ) Å R 20 ∑ [d(ui , £: au )] 2 ,
[5]
i Å1
where au denotes the parameter vector. The operator d represents the minimum nondimensional distance between an experimental point ui Å ( ji , zi ) and the theoretical curve £ Å [x(s), z(s)] determined from Eq. [3]. The objective function for the pendant drop is minimized using the same general purpose minimization algorithm as is used in the calibration procedure described above. After the objective function has been minimized, the surface tension is calculated from the optimized values of the Bond number and radius of curvature at the drop apex, R0 , according to the following identity: gå
DrgR 20 . b
[6]
The minimum value of the objective function is recorded along with the surface tension in the output from the program. Efficient methods of minimizing objective functions of the form given by Eq. [5] are based on a Taylor series expansion of Fd about a given point in parameter space, au 0 : Fd (au ) É Fd (au 0 ) / (au 0 au 0 )rÇ Fd (au 0 ) /
1 (au 0 au 0 )(au 0 au 0 ): ÇÇ Fd (au 0 ). 2
[7]
Assuming that Fd is well approximated by such a quadratic form near the minimum, and that the derivatives can be computed, an iterative scheme for optimizing the parameter vector au is given by au k /1 Å au k 0 [ÇÇ Fd (au k )] 01rÇ Fd (au k ),
[8]
where [ÇÇ Fd (au k )] 01 is known as the inverse Hessian matrix. This is the well-known Newton–Raphson method, which is very efficient when the objective function is of quadratic form but is susceptible to nonconvergence in cases where it is not. At the opposite extreme are highly robust algorithms which make no assumptions about the form of the objective function but simply march steadily downhill toward the minimum. A number of standard algorithms exist which constitute a compromise between these two extremes. Minimization algorithms from three categories have been evaluated in this work, including algorithms which utilize: (1) function evaluations only, (2) function evaluations and first partial derivatives of Fd with respect to the parameters, (3) function evaluations and first and second partial derivatives of Fd with respect to the parameters.
/ m4585$3915
12-22-95 21:57:10
coida
Computation speed is found to be greatly enhanced by using exact derivative evaluations (6). Although algorithms from the first category are generally slower than those from the latter two categories, they do have significant advantages in terms of flexibility and robustness. In particular, it is sometimes advantageous to use a ‘‘robust’’ fitting function as opposed to the sum of squared residuals given by Eq. [5]. Robust objective functions generally are not as smooth as that given by Eq. [5] and may have discontinuities in the derivatives which preclude the use of minimization algorithms in categories 2 and 3. A modular program has been created in order to provide maximum flexibility in image analysis. The minimization algorithms are standard off-the-shelf routines which are very efficient and reliable. The implementation of these algorithms is that given by Press et al. (21) with only minor modifications. These algorithms are completely generalized so that they can be passed the name of the subroutine which evaluates the objective function along with the number of parameters on which to optimize. An algorithm known as Powell’s method belongs to the first category of algorithms listed above. When Fd is well approximated by a quadratic form then surfaces of constant Fd near the minimum are hyperellipsoids. When searching for the minimum, it is most efficient to search along the directions of the principal axes of the ellipsoid. For a quadratic function of N variables the minimum can be found with only N line minimizations in this manner. The directions of the principal axes thus are called ‘‘noninterfering’’ or ‘‘conjugate’’ directions. In the more general case where Fd is not necessarily of quadratic form, Powell’s method attempts to find noninterfering directions and repeats the cycle of N line minimizations until it converges to a minimum. This algorithm has never failed to converge to a minimum of the objective function in the many cases which were run. A method which uses first derivatives (category 2) is the Broyden–Fletcher–Goldfarb–Shanno (BFGS) variable metric method. This method is a quasi-Newton method in that it uses gradient information to iteratively approximate the inverse Hessian matrix from which new estimates of the optimal parameter vector are computed according to a Newton-type formula (Eq. [8]). The first derivatives of the objective function are given by ÌFd Å2 Ìaj
M
∑ ri i Å1
Ìri Ìaj
j Å 1, . . . , N,
[9]
where ri å R0 d(ui ,£ : au ). The partial derivatives of the individual residuals were evaluated in the manner of Jennings and Pallas (6). The algorithm uses a backtracking scheme to avoid stepping too far in the Newton direction for cases where the function is poorly approximated by a quadratic form.
AP: Colloid
663
DIGITAL IMAGE ANALYSIS OF PENDANT DROPS
The third method which was evaluated is specifically intended for nonlinear least squares minimization. The Levenberg–Marquardt method uses evaluations of the Hessian matrix (category 3) as well as the gradient to perform minimization. The Hessian matrix is approximated as Ì 2 Fd É2 Ìaj Ìak
M
∑ i Å1
Ìri Ìri Ìaj Ìak
j, k Å 1, . . . , N.
[10]
This approximation neglects a term with second derivatives of ri . Inclusion of the second derivative term can actually cause instability in the procedure (21). This algorithm uses a combination of Newton’s method (Eqs. [8–10]) and the steepest descent method. The steepest descent method consists of taking a step in the direction of the negative gradient of Fd . The algorithm uses the steepest descent method far from the minimum and switches smoothly to Newton’s method as the minimum is approached. The computation speed of the three algorithms was compared for a typical pendant drop with 934 edge points. The CPU times (on a Power Macintosh 7100/80 AV) necessary to minimize the objective function were 66 seconds for Powell’s method, 4.2 seconds for the variable metric method, and 0.57 seconds for the Levenberg–Marquardt method. Limited tests indicate that each of these methods converge reliably even in cases where errors as high as 100% in the initial parameter estimates were introduced. Following minimization, it is possible to estimate the uncertainty in the fit parameters from the Hessian matrix evaluated at the minimum of the objective function (6). The confidence region is bounded by the ellipsoid given by (au 0 au m )(au 0 au m ): ÇÇ Fd (au m ) Å DFd ,
[11]
where au m is the optimal parameter vector and DFd is given by DFd Å 2s 2NF(N, M 0 N)
where s 2 is the variance [Fd (au m )/(M 0 N)], and F(N, M 0 N) is the F-test statistic. The application of Eq. [11] to determine the confidence ellipsoid is strictly valid only in the case where the errors are Gaussian distributed. Taking the confidence intervals of the individual parameters to be the extreme values of the ellipsoid, the confidence interval for the surface tension can be determined from Eq. [6]. At the 95% confidence level the confidence interval for the surface tension of a typical water drop was {0.13 dyn/cm. This should be a conservative estimate of the error in these measurements. 6. RESULTS
Surface-tension measurements were made on methanol and toluene at 207C in air saturated with the liquid’s vapor.
/ m4585$3915
12-22-95 21:57:10
coida
TABLE 1 Surface Tension Measurements for Methanol in Air at 207C Surface tension (dyn/cm)
Drop number 1 2 3 4 5 6 7 8 9 10 Mean Standard deviation Percent deviation
Selected plane
Profile fitting
Dimensionless objective function
22.314 22.312 22.243 22.226 22.249 22.222 22.263 22.217 22.231 22.231 22.251 0.036 0.16
22.452 22.444 22.454 22.457 22.457 22.434 22.422 22.456 22.425 22.418 22.439 0.016 0.069
0.0036 0.0035 0.0041 0.0045 0.0045 0.0046 0.0040 0.0044 0.0036 0.0036 0.0041 0.0004 10.5
The density difference between the liquid and gas phases was taken as 0.7904 g/cm3 for methanol and 0.8656 g/cm3 for toluene. The vapor-phase density was calculated from the Ideal Gas Law using Raoult’s Law to estimate the composition of the gas phase. Table 1 gives the measured values of surface tension for ten separate drops of methanol using both the selected plane and profile-fitting methods. The average value of the surface tension determined by the profile-fitting method is 22.44 { 0.02 dyn/cm. This value falls within the uncertainty limits of the accepted literature value of 22.5 { 0.1 dyn/cm which was obtained by the capillary-rise method (22). The selected plane method yields a value almost 1% lower than the profile-fitting method and gives a larger standard deviation [22.25 { 0.04 dyn/cm]. The objective function values reported in Tables 1 and 2 were made dimensionless by dividing Eq. [5] by R 20 . The results for toluene are presented in Table 2. Again, ten separate drops were measured by the two methods mentioned above. The average value of the surface tension determined by the profile-fitting method is 28.20 { 0.07 dyn/cm. This is 1% lower than the accepted literature value of 28.5 { 0.1 dyn/cm obtained by the maximum bubble-pressure method (22). The selected plane method gives a surface-tension value 0.6% lower than the profile-fitting method [28.04 { 0.07 dyn/cm], however the two methods have the same standard deviation for the toluene drops. The standard deviation is considerably larger for the measurements of toluene drops than for the methanol drops. It is proposed here that the minimized value of the objective function be used to indicate the quality of a surfacetension measurement. One would like to know whether a given value of the objective function indicates a mismatch between the theoretical and experimental drop shapes or
AP: Colloid
664
THIESSEN ET AL.
TABLE 2 Surface Tension Measurements for Toluene in Air at 207C Surface tension (dyn/cm)
Drop number 1 2 3 4 5 6 7 8 9 10 Mean Standard deviation Percent deviation
Selected plane
Profile fitting
Dimensionless objective function
27.949 28.011 28.053 28.099 28.149 28.089 27.903 28.082 28.003 28.034 28.037 0.074 0.263
28.105 28.173 28.149 28.132 28.304 28.277 28.146 28.297 28.148 28.230 28.196 0.074 0.263
0.0049 0.0052 0.0055 0.0054 0.0071 0.0063 0.0051 0.0060 0.0052 0.0057 0.0056 0.0007 11.67
whether the value arises strictly from ‘‘noise’’ in the digital imaging system. The noise level can be seen clearly by examining the located edge points from an image of a plumb wire to which the camera has been carefully aligned (Fig. 1). The horizontal axis of the plot in Fig. 1 is greatly expanded in order to show the noise level in the data. The x-coordinates of the edge points have a range of 0.225 pixels and a standard deviation of 0.04 pixels. The calibration procedure described in Section 4 involves imaging a precision sphere, locating the edge, and minimizing the objective function defined by Eq. [1]. The minimized objective function found for a precision sphere can serve as a reference value to compare with the objective function values obtained for drops. By using a precision sphere with a diameter close to the equatorial diameter of a given pendant drop, roughly the same number of image points are found for the sphere and the drop. This allows for a direct comparison of the objective-function values. The magnitude of the objective function for a precision sphere will indicate the noise level of the imaging system. A typical value of the minimized objective function for a calibration sphere (nondimensionalized with respect to R 2s ) is 1.6 1 10 04 . This value is approximately 30 times smaller than the objective function values obtained for drops in this work. The objective-function values for drops in this work are obviously above the noise level, indicating a mismatch between the experimental and theoretical drop shapes. This mismatch could arise from camera misalignment, drop vibration, needle-tip imperfections, or surface-tension gradients on the drop. Surface-tension gradients could arise from nonuniform temperatures or a nonuniform distribution of impurities on the drop surface. One significant source of error in the pendant drop measurements arises from vibrations; although the vibration table greatly helps in isolating the pendant drop apparatus from room vibrations, it is not
/ m4585$3915
12-22-95 21:57:10
coida
massive enough to eliminate all vibrations. The pendant drop apparatus should be mounted on a far more massive vibration isolation table to minimize this source of error. Future work will concentrate on improving the objective-function values by identifying and correcting other causes of error in the measurements. It is of interest to examine whether a relationship exists between the minimized value of the objective function and the corresponding surface-tension measurement. A correlation between these two quantities may indicate the presence of systematic error. Figure 2 shows surface-tension measurements plotted against the corresponding minimized objective-function values. For the methanol data (Fig. 2a) the objective function appears to be uncorrelated with the surface-tension measurements. The toluene data (Fig. 2b), on the other hand, do show a correlation, with low surfacetension measurements corresponding to low values of the objective function. This suggests that these measurements are systematically too high, although the measurements are already slightly lower than the literature value. The specific cause of the systematic error is not known for these measurements. The cause of the discrepancy between these measurements and the literature value could be an impurity in our samples or an error in the cited literature value. An impurity which is uniformly distributed in the sample fluid might cause a reduction in the surface tension. The presence of such an impurity could not be detected from the nonuniformity of drop shape (as indicated by the objective function value) since it would affect the surface tension uniformly over the entire drop surface. Furthermore, if surface equilibration is rapid, the impurity would have no effect on the standard deviation of the measurements.
FIG. 1. Edge coordinates from one edge of a vertical wire. The horizontal axis is greatly expanded to illustrate the noise level in the data.
AP: Colloid
665
DIGITAL IMAGE ANALYSIS OF PENDANT DROPS
tion routines and a modular program structure allows for flexibility in data analysis. The magnitude of the objective function can be used not only as a measure of the noise in fitting the experimental and predicted drop profiles, but also to determine whether there is any systematic error in the system. The improvements in the experimental apparatus and the digital image analysis discussed here permit obtaining surface tensions from the pendant drop technique with an accuracy equal to or better than other techniques such as capillary rise and the maximum bubble-pressure methods. ACKNOWLEDGMENTS The authors acknowledge Professor Charles M. Maldarelli from the City College of The City University of New York and Professor Kate J. Stebe and David O. Johnson from The Johns Hopkins University for helpful discussions and for providing a copy of the software they developed for digital image analysis of pendant drop profiles. They also acknowledge financial support for this research via NASA Grant NAG3-1278 and NSF Grant CTS-9103058.
REFERENCES 1. 2. 3. 4. 5.
FIG. 2. Surface tension measurements plotted against the corresponding minimized objective function values: (a) methanol drops, (b) toluene drops.
7. CONCLUSIONS
This study indicates that the performance of the commercial pendant drop apparatus can be improved significantly by employing a gas-tight syringe, providing for more uniform temperature control of the environmental chamber and syringe, using a regulated power supply for the light source, and mounting the apparatus on a vibration-isolation table. Convergence problems in fitting the drop profiles can be eliminated by using a highly robust minimization algorithm such as Powell’s method. Several standard algorithms which utilize derivative evaluations are found to be much faster than Powell’s method, however, they are not expected to be as robust or flexible. The use of general purpose minimiza-
/ m4585$3915
12-22-95 21:57:10
coida
6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22.
Ambwani, D., and Fort, T., Surf. Colloid Sci. 11, 93 (1979). Allen, B., Trans. Metall. Soc. AIME 227, 1175 (1963). Nagarajan, N., and Robinson, R., J. Chem. Eng. Data 31, 168 (1986). Rotenberg, Y., Boruvka, L., and Neumann, A., J. Colloid Interface Sci. 93, 169 (1983). Anastasiadis, S., Chen, J., Koberstein, J., Siegel, A., Sohn, J., and Emerson, J., J. Colloid Interface Sci. 119, 55 (1987). Jennings, J. W., and Pallas, N. R., Langmuir 4, 959 (1988). Lin, S., McKeigue, K., and Maldarelli, C., AIChE J. 36, 1785 (1990). Cheng, P., Li, D., Boruvka, L., Rotenberg, Y., and Neumann, A. W., Colloids Surf. 43, 151 (1990). Pallas, N. R., and Harrison, Y., Colloids Surf. 43, 169 (1990). Hansen, F., and Rodsrud, G., J. Colloid Interface Sci. 141, 1 (1991). Hansen, F. K., J. Colloid Interface Sci. 160, 209 (1993). Andreas, J. M., Hauser, E. A., and Tucker, W. B., J. Phys. Chem. 42, 1001 (1938). Fordham, S., Proc. R. Soc. London A 194, 1 (1948). Mills, O. S., Br. J. Appl. Phys. 3, 358 (1952). Stauffer, C. E., J. Phys. Chem. 69, 1933 (1965). Misak, M. D., J. Colloid Interface Sci. 27, 141 (1968). Padday, J. F., Philos. Trans. R. Soc. London A 269, 265 (1971). Huh, C., and Reed, R. L., J. Colloid Interface Sci. 91, 472 (1983). Girault, H. H. J., Schiffrin, D. J., and B. D. V. Smith, J. Electroanal. Chem. 137, 207 (1982). Girault, H. H. J., Schiffrin, D. J., and B. D. V. Smith, J. Colloid Interface Sci. 101, 257 (1984). Press, W., Teukolsky, S., Vetterling, W., and Flannery, B., ‘‘Numerical Recipes in C.’’ Cambridge Univ. Press, Cambridge, UK, 1992. Jasper, J. J., J. Phys. Chem. Ref. Data 1, 841 (1972).
AP: Colloid
177, 658–665 (1996)
0080
Robust Digital Image Analysis of Pendant Drop Shapes DAVID B. THIESSEN, DOMINIC J. CHIONE, CLINTON B. MCCREARY,
AND
WILLIAM B. KRANTZ 1
Department of Chemical Engineering, University of Colorado, Boulder, Colorado 80309-0424 Received May 10, 1994; accepted June 23, 1995
The advent of digital image analysis permits the use of several hundred data points defining the profile of a pendant drop in order to determine the surface tension from a fit to a solution to the Young–Laplace equation. This paper discusses several improvements to both the experimental apparatus for carrying out pendant drop measurements and to the methodology for analyzing the drop profiles in order to determine the surface tension. Several standard minimization algorithms were evaluated for fitting the theoretical and experimental drop profiles. These algorithms are very reliable, although there is a tradeoff between robustness and speed of computation. The most robust algorithm was found to converge irrespective of the ‘‘roughness’’ of the experimental drop profile or the accuracy of the initial estimates of the fitting parameters. Whereas the standard deviation of the surface-tension data provides a measure of the precision, the objective function is shown to be an indicator of systematic error in the data. Surface-tension measurements are reported for two standard liquids, methanol and toluene. A standard deviation of 0.016 dyn/cm is obtained for methanol using the full drop profile; this is significantly below the value of 0.036 dyn/cm obtained using the selected plane method, thus indicating that this improved pendant drop technique can give very precise measurements of surface tension. q 1996 Academic Press, Inc.
1. INTRODUCTION
The analysis of the shape of a pendant liquid drop is the basis for a very accurate method of measuring surface tension (1). Moreover, it offers the advantages of requiring only very small samples and being amenable to measuring dynamic interfacial tensions. This method has been applied to common liquids under ambient conditions as well as to high-temperature molten metals (2) and high-pressure liquid/vapor systems near the critical point (3). In recent years a number of groups have developed techniques to automate the pendant drop method using digital image analysis (4– 11). Rapid analysis of surface tension by computer has been applied to tracking the dynamically changing surface tension in a nonequilibrium system undergoing mass transfer (7). Computer analysis of surface tension by the pendant drop 1
To whom correspondence should be addressed.
0021-9797/96 $12.00 Copyright q 1996 by Academic Press, Inc. All rights of reproduction in any form reserved.
/ m4585$3915
12-22-95 21:57:10
method involves first locating the edge of the drop from a digital image which yields a set of points spaced along the drop profile. Next, the Young–Laplace equation for an axisymmetric interface is solved for an initial estimate of the parameters on which it depends. The parameters which govern the theoretical drop profile are adjusted and the equation is solved repeatedly until a best fit to the experimental points is obtained. The best fit is determined by minimizing an objective function which is typically a summation of the squared distances between the experimental points and the theoretical drop profile. Several problems are still encountered when using the pendant drop method to determine surface tension. The apparatus which is commercially available has some shortcomings related to operation at pressures different from ambient, temperature control, and stability of the light source used for illuminating the drop. Some methods used for determining the experimental drop profile yield ‘‘rough’’ shapes which cause problems in fitting the Young–Laplace equation. Moreover, some choices employed for the objective function used in fitting the theoretical and experimental drop profiles are not optimal. In addition, the Newton–Raphson numerical technique used to determine the fitting parameters is prone to nonconvergence for ‘‘rough’’ drop profiles or when the initial estimates of the parameters are not sufficiently accurate. Finally, although the standard deviation provides a measure of the precision of the data, no attempts have been made to assess any systematic error which might be present in the fit between the theoretical and experimental drop profiles. This paper addresses all of the aforementioned problems in an attempt to improve the precision of the pendant drop method for measuring surface tension. Moreover, this presents the first implementation of software to analyze pendant drop data using the Apple/Macintosh operating system. This paper is organized as follows. First, a brief review of prior studies will be given along with a critique of these studies which elaborates on some of the above mentioned problems. This is followed by a summary of the experimental apparatus which incorporates improvements which have been made to the commercially available instrument. The experimental procedure and data analysis then are summa-
658
coida
AP: Colloid
DIGITAL IMAGE ANALYSIS OF PENDANT DROPS
rized; the latter incorporates several improvements which have been made in the edge-detection method, calibration, and fitting algorithm; in particular, several standard minimization algorithms for determining the fitting parameters (which include the surface tension) are evaluated. The most robust algorithm was found to converge irrespective of the roughness of the drop profile or the accuracy of the initial guesses for the fitting parameters. 2. PRIOR STUDIES
659
to the fitting parameters. The method also prevented large oscillations in the parameters as the minimum was approached. Lin et al. (7) used a system having a resolution of 480 1 512 pixels with 256 levels of gray. These investigators fully utilized this high resolution in their edge detection routine and employed a local threshold and interpolation technique to provide smoother drop profiles; the local threshold was chosen to be halfway between the local maximum and minimum light intensity near the drop interface. However, they still encountered convergence problems owing to their use of a Newton–Raphson method to determine the four parameters involved in fitting the Young–Laplace equation to the experimental drop profile. The recent work of Hansen and Rodsrud (10) employed a system with a resolution of 256 1 256 pixels with 64 levels of gray. Their edge detection routine was similar to that used by Lin et al. (7). However, the objective function in their minimization routine used the square of the horizontal rather than the minimum distance between the theoretical and experimental drop profiles; as such, this method incurred more error near the apex of the drop where the horizontal distance can be quite significant. Nonetheless, their standard deviations were quite small, 0.01 to 0.06 dyn/ cm although some of their surface tensions differed significantly from currently accepted values; this could indicate some type of systematic error. Hansen (11) has improved on the method of Hansen and Rodsrud by utilizing the sum of squared normal distances as an objective function. This brief review indicates that image digitization can significantly increase the accuracy of the pendant drop technique by employing many points to determine the experimental drop shape from which the surface tension can be extracted. Efficient routines for determining the fitting parameters are based on least-squares minimization employing an objective function based on the square of the minimum distance between the experimental and predicted drop profile. Determination of the fitting parameters which include the surface tension via Newton–Raphson methods is complicated by convergence problems owing to ‘‘rough’’ drop profiles or inaccurate estimates of the initial parameter values. Whereas the standard deviation provides a measure of the scatter in the data, systematic errors can result in anomalously high or low values of the surface tension determined by the pendant drop method; clearly it would be of value to develop some measure of any such systematic errors. Although personal computers are now used to implement the image digitization and profile fitting, no software has been developed for the Apple/Macintosh operating system which offers a very powerful user-friendly graphical interface.
Early applications of the pendant drop technique by Andreas et al. (12) employed either the inflection plane or the selected plane method; these methods correlated the surface tension in terms of a drop-shape factor. Subsequent studies extended the tables for the drop-shape factor (13–17). These methods employed relatively few data points taken from the drop shape and hence were less accurate than other methods for determining the surface tension. The advent of image digitization permitted more accurate fitting of the drop shape using mainframe computers (18), thereby allowing determination of the surface tension by fitting the Young–Laplace equation to the full drop profile; however, the digitization was done on enlarged photographs of the drop. The first use of direct digitization of the drop without using any photographs was done by Girault et al. (19). Rapid advances in image digitization permitted resolution as high as 1066 1 575 pixels (20). Rotenberg et al. (4) improved the method of fitting the Young–Laplace equation to the data by employing rigorous least-squares minimization of an objective function equal to the sum of the squares of the minimum distance between the predicted and experimental drop profile. This minimization involves four parameters, the coordinates and radius of curvature of the apex of the drop and the Bond number which contains the unknown surface tension. These parameters are determined by a Newton-Raphson method employing exact derivatives of the objective function and incremental loading. Anastasiadis et al. (5) carried out image digitization and fitting of the Young–Laplace equation for the drop shape using a personal computer. They used a global thresholding method together with a smoothing algorithm for edge detection. They utilized a robust shape comparison method for fitting the Young–Laplace equation to the drop profile. This method is less sensitive to large errors in individual edge points. Optimization was reduced to a single shape parameter and carried out by exhaustive search. The standard deviation of their data, 0.3 dyn/cm, was larger than that typically obtained for the generally less accurate method of selected plane (0.1 dyn/cm). Jennings and Pallas (6) made a significant improvement in 3. EXPERIMENTAL APPARATUS the efficiency of the fitting method by utilizing a modified Newton–Raphson method which did not require evaluation of A commercial pendant drop apparatus (Rame-Hart, Inc.) the second derivatives of the objective function with respect with several modifications was used in this work. The appa-
/ m4585$3915
12-22-95 21:57:10
coida
AP: Colloid
660
THIESSEN ET AL.
ratus as purchased consisted of an optical rail on which was mounted a three-axis stage supporting an environmental chamber and heated syringe assembly, a three-axis camera support stage, and a tungsten filament lamp with a variac power supply. Temperature control of the environmental chamber could be accomplished by circulation of water at a controlled temperature through a jacket or by resistive heating with feedback control. Temperature control of the syringe could be accomplished only by resistive heating regulated by a separate controller. This resulted in temperature differences between the environmental chamber and the syringe. Consequently the syringe housing was modified by replacing the resistive heater with a water jacket and circulating water from the same temperature-controlled bath through both the chamber jacket as well as this syringe jacket. This resulted in a negligible temperature difference between the pendant drop and the chamber walls and maintained the temperature constant to within {0.17C. The syringe housing also was modified to accommodate a Hamilton (Model 1001) gas-tight syringe. A gas-tight syringe is necessary when working with reduced chamber pressures and to eliminate cooling effects owing to evaporation of volatile liquids. The final modification of the apparatus involved replacing the lamp variac with a regulated power supply which reduced voltage drift from 0.2 V to less than 0.05 V. The pendant-drop apparatus rests on a bench-top vibration-isolation table (Ehrenreich Photo-Optical Industries Model 78221), which consists of massive steel plates floating on three air bladders. Water hoses and electrical cords which are attached to the apparatus are supported from springs suspended from the lab ceiling to reduce vibrations transmitted to the apparatus. A Panasonic CCD camera with a Nikon lens (micro NIKKOR 55 mm) and a 27.5-mm extension tube are used for imaging the drop. A digital image of a pendant drop is acquired using an image-capture card and associated software (Video Image 1200, Scion Corp.) on a Macintosh IIfx personal computer. 4. EXPERIMENTAL PROCEDURE
First, the air table and the optical stand must be adjusted to ensure that a pendant drop will hang exactly perpendicular with respect to the camera and the light source. Analysis via the technique developed here then requires that a column of pixels in the image be parallel to the gravity vector. The apparatus is aligned using a plumb bob consisting of a weight hanging from a fine copper wire attached to the apparatus. The weight is submerged in a beaker of water to damp oscillations. The plumb bob wire is imaged with the CCD camera. The level of the apparatus is adjusted until the edge of the wire in the image aligns with a column of pixels. This leveling procedure must be carried out periodically because
/ m4585$3915
12-22-95 21:57:10
coida
of changes in the volume of the air bladders that support the vibration-isolation table. Before surface-tension measurements are made, a calibration image of a 3.0 { 0.0025 mm precision ball bearing (Industrial Tectonics) is captured. This calibration yields the dimensions of a pixel. The horizontal dimension of a pixel, lx , is not necessarily the same as its vertical dimension, lz , owing to an aspect-ratio error which could be caused by the camera and/or the image-capture system. This calibration involves suspending the ball bearing at the same location as the pendant drop; the camera adjustments and light-source settings are not changed after the calibration image is made. Surface-tension measurements were made on HPLC-grade methanol and toluene without further purification. The syringe was cleaned and dried and then rinsed several times with the liquid to be measured. A pool is made on the bottom of the environmental chamber with liquid from the final rinse. The syringe then is filled and placed in the syringe housing. Water from a constant temperature bath then is circulated through the chamber and syringe housing jackets; the liquid and air in the chamber are allowed to equilibrate for 1–2 h. The circulating bath is turned off for 60 s before a pendant drop image is captured to allow vibrations to damp out. The room lights are then turned off and the image is captured by the computer directly from the CCD camera; this yields considerably improved image quality over capturing a videotaped image. After a drop has been imaged, it is allowed to fall to the bottom of the chamber and a new drop is formed. The bath is turned on for 60 s and then the image capture sequence is repeated. Typically ten different drops are imaged in a single run. The digital images of the calibration sphere and pendant drops are stored as graphics files on the computer hard drive until they are analyzed. 5. DATA ANALYSIS
A typical experimental run yields one graphics file containing an image of the calibration sphere and a number of graphics files containing pendant drop images. The images are in 560 1 420 pixel format with 8 bits per pixel gray scale resolution (256 shades of gray). Image-analysis software was developed on the Macintosh IIfx computer in the C programming language. Picture files stored in the PICT, TIFF, or Encapsulated PostScript (EPSF) formats can be imported by the program. The precision sphere image is analyzed to obtain calibration parameters which are used in the subsequent analysis of the pendant drop images for surface-tension calculation. A subroutine was written to open the graphics file containing the image data, read the size of the image, read the pixel intensities, and store these values in a large two-dimensional array. The edge of an object defined by a digital image is found by searching along rows or columns of pixels looking for an abrupt intensity change. The sides of the sphere or drop
AP: Colloid
661
DIGITAL IMAGE ANALYSIS OF PENDANT DROPS
are found by searching along rows while the bottom is found by searching along columns. An intensity difference between adjacent pixels greater than some specified value indicates the vicinity of the interface. Once the approximate position of the interface has been located, a local threshold intensity is determined, which is taken to define the ‘‘true’’ edge. The local threshold is defined to be the average of 20 pixels in the same row or column, 10 on each side of the interface, starting 5 pixels away from the approximate edge position. The fractional coordinate defining the edge of the object along the search direction is then found by linear interpolation between adjacent pixels which bound the threshold intensity. Successive edge coordinates defining both sides of the axisymmetric object are then stored in an array. This technique locates the interface within less than one pixel with floating point/integer profile coordinates rather than the integer/integer profile coordinates used in earlier less accurate techniques. A subroutine called CALIBRATE analyzes the image of a precision sphere with known radius, Rs , in order to determine the dimensions of a pixel. The horizontal and vertical pixel dimensions, lx and lz , are determined by fitting a circle to the edge coordinates (xi , zi ). An objective function defined as M
Fc Å ∑ {[l 2x (xi 0 x0 ) 2 / l 2z (zi 0 z0 ) 2 ] 1 / 2 0 Rs } 2
[1]
i Å1
is minimized with respect to the four parameters lx , lz , x0 , and z0 , where the latter two are the pixel coordinates of the center of the circle being fit to the data. The parameters are initialized from the analytical solution for a circle which passes through three widely spaced edge points. A minimization subroutine then is called to minimize the objective function given by Eq. [1]; these subroutines will be discussed following a description of how the surface tension is determined. To analyze for the surface tension of a pendant drop, the edge of the drop is located as described above. In order to provide good initial estimates of the parameters, the surface tension is first estimated using the method of selected plane. The equatorial diameter, de , of the pendant drop is found by searching for a maximum diameter starting from the top of the drop. The selected plane diameter, ds , corresponding to the diameter of the drop at a location equal to one equatorial diameter above the drop apex, is then determined. These measurements are corrected for the aspect ratio error determined from the calibration using the precision sphere. The ratio S å ds /de is correlated to a shape factor, 1/H, using Misak’s (16) equation. The surface tension, g, is determined from the equation gÅ
/ m4585$3915
Drgd 2e , H
12-22-95 21:57:10
[2]
coida
where Dr is the density difference between the drop and its surrounding fluid and g is the gravitational acceleration. This initial estimate of the surface tension then is used in a procedure to determine the surface tension by fitting a solution to the Young–Laplace equation to the experimentally determined profile coordinates; approximately 1000 points were used to determine the drop profile. The Young– Laplace equation relates the curvature of an interface to the pressure difference across that interface. In a coordinate system whose origin is at the drop apex, the Young–Laplace equation for an axisymmetric pendant drop can be recast as three coupled first-order differential equations [Rotenberg et al. (4)], df sin f Å 2 0 bz 0 ds x dx Å cos f ds dz Å sin f, ds
[3]
where b å DrgR 20 / g is the Bond number; x and z are the horizontal and vertical coordinates, respectively, nondimensionalized with respect to R0 , the radius of curvature at the drop apex; s is the arc length along the drop surface measured from the drop apex nondimensionalized with respect to R0 ; and f is the angle between the surface tangent and the horizontal. The differential equations described by [3] must satisfy the conditions that f Å x Å z Å 0 at s Å 0. An initial estimate of b is obtained from a correlation between b and the diameter ratio S [Hansen and Rodsrud (10)], the latter being obtained from the method of selected plane. The differential equations given by [3] are integrated by expanding each of the variables f, x, and z in a four-term Taylor series starting at s Å 0 and propagating forward in s; the values of the first, second, and third derivatives in this Taylor series are obtained from Eq. [3], evaluating them at the point of interest. The experimental edge coordinates in terms of pixels, (xi , zi ), are converted to dimensionless form as ji Å
lx (xi 0 x0 ) , R0
zi Å
lz (zi 0 z0 ) , R0
[4]
where lx and lz are the horizontal and vertical pixel dimensions, respectively, determined from the calibration; x0 and z0 now denote the pixel coordinates of the drop apex; and R0 is the radius of curvature at the drop apex. The four parameters b, x0 , z0 , and R0 are varied in order to minimize the difference between the theoretical and experimental drop shapes. An objective function of the form used by Rotenberg et al. (4) is minimized,
AP: Colloid
662
THIESSEN ET AL. M
Fd (au ) Å R 20 ∑ [d(ui , £: au )] 2 ,
[5]
i Å1
where au denotes the parameter vector. The operator d represents the minimum nondimensional distance between an experimental point ui Å ( ji , zi ) and the theoretical curve £ Å [x(s), z(s)] determined from Eq. [3]. The objective function for the pendant drop is minimized using the same general purpose minimization algorithm as is used in the calibration procedure described above. After the objective function has been minimized, the surface tension is calculated from the optimized values of the Bond number and radius of curvature at the drop apex, R0 , according to the following identity: gå
DrgR 20 . b
[6]
The minimum value of the objective function is recorded along with the surface tension in the output from the program. Efficient methods of minimizing objective functions of the form given by Eq. [5] are based on a Taylor series expansion of Fd about a given point in parameter space, au 0 : Fd (au ) É Fd (au 0 ) / (au 0 au 0 )rÇ Fd (au 0 ) /
1 (au 0 au 0 )(au 0 au 0 ): ÇÇ Fd (au 0 ). 2
[7]
Assuming that Fd is well approximated by such a quadratic form near the minimum, and that the derivatives can be computed, an iterative scheme for optimizing the parameter vector au is given by au k /1 Å au k 0 [ÇÇ Fd (au k )] 01rÇ Fd (au k ),
[8]
where [ÇÇ Fd (au k )] 01 is known as the inverse Hessian matrix. This is the well-known Newton–Raphson method, which is very efficient when the objective function is of quadratic form but is susceptible to nonconvergence in cases where it is not. At the opposite extreme are highly robust algorithms which make no assumptions about the form of the objective function but simply march steadily downhill toward the minimum. A number of standard algorithms exist which constitute a compromise between these two extremes. Minimization algorithms from three categories have been evaluated in this work, including algorithms which utilize: (1) function evaluations only, (2) function evaluations and first partial derivatives of Fd with respect to the parameters, (3) function evaluations and first and second partial derivatives of Fd with respect to the parameters.
/ m4585$3915
12-22-95 21:57:10
coida
Computation speed is found to be greatly enhanced by using exact derivative evaluations (6). Although algorithms from the first category are generally slower than those from the latter two categories, they do have significant advantages in terms of flexibility and robustness. In particular, it is sometimes advantageous to use a ‘‘robust’’ fitting function as opposed to the sum of squared residuals given by Eq. [5]. Robust objective functions generally are not as smooth as that given by Eq. [5] and may have discontinuities in the derivatives which preclude the use of minimization algorithms in categories 2 and 3. A modular program has been created in order to provide maximum flexibility in image analysis. The minimization algorithms are standard off-the-shelf routines which are very efficient and reliable. The implementation of these algorithms is that given by Press et al. (21) with only minor modifications. These algorithms are completely generalized so that they can be passed the name of the subroutine which evaluates the objective function along with the number of parameters on which to optimize. An algorithm known as Powell’s method belongs to the first category of algorithms listed above. When Fd is well approximated by a quadratic form then surfaces of constant Fd near the minimum are hyperellipsoids. When searching for the minimum, it is most efficient to search along the directions of the principal axes of the ellipsoid. For a quadratic function of N variables the minimum can be found with only N line minimizations in this manner. The directions of the principal axes thus are called ‘‘noninterfering’’ or ‘‘conjugate’’ directions. In the more general case where Fd is not necessarily of quadratic form, Powell’s method attempts to find noninterfering directions and repeats the cycle of N line minimizations until it converges to a minimum. This algorithm has never failed to converge to a minimum of the objective function in the many cases which were run. A method which uses first derivatives (category 2) is the Broyden–Fletcher–Goldfarb–Shanno (BFGS) variable metric method. This method is a quasi-Newton method in that it uses gradient information to iteratively approximate the inverse Hessian matrix from which new estimates of the optimal parameter vector are computed according to a Newton-type formula (Eq. [8]). The first derivatives of the objective function are given by ÌFd Å2 Ìaj
M
∑ ri i Å1
Ìri Ìaj
j Å 1, . . . , N,
[9]
where ri å R0 d(ui ,£ : au ). The partial derivatives of the individual residuals were evaluated in the manner of Jennings and Pallas (6). The algorithm uses a backtracking scheme to avoid stepping too far in the Newton direction for cases where the function is poorly approximated by a quadratic form.
AP: Colloid
663
DIGITAL IMAGE ANALYSIS OF PENDANT DROPS
The third method which was evaluated is specifically intended for nonlinear least squares minimization. The Levenberg–Marquardt method uses evaluations of the Hessian matrix (category 3) as well as the gradient to perform minimization. The Hessian matrix is approximated as Ì 2 Fd É2 Ìaj Ìak
M
∑ i Å1
Ìri Ìri Ìaj Ìak
j, k Å 1, . . . , N.
[10]
This approximation neglects a term with second derivatives of ri . Inclusion of the second derivative term can actually cause instability in the procedure (21). This algorithm uses a combination of Newton’s method (Eqs. [8–10]) and the steepest descent method. The steepest descent method consists of taking a step in the direction of the negative gradient of Fd . The algorithm uses the steepest descent method far from the minimum and switches smoothly to Newton’s method as the minimum is approached. The computation speed of the three algorithms was compared for a typical pendant drop with 934 edge points. The CPU times (on a Power Macintosh 7100/80 AV) necessary to minimize the objective function were 66 seconds for Powell’s method, 4.2 seconds for the variable metric method, and 0.57 seconds for the Levenberg–Marquardt method. Limited tests indicate that each of these methods converge reliably even in cases where errors as high as 100% in the initial parameter estimates were introduced. Following minimization, it is possible to estimate the uncertainty in the fit parameters from the Hessian matrix evaluated at the minimum of the objective function (6). The confidence region is bounded by the ellipsoid given by (au 0 au m )(au 0 au m ): ÇÇ Fd (au m ) Å DFd ,
[11]
where au m is the optimal parameter vector and DFd is given by DFd Å 2s 2NF(N, M 0 N)
where s 2 is the variance [Fd (au m )/(M 0 N)], and F(N, M 0 N) is the F-test statistic. The application of Eq. [11] to determine the confidence ellipsoid is strictly valid only in the case where the errors are Gaussian distributed. Taking the confidence intervals of the individual parameters to be the extreme values of the ellipsoid, the confidence interval for the surface tension can be determined from Eq. [6]. At the 95% confidence level the confidence interval for the surface tension of a typical water drop was {0.13 dyn/cm. This should be a conservative estimate of the error in these measurements. 6. RESULTS
Surface-tension measurements were made on methanol and toluene at 207C in air saturated with the liquid’s vapor.
/ m4585$3915
12-22-95 21:57:10
coida
TABLE 1 Surface Tension Measurements for Methanol in Air at 207C Surface tension (dyn/cm)
Drop number 1 2 3 4 5 6 7 8 9 10 Mean Standard deviation Percent deviation
Selected plane
Profile fitting
Dimensionless objective function
22.314 22.312 22.243 22.226 22.249 22.222 22.263 22.217 22.231 22.231 22.251 0.036 0.16
22.452 22.444 22.454 22.457 22.457 22.434 22.422 22.456 22.425 22.418 22.439 0.016 0.069
0.0036 0.0035 0.0041 0.0045 0.0045 0.0046 0.0040 0.0044 0.0036 0.0036 0.0041 0.0004 10.5
The density difference between the liquid and gas phases was taken as 0.7904 g/cm3 for methanol and 0.8656 g/cm3 for toluene. The vapor-phase density was calculated from the Ideal Gas Law using Raoult’s Law to estimate the composition of the gas phase. Table 1 gives the measured values of surface tension for ten separate drops of methanol using both the selected plane and profile-fitting methods. The average value of the surface tension determined by the profile-fitting method is 22.44 { 0.02 dyn/cm. This value falls within the uncertainty limits of the accepted literature value of 22.5 { 0.1 dyn/cm which was obtained by the capillary-rise method (22). The selected plane method yields a value almost 1% lower than the profile-fitting method and gives a larger standard deviation [22.25 { 0.04 dyn/cm]. The objective function values reported in Tables 1 and 2 were made dimensionless by dividing Eq. [5] by R 20 . The results for toluene are presented in Table 2. Again, ten separate drops were measured by the two methods mentioned above. The average value of the surface tension determined by the profile-fitting method is 28.20 { 0.07 dyn/cm. This is 1% lower than the accepted literature value of 28.5 { 0.1 dyn/cm obtained by the maximum bubble-pressure method (22). The selected plane method gives a surface-tension value 0.6% lower than the profile-fitting method [28.04 { 0.07 dyn/cm], however the two methods have the same standard deviation for the toluene drops. The standard deviation is considerably larger for the measurements of toluene drops than for the methanol drops. It is proposed here that the minimized value of the objective function be used to indicate the quality of a surfacetension measurement. One would like to know whether a given value of the objective function indicates a mismatch between the theoretical and experimental drop shapes or
AP: Colloid
664
THIESSEN ET AL.
TABLE 2 Surface Tension Measurements for Toluene in Air at 207C Surface tension (dyn/cm)
Drop number 1 2 3 4 5 6 7 8 9 10 Mean Standard deviation Percent deviation
Selected plane
Profile fitting
Dimensionless objective function
27.949 28.011 28.053 28.099 28.149 28.089 27.903 28.082 28.003 28.034 28.037 0.074 0.263
28.105 28.173 28.149 28.132 28.304 28.277 28.146 28.297 28.148 28.230 28.196 0.074 0.263
0.0049 0.0052 0.0055 0.0054 0.0071 0.0063 0.0051 0.0060 0.0052 0.0057 0.0056 0.0007 11.67
whether the value arises strictly from ‘‘noise’’ in the digital imaging system. The noise level can be seen clearly by examining the located edge points from an image of a plumb wire to which the camera has been carefully aligned (Fig. 1). The horizontal axis of the plot in Fig. 1 is greatly expanded in order to show the noise level in the data. The x-coordinates of the edge points have a range of 0.225 pixels and a standard deviation of 0.04 pixels. The calibration procedure described in Section 4 involves imaging a precision sphere, locating the edge, and minimizing the objective function defined by Eq. [1]. The minimized objective function found for a precision sphere can serve as a reference value to compare with the objective function values obtained for drops. By using a precision sphere with a diameter close to the equatorial diameter of a given pendant drop, roughly the same number of image points are found for the sphere and the drop. This allows for a direct comparison of the objective-function values. The magnitude of the objective function for a precision sphere will indicate the noise level of the imaging system. A typical value of the minimized objective function for a calibration sphere (nondimensionalized with respect to R 2s ) is 1.6 1 10 04 . This value is approximately 30 times smaller than the objective function values obtained for drops in this work. The objective-function values for drops in this work are obviously above the noise level, indicating a mismatch between the experimental and theoretical drop shapes. This mismatch could arise from camera misalignment, drop vibration, needle-tip imperfections, or surface-tension gradients on the drop. Surface-tension gradients could arise from nonuniform temperatures or a nonuniform distribution of impurities on the drop surface. One significant source of error in the pendant drop measurements arises from vibrations; although the vibration table greatly helps in isolating the pendant drop apparatus from room vibrations, it is not
/ m4585$3915
12-22-95 21:57:10
coida
massive enough to eliminate all vibrations. The pendant drop apparatus should be mounted on a far more massive vibration isolation table to minimize this source of error. Future work will concentrate on improving the objective-function values by identifying and correcting other causes of error in the measurements. It is of interest to examine whether a relationship exists between the minimized value of the objective function and the corresponding surface-tension measurement. A correlation between these two quantities may indicate the presence of systematic error. Figure 2 shows surface-tension measurements plotted against the corresponding minimized objective-function values. For the methanol data (Fig. 2a) the objective function appears to be uncorrelated with the surface-tension measurements. The toluene data (Fig. 2b), on the other hand, do show a correlation, with low surfacetension measurements corresponding to low values of the objective function. This suggests that these measurements are systematically too high, although the measurements are already slightly lower than the literature value. The specific cause of the systematic error is not known for these measurements. The cause of the discrepancy between these measurements and the literature value could be an impurity in our samples or an error in the cited literature value. An impurity which is uniformly distributed in the sample fluid might cause a reduction in the surface tension. The presence of such an impurity could not be detected from the nonuniformity of drop shape (as indicated by the objective function value) since it would affect the surface tension uniformly over the entire drop surface. Furthermore, if surface equilibration is rapid, the impurity would have no effect on the standard deviation of the measurements.
FIG. 1. Edge coordinates from one edge of a vertical wire. The horizontal axis is greatly expanded to illustrate the noise level in the data.
AP: Colloid
665
DIGITAL IMAGE ANALYSIS OF PENDANT DROPS
tion routines and a modular program structure allows for flexibility in data analysis. The magnitude of the objective function can be used not only as a measure of the noise in fitting the experimental and predicted drop profiles, but also to determine whether there is any systematic error in the system. The improvements in the experimental apparatus and the digital image analysis discussed here permit obtaining surface tensions from the pendant drop technique with an accuracy equal to or better than other techniques such as capillary rise and the maximum bubble-pressure methods. ACKNOWLEDGMENTS The authors acknowledge Professor Charles M. Maldarelli from the City College of The City University of New York and Professor Kate J. Stebe and David O. Johnson from The Johns Hopkins University for helpful discussions and for providing a copy of the software they developed for digital image analysis of pendant drop profiles. They also acknowledge financial support for this research via NASA Grant NAG3-1278 and NSF Grant CTS-9103058.
REFERENCES 1. 2. 3. 4. 5.
FIG. 2. Surface tension measurements plotted against the corresponding minimized objective function values: (a) methanol drops, (b) toluene drops.
7. CONCLUSIONS
This study indicates that the performance of the commercial pendant drop apparatus can be improved significantly by employing a gas-tight syringe, providing for more uniform temperature control of the environmental chamber and syringe, using a regulated power supply for the light source, and mounting the apparatus on a vibration-isolation table. Convergence problems in fitting the drop profiles can be eliminated by using a highly robust minimization algorithm such as Powell’s method. Several standard algorithms which utilize derivative evaluations are found to be much faster than Powell’s method, however, they are not expected to be as robust or flexible. The use of general purpose minimiza-
/ m4585$3915
12-22-95 21:57:10
coida
6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22.
Ambwani, D., and Fort, T., Surf. Colloid Sci. 11, 93 (1979). Allen, B., Trans. Metall. Soc. AIME 227, 1175 (1963). Nagarajan, N., and Robinson, R., J. Chem. Eng. Data 31, 168 (1986). Rotenberg, Y., Boruvka, L., and Neumann, A., J. Colloid Interface Sci. 93, 169 (1983). Anastasiadis, S., Chen, J., Koberstein, J., Siegel, A., Sohn, J., and Emerson, J., J. Colloid Interface Sci. 119, 55 (1987). Jennings, J. W., and Pallas, N. R., Langmuir 4, 959 (1988). Lin, S., McKeigue, K., and Maldarelli, C., AIChE J. 36, 1785 (1990). Cheng, P., Li, D., Boruvka, L., Rotenberg, Y., and Neumann, A. W., Colloids Surf. 43, 151 (1990). Pallas, N. R., and Harrison, Y., Colloids Surf. 43, 169 (1990). Hansen, F., and Rodsrud, G., J. Colloid Interface Sci. 141, 1 (1991). Hansen, F. K., J. Colloid Interface Sci. 160, 209 (1993). Andreas, J. M., Hauser, E. A., and Tucker, W. B., J. Phys. Chem. 42, 1001 (1938). Fordham, S., Proc. R. Soc. London A 194, 1 (1948). Mills, O. S., Br. J. Appl. Phys. 3, 358 (1952). Stauffer, C. E., J. Phys. Chem. 69, 1933 (1965). Misak, M. D., J. Colloid Interface Sci. 27, 141 (1968). Padday, J. F., Philos. Trans. R. Soc. London A 269, 265 (1971). Huh, C., and Reed, R. L., J. Colloid Interface Sci. 91, 472 (1983). Girault, H. H. J., Schiffrin, D. J., and B. D. V. Smith, J. Electroanal. Chem. 137, 207 (1982). Girault, H. H. J., Schiffrin, D. J., and B. D. V. Smith, J. Colloid Interface Sci. 101, 257 (1984). Press, W., Teukolsky, S., Vetterling, W., and Flannery, B., ‘‘Numerical Recipes in C.’’ Cambridge Univ. Press, Cambridge, UK, 1992. Jasper, J. J., J. Phys. Chem. Ref. Data 1, 841 (1972).
AP: Colloid