Langmuir 2005, 21, 11241-11250
11241
Applying Tikhonov Regularization to Process Pendant Droplet Tensiometry Data Y. Leong Yeow,* Krista L. McGowan, and S. Ian Ong Department of Chemical and Biomolecular Engineering, The University of Melbourne, Victoria 3010, Australia
Yee-Kwong Leong School of Engineering, James Cook University, Townsville, Queensland 4811, Australia Received June 20, 2005. In Final Form: September 6, 2005 The problem of obtaining the first and second derivatives of the profile of a pendant droplet is formulated as an integral equation of the first kind. This equation is solved by Tikhonov regularization in which the method of general cross validation is used to guide the selection of the regularization parameter. These derivatives are converted into mean curvature as a function of droplet height. Surface tension is then obtained by regression computation between the mean curvature and two possible algebraic expressions suggested by the Laplace-Young equation. This way of obtaining surface tension is demonstrated by applying it to a number of published droplet profiles. Some of the problems encountered are discussed and solutions suggested.
1. Introduction Pendant droplet tensiometry is a well-established method for measuring surface and interfacial tension.1-5 The general physical principle behind the method is comparatively simple. A single droplet of the liquid under test is carefully formed at the tip of a vertical capillary. The equilibrium shape of the resulting axisymmertic droplet is determined by the balance between surface tension and gravitational body forces. By measuring the shape of the droplet and knowing the density of the liquid it is then possible, by solving the differential equation that describes the droplet shape, to back-calculate the surface or interfacial tension.2 CCD cameras interfaced with microcomputers are now routinely used to capture and store droplet images. Image processing software is then used to convert these images into droplet profiles, i.e., plot of the variation of droplet radius, rM(z), with vertical distance, z. Superscript M is used here to distinguish the measured droplet profile from the computed profile, rC(z). Reliable experimental techniques and efficient data handling procedures have been developed for these new generation of computer-based tensiometers.6-11 There are now tensiometers on the * To whom correspondence should be addressed. E-mail: yly@ unimelb.edu.au. Fax: 61 3 8344 4153. (1) Bashforth, F.; Adams, J. C. An attempt to test the theories of Capillary Action; Cambridge University Press: Cambridge, 1883. (2) Padday, J. F. Surface Colloid Sci. 1969, 1, 101. (3) Rotenberg, Y.; Boruvka, L.; Neumann, A. W. J. Colloid Interface Sci. 1983, 93, 169. (4) Adamson, A. W.; Gast, A. P. Physical Chemistry of Surfaces, 6th ed.; John Wiley: New York, 1997. (5) Lin, S.-Y.; Chen, L.-J.; Xyu, J.-W.; Wang, W.-J. Langmuir 1995, 11, 4159. (6) Satherley, J.; Girault, H. H. J.; Schiffrin, D. J. J. Colloid Interface Sci. 1990, 136, 574. (7) Thiessen, D. B.; Chione, D. J.; McCreary, C. B.; Krantz, W. B. J. Colloid Interface Sci. 1996, 177, 658. (8) Song, B.; Springer, J. J. Colloid Interface Sci. 1996, 184, 64. (9) Song, B.; Springer, J. J. Colloid Interface Sci. 1996, 184, 77. (10) Emelyanenko, A. M.; Boinovich, L. B. Colloids Surf., A 2001, 189, 197. (11) Zhou, Y. Z.; Gaydos, J. J. Adhesion 2004, 80, 1017.
market that will perform all these steps almost automatically and generate rM(z) at the press of a button. A number of computational schemes have been developed to convert the rM(z) generated by computer-based tensiometers to surface or interfacial tension. The various schemes currently in common use may differ in details; the essential principles behind them are, however, similar. The governing equation is the Laplace-Young equation, which describes the force balance across the droplet surface2
(
) ( )
1 2 1 ∆Fg + ) z RApex R1 R2 γ
(1)
or
∆Fg 2 z RApex γ
( )
2H(z) )
(1b)
R1 and R2 are the local principal radii of curvature of the surface and H(z) ≡ 1/2[1/R1 + 1/R2] is the mean curvature. R1 and R2 vary with the vertical position, z, in the droplet and are given by12
[ ( )] [ ( )]
R 1 ) rC × 1 + 1+
R2 ) -
drC dz
drC dz d2rC dz2
2 1/2
(2a)
2 3/2
(2b)
These are standard results of differential geometry independent of eq 1. At the apex of the droplet z ) 0 and R1 ) R2 ) RApex. ∆Fg is the known specific gravitational body force acting on the droplet. Surface tension, γ, is the only unknown material property in eq 1. (12) Pearson, C. E. Handbook of Applied Mathematics, 2nd ed.; van Nostrand Reinhold: New York, 1983.
10.1021/la051650p CCC: $30.25 © 2005 American Chemical Society Published on Web 10/14/2005
11242
Langmuir, Vol. 21, No. 24, 2005
Combining eqs 1 and 2 results in a differential equation for the computed profile rC(z) with ∆Fg/γ appearing as an unknown but constant coefficient. This differential equation is nonlinear and does not have a closed-form solution, but for any assumed value of ∆Fg/γ together with the relevant boundary conditions, it can be integrated numerically using a variety of techniques. These include the Adams-Bashforth method used to solve this problem in the original investigation of Bashforth and Adams1 and the different variants of the Runge-Kutta procedure.3,7,13 More recently, the Galerkin finite element method has also been applied to solve this equation.14 Surface tension is then determined by iteratively adjusting the parameter ∆Fg/γ so that the computed radius, rC(z), closely mirrors its measured counterpart, rM(z). Here again, different numerical techniques have been developed to locate the γ that minimizes the difference between rM(z) and rC(z). These methods may differ in the way the difference between rM(z) and rC(z) is quantified and the strategy adopted to minimize this difference, but most can be reduced to some form of iterative nonlinear least-squares regression computation.3,7-9,13-15 As can be seen from the above general description, irrespective of way the computation is performed, the determination of surface tension is treated as a problem of determining the unknown coefficient of a differential equation so that its solution matches with a set of experimental data. Such a problem is often referred to as an inverse problem to distinguish it from the closely related but simpler problem, known as the direct problem, where the coefficient is known and the task is to solve the equation for the dependent variable.16 For the specific case of eq 1, the direct problem is to solve this equation to determine the shape of the pendant droplet for a liquid with known surface tension and density. The corresponding inverse problem is to back-calculate the surface tension from a measured droplet profile. The inverse problem is ill-posed in the sense that two nearby measured droplet profiles of a given liquid, differing from each other possibly only at the level of measurement error bars, if inappropriate methods are used to handle the inverse problem the end result may be two significantly different surface tensions, γ. This is in complete contrast with the well-posed direct problem, where for two liquids with material property group ∆Fg/γ close to one another and the same boundary conditions, the solutions of eq 1 will result in two droplet shapes that are also close to one another.16,17 Since the mathematical nature of an inverse problem is very different from its direct counterpart, a specialized class of numerical procedures has been developed to deal with inverse problems. The Tikhonov regularization procedure adopted in the present investigation is a wellproven method from this class of specialized procedures. This represents a significant new development in pendant droplet tensiometry. Since the method allows for the inverse nature of the pendant droplet problem, it can be expected to lead to more reliable results with less computational efforts.16,17 The main source of the difficulty with the inverse problem can be traced to the amplification of the unavoidable noise in the measured droplet profile. The (13) de Ramos, A. L.; Redner, R. A.; Cerro, R. L. Langmuir 1993, 9, 3691. (14) Dingle, N. M.; Tjiptowidjojo, K.; Basaran, O. A.; Harris, M. T. J. Colloid Interface Sci. 2005, 286, 647. (15) Jennings, J. W.; Pallas, N. R. Langmuir 1988, 4, 959. (16) Groetsch, C. W. Inverse Problems in Mathematical Sciences; Vieweg: Braunschweig, 1993. (17) Engl, H. W.; Hanke, M.; Neubauer, A. Regularization of Inverse Problems; Kluwer: Dordrecht, 2000.
Yeow et al.
standard procedure for dealing with such noise is to apply various types of smoothing, such as least-squares splines, to smooth out the noise.13 While this will remove noise, it is just as likely to filter out some of the essential information embedded in the data. Tikhonov regularization, with its built-in regularization parameter, λ, is able to achieve an appropriate balance between data smoothing and information retention.17 It will be applied to solve the ill-posed problem of obtaining d2rC/dz2 and drC/dz from the measured rM(z) while keeping noise amplification under control. These derivatives are then converted into R1 and R2 via eq 2. Finally, they are used, in conjunction with eq 1, to give an estimate of the surface tension.13 The basic steps in applying Tikhonov regularization to compute d2rC/dz2 and drC/dz will be explained in the next section and all the required equations derived. The method of generalized cross validation (GCV) used to guide the choice of the regularization parameter, λ, will also be described.18 This is followed by a series of examples that demonstrate the performance of this route to surface tension. Some of the problems encountered will be highlighted, and possible ways of overcoming them will be discussed. 2. Governing Equations Integral Equation for Second Derivative. The direct way of obtaining the first and second derivatives of a set of measured droplet radius, rM(z), data is to fit local spline curves through these data points and then differentiating these fitted curves analytically to yield the required derivatives.13 Since locally fitted curves do not generally have continuous derivatives, this way of computing the derivatives may not lead to reliable results, particularly at the level of the second derivative. The starting point of the present investigation is the following exact equation that relates droplet radius at any height, z, to the second derivative
rC(z) )
z ∫z′)z
0
(z - z′)f(z′) dz′ + r0 + (z - z0)g0 (3)
In this equation and in subsequent derivations, f(z) ≡ d2rC/dz2 and g(z) ≡ drC/dz. Equation 3 is a two-term Taylor series expansion of the droplet radius rC(z) about an arbitrary reference point, z0. The first term on the RHS is the remainder term of the expansion written in integral form and r0 ≡ rC(z0) and g0 ≡ g(z0).19 Equation 3 is in the form of an integral equation of the first kind and is known to be ill-posed.17 It will be solved simultaneously, by Tikhonov regularization, for the second derivative f(z) and the unknown constants r0 and g0.17 Discretized Equation. In the derivation of the equations of Tikhonov regularization, the following notations for experimental data and discretized variables will be adopted. The locations where the droplet radius is measured will be denoted by the vector zM ) (zM 1 ) z0, M M zM 2 ,..., zi , ..., zND) and the droplet radius at these points M M M by rM ) (rM 1 ) r0, r2 , ..., ri , ..., rND). The number of data point, ND, typically varies from as low as 30 to as high as M - zM well over 100. The height of the pendant droplet zN 1 D is discretized into NK uniformly spaced points zC ) (zC1 ) C M M C C M zM 1 , z2 , ..., zj , ..., zNK ) zND) at ∆ ) (zND - z1 )/(NK - 1) distance apart. The values of the unknown second derivative f(z) at these discretization points will be denoted (18) Wahba, G. Spline Models for Observational Data; SIAM: Philadelphia, 1990. (19) Burkill, J. C. A First Course in Mathematical Analysis; Cambridge University Press: Cambridge, 1978.
Tikhonov Regularization for Tensiometry Data
Langmuir, Vol. 21, No. 24, 2005 11243
by the vector f ) (f1, f2, f3, ..., fNK). To ensure accurate representation of f(z), the number of discretization points, NK, is usually set to be much larger than ND, typically NK ) 301-601. In terms of the discretized variables, eq 3 becomes NK
rCi )
∑ j)1
Bijfj + r0 + (zM i - z0)g0 for i ) 1, 2, ..., ND (4a)
or in matrix notation
rC ) Bf + 1r0 + (zM - 1z0)g0
(4b)
B is a ND × NK matrix of known numerical coefficients that arise from the approximation of the integral in eq 3 by numerical quadrature such as the trapezoidal or the Simpson’s rule. 1 is a column vector of 1. The unknowns f1, f2, f3, ..., fNK and r0, g0 are required to minimize ND
(i) S1 )
C 2 (rM ∑ i - ri ) i)1
and NK-1
(ii) S2 )
∑ j)2
(5a)
( ) d2f
2
dz2
j
(5b)
Condition (i) ensures that the computed profile approximates the measured profile closely, and condition (ii) ensures that the second derivative of the radius does not show excessive physically unreal fluctuations.17 Tikhonov Regularization. In Tikhonov regularization, instead of minimizing (i) and (ii) separately, a linear combination, R ) S1 + λS2, is minimized.17 λ is the weighting or regularization parameter that balances these two requirements. A large λ favors (ii) the smoothness condition, while a small λ favors (i) the accuracy condition. The value of λ clearly depends on the noise level in the measured profile data, the number of data points, ND, and that of discretization points, NK; its numerical value therefore does not have any direct physical significance. Methods such as the Morozov Principle,17 the L-curve method20 and the method of GCV18 can be used to guide the choice of this parameter. In this investigation, the choice of λ will be based on GCV. For any λ, the unknowns f1, f2, f3, ..., fNK, r0, and g0 that minimize R are given by17
(
f ′ ) B′TB′ +
)
λ T β β ∆4
-1
B′TrM
(6)
[
0 0 0 0
0 0 0 0
]
(7)
is a modified tri-diagonal matrix arising from standard finite difference approximation of the d2f(z)/dz2 at the uniformly distributed discretization points. The two extra (20) Hansen, P. C. SIAM Rev. 1992, 34, 561.
rC ) B′f ′
(8)
eq 8 is used in the computation of the GCV function below. Equation 6 is the linear algebraic equation that converts a set of droplet profile data, rM, into its second derivative, f(z), described by f. Since f(z) is known at a large number of closely and uniformly spaced points, it can be integrated successively and with relative ease to give the first derivative, g(z), and a back-calculated version, rC(z), of the droplet radius. In the integration operations, the g0 and r0 given by eq 6 are used as boundary conditions. All the integration will be performed using commercial software independent of the computer code developed to solve for f(z). Thus, the comparison of the back-calculated rC(z) with the measured data, rM, serves as an independent check of the reliability of the derivatives given by eq 6. Such comparisons will be performed in all the examples described below. It should be noted that the backcalculated profile obtained by integrating f(z) is independent of the Laplace-Young equation and does not involve the surface tension. In all the droplet profile data to be analyzed, the solution obtained this way will also be compared against that generated by the solution of the Laplace-Young equation as a direct problem using the surface tension deduced from the profile data. Generalized Cross Validation. GCV used to guide the choice of λ is based on the “leaving-out-one” principle.18 In principle, the computation described by eq 6 can be repeated ND times, each time leaving out one data point. The optimum λ is taken to be the value that minimizes the sum of squares, V(λ), of the difference between the predicted value and the actual value of each of the left-out data point. It can be shown that, in the GCV implementation of the “leaving-out-one” principle, V(λ), is given by18
V(λ) )
(rM - rC)T(rM - rC)/ND (1 - Tr[E]/ND)2
(9)
Tr[E] denotes the trace of the square matrix E, known as the influence matrix, defined by18
(
E ) B′ B′TB′ +
)
λ T β β ∆4
-1
(B′)T
(10)
Equations 6 and 8, together with the definition of B′, β, and E, allow V(λ) to be evaluated and plotted against λ/∆.4 Minimization of V(λ) is used to locate the optimum λ. 3. Data and Results
where f ′ ) (f1, f2, f3, ..., fNK, r0, g0) is used to denote all the unknowns and B′ is the matrix B with the column vectors 1, zM - 1z0 added to it to reflect the incorporation of r0, and g0 into f ′. The matrix β defined by
1 -2 1 1 -2 1 β) · · · 1 -2 1
columns of 0 in eq 7 are again the consequence of the incorporation of r0 and g0, which do not feature in the smoothness condition, into f ′. In terms of f′ and B′
Equations 9 and 6 will now be used to compute the derivatives f(z), g(z), and rC(z) of a number of sets of droplet profile data. These computed functions are then converted into curvatures R1 and R2 and used in the computation of surface tension, γ. The ways γ is computed will be described with reference to a specific example. The relevant experimental conditions of the droplet data sets investigated can be found in the original papers and will be omitted here. Simulated Water Droplet. The discrete points in Figure 1a are the simulated profile of water (∆F ) 998.2 kg m-3, γ ) 72.79 mN m-1) presented in the form of rM as a function of the vertical position, zM, in the droplet. These data were generated by solving the direct problem associated with eq 1. There are 68 points in this data set
11244
Langmuir, Vol. 21, No. 24, 2005
Yeow et al.
Figure 1. Simulated water droplet. (a) Droplet Profile. 2, simulated profile; 9, excluded point; dark continuous curve, backcalculated profile from eq 6; light continuous curve, solution of eq 1 with γ ) 72.99 mN m-1. (b) Second derivative. (c) First derivative. (d) GCV plot. (e) Curvatures. Dark continuous curves, curvature given by eq 6. Light straight line, best-fit eq 11. (f) Radius ratio. Dark continuous curve, 2H × rC from eq 6. Light continuous curve, best-fit eq 12. (g) Curvatures with 1.5% random noise. Dark continuous curves, curvature given by eq 6. Light straight line, best-fit eq 11. (h) Radius ratio the with 1.5% random noise. Dark continuous curve, 2H × rC from eq 6. Light continuous curve, best-fit eq 12.
inclusive of the apex of the droplet, shown as a filled square, located at the origin of the coordinate system. At the apex, drC/dz f +∞ while d2rC/dz2 f -∞, and all the higher derivatives, alternating between (∞, are also singular at this point. To avoid the computational difficulties caused by this singularity, the apex will be excluded in this and also in all the subsequent profiles to be analyzed. For profiles where the gradient, drC/dz, close to the apex is large, additional data points will also be excluded. Excluded data points will be indicated by filled squares. The second derivative, f(z), and the first derivative, g(z), of the simulated profile given by eqs 9 and 6 are shown in parts b and c of Figure 1, respectively. The droplet radius, rC(z), obtained by further integration of g(z) is shown as a dark continuous curve in Figure 1a. Even though there is no explicitly imposed noise in the simulated data, noise can clearly be seen in these derivatives, particularly in f(z). Despite the noise, there is a very satisfactory agreement between the back-calculated profile and the simulated data, confirming the general reliability
of the derivatives generated by eq 6. This is also observed for all the droplet profiles investigated. The regularization parameter used to arrive at these results is λ/∆4 ≈ 9.40 × 10-6 (dimensionless) given by the minimization of the GCV function shown in Figure 1d. The profile data have been converted into dimensionless form prior to the Tikhonov regularization computation, so the GCV function plot in Figure 1d is in a dimensionless form. As already mentioned, no physical significance should be read into the value of λ used. The behavior of the function in Figure 1d is representative of the GCV functions for all the droplet profiles investigated. Unless they feature explicitly in discussion, the corresponding GCV function in the subsequent examples and the associated f(z) and g(z) will not be presented. The computed f(z), g(z), and rC(z) are substituted into eq 2 to give the principal radii of curvature, R1(z) and R2(z), of the droplet. Two methods will be used to deduce the surface tension, γ, from these radii of curvature. In the first method, these are plotted as 1/R1(z), 1/R2(z), and
Tikhonov Regularization for Tensiometry Data
Langmuir, Vol. 21, No. 24, 2005 11245
Table 1. Regression Parameters of Equation 12 and Surface Tension figure
R (mm-1)
β (mm-2)
γ (mN m-1)
1f 2c 3c 4d 5g 5h 5i
1.73505 1.27323 1.70943 1.7846 1.73318 1.9350 2.82025
-0.134168 -0.132267 -0.319228 -0.330876 -0.134073 -0.132473 -0.132598
72.99 73.95 22.39 23.40 73.02 73.90 73.84
2H(z) in Figure 1e. According to eq 1b, 2H(z) varies linearly with of z, i.e.
2H(z) ) R + βz
(11)
where R ) 2/RApex and β ) -∆Fg/γ. It can be seen that, apart form the neighborhood of the apex, the computed 2H(z) does indeed vary linearly with z. Close to the apex, because of the ill-posed nature of integral equation of the first kind, 1/R1(z) and 1/R2(z) show significant fluctuations, and consequently, 2H(z) deviated from linearity. As these fluctuations are physically unreal, they will be ignored.13 The least squares best-fit of eq 11 through the linear portion of the computed 2H(z) data, for z between 0.8 mm e z e 2.35 mm, is shown as a light straight line in Figure 1e. The 2H(z) results for z < 0.8 mm have been ignored. Following de Ramos et al.,13 the slope of the best-fit line is used to compute the surface tension. For the straight line in Figure 1e, β ) -0.133194 mm-2, giving a surface tension γ ) 73.52 mN m-1. Similar to the data in Figure 1a, all the other droplet profiles analyzed in this investigation also showed significant fluctuations in the 2H(z) versus z plots. For some of these experimental data, the fluctuations are not confined to the neighborhood of the apex. Such fluctuations have been reported by de Ramos et al.13 Least-squares straight lines through the computed 2H(z), even for a restricted range of z, therefore may not provide a reliable estimate of γ. To overcome this problem, a second method that makes use of not only the 2H(z) given by eq 6 but also of the rC(z) will be explored. In this new method, eq 1b is multiplied by rC(z), and this leads to
2H(z) × rC(z) )
( )
( )
2 ∆Fg rC(z) z × rC(z) ) RApex γ RrC(z) + βz × rC(z) (12)
R and β are the same two parameters in eq 11. The product H(z) × rC(z) can be viewed as the ratio rC(z):mean radius of curvature, 1/H(z). From the expression for H(z) in eq 2, this ratio can be expected to be better behaved numerically than the computed H(z). In particular, as z becomes small, this ratio is not expected to exhibit the kind of large amplitude fluctuations exhibited by H(z). Since H(zj) and rC(zj) are known for j ) 1 to NK, R and β can be determined by standard least-squares regression between the computed 2H(z) × rC(z) from eq 6 and the simple expression on the second line of eq 12. For the simulated water droplet, the computed 2H(z) × rC(z) and the best-fit eq 12 are shown in Figure 1f as a continuous dark curve and a lighter curve, respectively. The numerical values of the regression parameters R and β are summarized in Table 1. This regression covers the entire range of z. Apart from some local fluctuations, 2H(z) × rC(z) and the best-fit eq 12 are in satisfactory agreement. On the basis of this new regression, γ ) -∆Fg/β ) 72.99 mN m-1, which is larger than the value used in the generation of
the profile by approximately 0.3% and is closer to the expected value than that given by restricted regression of eq 11. It should be mentioned that this level of accuracy is routinely attainable with eq 12 and not a selected best outcome. The reliability of this γ can be verified by using it, together with the r0 and g0 given by eq 6, to solve the direct problem of eq 1. The resulting droplet profile is shown as a light continuous curve in Figure 1a and is in satisfactory agreement with the other two profiles shown on this plot. To simulate the effects of experimental noise, small random noise is added to the z coordinates of the synthetic data in Figure 1a and the procedures for obtaining surface tension described above are repeated. Parts g and h of Figure 1 are the analogues of parts e and f of Figure 1, respectively. The level of the added random noise is a relatively small (1.5%. This noise showed up as amplified fluctuations in the first and second derivatives of R(z). While these derivatives (not shown) remain relatively smooth, when they are used to compute H(z) and the ratio 2H(z) × rC(z) the noise are further amplified and these can be see in Figure 1g and h. Even though the regularization parameter has been adjusted so that the average deviation between the noisy data and the back-calculated profile based on Tikhonov regularization is consistent with the level of the added noise, the fluctuations in H(z) and 2H(z) × rC(z) become more widespread. The resulting surface tension from the best-fit straight line in Figure 1g is 61.01 mN m-1, and the corresponding result from the best-fit radius ratio curve in Figure 1h is 74.81 mN m-1. Attempts to increase the level of the added random noise to 5% led to surface tensions that are not in acceptable agreement of the expected value. In a recent paper, Lubansky et al.21 showed Tikhonov regularization can cope with random noise and is able to generate reliable first and second derivatives for data with noise levels as high as (20%. Thus, the inability of Tikhonov regularization to obtain acceptable surface tension at similar levels of added noise is likely to be caused by the further amplification of noise in the computation of R1(z), R2(z), and H(z). This explains the extreme sensitivity of the pendant drop method to experimental noise, such as those generated by lighting and the digitizing process and also by the distortions introduced by the optical system employed. Water Droplet in Saturated Nitrogen. To test their method of processing pendant droplet data, Jennings and Pallas15 tabulated the profile of a water droplet at 22 °C suspended in a nitrogen atmosphere saturated with water vapor. Their data are replotted as discrete points in Figure 2a. The LHS and RHS profiles, reported separately by these authors, have been combined and treated as a single set. There are only 28 points in this combined set. This is an extremely small data set compared to the large data sets routinely extracted from CCD images of pendant droplets by the current generation of edge detection software. Since there are very perceptible differences between the LHS and RHS profiles, additional noise will be introduced by treating them as a single set. All these factors should be taken into consideration when assessing the reliability of the surface tension based on the profile data in Figure 2a. The apex is again omitted from the data processing. As in the previous example, eqs 9 and 6 are used to compute f(z) and g(z). The rC(z) back-calculated from these derivatives is shown as a dark continuous curve in Figure (21) Lubansky, A. S.; Yeow, Y. L.; Leong, Y.-K.; Wickramasinghe, S. R.; Han, B. AIChE J. 2005, in press.
11246
Langmuir, Vol. 21, No. 24, 2005
Figure 2. Experimental water droplet. (a) Droplet Profile. 2, experimental data;15 9, excluded point, dark continuous curve, back-calculated profile from eq 6; light continuous curve, solution of eq 1 with γ ) 73.95 mN m-1. (b) Curvatures. Dark continuous curves, curvature given by eq 6. Light straight line, best-fit eq 11. (c) Radius ratio. Dark continuous curve, 2H × rC from eq 6. Light continuous curve, best-fit eq 12.
2a. These computed functions are also used to construct the 2H(z) versus z plot shown in Figure 2b and the 2H(z) × rC(z) versus z plot shown in Figure 2c. It is noted that the fluctuations in 2H(z) × rC(z) are reduced in magnitude compared to that exhibited by 2H(z). The best-fit eqs 11 and 12 are shown as lighter curves on the respective figures. In this and in subsequent regression of eqs 11 and 12, because of the generally observed slightly irregular behavior of f(z) and hence of H(z) in the neighborhood of M zN , the regression computation will exclude this neighD borhood. The excluded region is arbitrarily fixed at approximately 10% of the droplet height, which correspond roughly to the start of the observed irregularities. Another justification for excluding this region is to isolate the perturbation caused by the attachment point of the droplet to the capillary wall. In Figure 2b, because of the spread of the fluctuations in the computed data, it is no longer practical to completely isolate these fluctuations, and consequently, the best-fit straight line resulted in a γ ()87.28 mN m-1) that is significantly larger than the expected value. In Figure 2c, the regression parameters for eq 12 are shown in Table 1. Despite the fluctuations, this regression gave γ ) 73.95 mN m-1. This figure should be compared against the values 71.79 mN m-1 e γ e 72.40 mN m-1 obtained by Jennings and Pallas.15 In view of the small number of data points and the additional noise introduced by the combination of the LHS and RHS profiles into a single set, the result of γ ) 73.95 mN m-1 can be regarded as very satisfactory. This result again represents what can be routinely achieved and not the best performance under ideal conditions. This remark applies to all the results in this paper. The back-calculated droplet profile based on γ ) 73.95 mN m-1 and the r0 and g0 given by eq 6, is shown as a light continuous curve in Figure 2a. It is suspected that the relatively poor fit of the regression curve in some part of Figure 2c is responsible for the deviation of the back-calculated profile curve from the experimental data for large z.
Yeow et al.
Figure 3. Experimental n-decane droplet. (a) Droplet Profile. 2, experimental;13 dark continuous curve, back-calculated profile from eq 6; light continuous curve, solution of eq 1 with γ ) 22.39 mN m-1. (b) Curvatures. Dark continuous curves, curvature given by eq 6. Light straight line, best-fit eq 11. (c) Radius ratio. Dark continuous curve, 2H × rC from eq 6. Light continuous curve, best-fit eq 12.
Figure 4. Ethanol droplet. (a) Droplet Profile. 2, experimental profile;14 9, excluded point; dark continuous curve, backcalculated profile from eq 6; light continuous curve, solution of eq 1 with γ ) 23.40 mN m-1. (b) Second derivative. (c) First derivative. (d) Radius ratio. Dark continuous curve, 2H × rC from eq 6. Light continuous curve, best-fit eq 12.
n-Decane Droplet. de Ramos et al.13 used locally fitted splines to evaluate the derivatives and the mean curvature, H(z) of pendant droplets. They then obtained the surface tension from the slope of the best-fit eq 11 through the 2H(z) versus z data. These authors used the profile of an n-decane droplet to demonstrate the performance of their technique. This profile, reconstructed by automatic scanning and digitizing from their Figure 4, is reproduced as discrete points in Figure 3a. As scanned, this profile has 587 data points. Only every fifth of the scanned points are shown in Figure 3a and used in the present analysis.
Tikhonov Regularization for Tensiometry Data
In should be noted that this data set does not include data points in the neighborhood of the apex, and therefore, exclusion of data points for small z in the regression computation of eqs 11 and 12 may not be justified. The back-calculated profile based on the derivatives given by eqs 9 and 6 is shown as a dark continuous curve in Figure 3a. The 2H(z) versus z data constructed from the results of these two equations is shown as a dark continuous curve in Figure 3b. There is significant departure from linearity, particularly at the two ends of this plot even though the apex region was not included in the original data set. The departure from linearity at the lower end is clearly a consequence of the ill-posed nature of evaluating the derivatives of experimental data. de Ramos et al.13 reported similar departure from linearity for their plot constructed from locally fitted splines. The best-fit eq 11 to this plot is shown as a straight line in Figure 3b. In the regression computation for this line, the lower 10% and the top 10% of the data have been arbitrarily excluded.13 The surface tension deduced from the slope of this best-fit straight line is γ ) 21.59 mN m-1. The computed 2H(z) × rC(z) versus z plot and the corresponding best-fit eq 12 are shown in Figure 3c. The parameters of this regression are tabulated in Table 1. The surface tension deduced from the best-fit β is γ ) 22.39 mN m-1. In this regression, only the top 10% of the data have been excluded from the regression computation, and this corresponds to the clearly unexpected change in trend of the data in Figure 3b and c. As already mentioned, since the data set does not include the apex of the droplet, exclusion of the lower region is not justified and was indeed found to be unnecessary in this case. If no data have been excluded at all in the regression of eq 12, the corresponding γ ) 23.37 mN m-1. Both of these γ are in acceptable agreement with the 22.2 mN m-1 < γ < 23.5 mN m-1 reported by de Ramos et al.13 As in the previous examples, the back-calculated droplet profile given by γ ) 22.39 mN m-1 is shown as a light continuous curve in Figure 3a. The deviation of the regression curve from the computed data in Figure 3c is again responsible for the deviation observed in this backcalculated droplet profile. Not withstanding these significant deviations, the regression parameter β given by eq 12 has again extracted a reasonably accurate surface tension from the droplet profile data. Ethanol Droplet. Dingle et al. applied Galerkin FEM to solve the Laplace-Young equation as a differential equation for rC(z).14 They then obtained the surface tension by performing a nonlinear regression computation between the droplet profiles given by FEM and by experimental measurements. The experimental data of an ethanol droplet used by these authors to demonstrate their procedure are shown as discrete points in Figure 4a. The data points from their LHS and RHS profiles have again been combined into a single set and processed together. Apart from the apex, two additional data points each of the LHS and RHS profiles have been excluded in the data processing. The excluded points are shown as filled squares in Figure 4a. As the LHS and RHS profiles are identical over these points, there are thus only three distinct filled squares in this figure. These points are collinear and have a relatively large gradient. If these points are not excluded, GCV will home on to a very small λ, leading to smallamplitude but high-frequency fluctuations in the second derivative f(z). These fluctuations become grossly amplified in the curvature H(z). In this and in the next example, obtaining γ via linear regression of H(z) versus z will not be attempted.
Langmuir, Vol. 21, No. 24, 2005 11247
The f(z), g(z), and 2H(z) × rC(z) versus z plots of the ethanol droplet given by eqs 9 and 6 are shown in Figure 4b-d. There are perceptible differences between the LHS and RHS data points of the ethanol droplet for z > 2.5 mm (approx), and these show up clearly as changes in the trend in f(z) and in g(z) in this region which appeared as large fluctuations in the 2H(z) × rC(z) versus z plot (see Figure 4d). In the regression computation of eq 12, this region of Figure 4d will therefore be excluded. The bestfit parameters are summarized in Table 1, and these resulted in γ ) 23.40 mN m-1. This should be compared against the average value of 22.60 ( 0.25 mN m-1, based on the results of five separate droplet profiles, reported by Dingle et al.14 Effects of Droplet Volume. It is well known that the reliability of pendant droplet tensiometry decreases with droplet volume. Lin et al.5 considered this problem, and as part of their investigation, they reported the droplet profiles of water with volume V ) 7.88, 5.34, and 1.58 mm3. Subsets of these profiles are plotted as discrete points in parts a, b, and c of Figure 5, respectively. The original data sets contain a large number of points. The three subsets reproduced here are based on every 5th, 3rd, and 2nd point, respectively, of the RHS profiles reported by these authors. The original rM are reported in pixels, and these have been converted into millimeters using the average conversion factor of 86.67 pixels mm-1. This linear conversion factor is based on the volumetric conversion factor deduced from the reported physical droplet volumes in cubic millimeters and the corresponding volume in cubic pixels obtained by numerical integration of the droplet profiles tabulated in pixels. The resulting conversion factors differed slightly, by about 1-2%, for the three different droplet volumes. The average linear conversion factor has also been checked against the specified OD of the capillaries shown on the profile plots in Lin et al.5 Plots of the second derivative, f(z), given by eqs 9 and 6, for the three water droplets are shown in Figure 5d-f. Fluctuations are observed in all three plots, but it is clear that they all approach a small constant value as z increases. The magnitude of the second-order derivatives and the frequency of the small fluctuations appear to increase without limit as the droplet size is decreased. For the V ) 1.58 mm3 droplet, the data points close to the apex, shown as filled squares in Figure 5c, have been excluded. This reduces the magnitude of the derivatives of the data close to the apex and is responsible for reducing the frequency of the fluctuations in Figure 5f. It is also noticed that, for all three droplets, close to the attachment point of the droplets, f(z) departed from the general trend of approaching a constant value. This is particularly clear in Figure 5f for the V ) 1.58 mm3 droplet. This change in trend is likely to be a consequence of the perturbations caused by the droplet attachment point or the end effects of the Tikhonov regularization computation. The curvatures of the droplets in this neighborhood will be excluded in the subsequent regression computation. The backcalculated droplet profiles based on the f(z) and g(z) (not shown), plotted as dark continuous curves in Figure 5a-c, confirm the general reliability of these derivatives. The computed 2H(z) × rC(z) versus z plots for the three water droplets are shown in Figure 5g-i. Amplified fluctuations are again observed in these plots. Because of the excluded points for the V ) 1.58 mm3 droplet, the fluctuations in Figure 5i are greatly reduced (see Figure 6). Similar reduction can be achieved for the larger droplets, but this has not been implemented in Figure 5g and h. The best-fit eq 12 to the 2H(z) × rC(z) versus z data are shown as lighter curves on the same figures, and the
11248
Langmuir, Vol. 21, No. 24, 2005
Yeow et al.
Figure 5. Water droplets of different volumes. (a), (b), and (c) Profiles for V ) 7.88, 5.34, and 1.58 mm3, respectively. 2, experimental profile;5 9, excluded point; dark continuous curve, back-calculated profile from eq 6; light continuous curve, solution of eq 1 with γ ) 73.02, 73.90, and 73.84 mN m-1, respectively. (d), (e), and (f) Second derivatives of (a), (b), and (c) from eq 6. (g), (h), and (i) Radius ratio of (a), (b), and (c). Dark continuous curve, 2H × rC from eq 6. Light continuous curves, best-fit eq 12.
Figure 6. Water droplet V ) 1.58 mm3. (a) Second derivative given by eqs 6 and 9 with only the apex excluded. (b) Radius ratio. Dark continuous curve, 2H × rC from eq 6. Light continuous curve, best-fit eq 12.
regression parameters are summarized in Table 1. As already mentioned, in these regression computations, the 2H(z) × rC(z) data in the neighborhood where f(z) deviated from the general expected trend have been excluded. The surface tensions based on the three regression parameters, β, are 73.02, 73.90 , and 73.84 mN m-1 for the large, medium, and small droplets, respectively. Despite the physically unreal fluctuations in the computed 2H(z) × rC(z), all these values are reasonably close to one another
and to the expected value of 72.67 mN m-1. The corresponding values read off from Figure 4 in Lin et al.5 are 72, 70, and 63 mN m-1. These authors attributed the effects of droplet volume on surface tension to the progressive deterioration of the regression computation between computed and measured droplet profiles, particularly in the neighborhoods of the attachment points, as the droplet volume is reduced. For the three values of γ just obtained together with the associated r0 and g0, the back-calculated profiles given by the solution of the direct problem of eq 1 are shown as light continuous curves in Figure 5a-c. For all three droplets, there is good agreement between the backcalculated profile based on the surface tension, the original data, and that given by integrating the derivatives g(z) and f(z). 4. Discussion In this investigation, Tikhonov regularization is used to obtain the first and second derivatives of pendant droplet profiles and hence the mean curvature as a function of height. The regression computation between the radius ratio 2H(z) × rC(z) and R rC(z) + βz × rC(z) is then used to obtain β, from which the surface tension is calculated.
Tikhonov Regularization for Tensiometry Data
For most of the droplet profiles analyzed, this route to surface tension proved to be reliable and remains applicable when linear regression between 2H(z) and R + βz, because of amplified noise in the mean curvature, becomes unreliable. The success of this method can be attributed partially to the fact that Tikhonov regularization allowed for the ill-posed nature of the profile differentiation step. Suppression of the noise in H(z) close to the apex by the product H(z) × rC(z) is another contributing factor. The key equation of Tikhonov regularization is eq 6. Derivation of this equation may appear lengthy, but the physical motivation behind it is direct. More importantly, the equation itself is also relatively straightforward to apply and does not involve iterative computation. All the numerical operations in setting up the matrixes and their subsequent modifications and inversion can be performed using standard commercial computer software. As already mentioned, most existing methods of processing pendant droplet data involve nonlinear regression computation between the computed drop profile and the measured profile. This requires repeated solution of the Laplace-Young equation, eq 1, as a differential equation for rC(z). For each adjustment of the regression parameter, the integration of this differential system has to be repeated. Some regression computations between rC(z) and rM(z) involve iterative adjustments of more than one parameter such as the surface tension, γ, and the mean radius of curvature, RApex, at the apex. The convergence of this nonlinear regression process depends on a proper initial guess of these parameters. Thus, this route to surface tension can become computationally intensive. In contrast, the present way of computing γ treats the Laplace-Young equation as an algebraic equation in the form shown in eq 12 and does not require the integration of a differential equation. The regression computation of eq 12 for the parameters R and β is essentially a problem in linear regression and therefore also does not require iterative computation. The only time-consuming computation is the evaluation of the GCV function, V(λ). This again does not involve iterative computation and therefore does not suffer from convergence difficulty. Thus, apart from its reliability, the route to surface tension via eq 12 is also computationally simpler. It was found that, for the regression computation of eq 12 to work properly, it is necessary to exclude some of the 2H(z) × rC(z) data. While it is possible to put forward physical justifications for this exclusion, at present, the ultimate decision of how large a region (or regions) has to be excluded requires case-by-case consideration. For general adoption of the present procedure for processing pendant droplet data, it is clearly desirable to have an automated procedure for making this decision. Such a procedure needs to be based on the analysis of a large number of case studies and therefore awaits the accumulated practical experience from different users. Similar to most existing methods, the present method also encountered difficulties when the droplet volume becomes small. For the present method, these show up as large amplitude fluctuations in the computed 2H(z) × rC(z) versus z data. This is typified by the V ) 1.58 mm3 profile in Figure 5c. If the data points close to the apex have not been excluded, the resulting second derivative, f(z), shows high-frequency fluctuations (see Figure 6a). These appeared as amplified fluctuations in the 2H(z) × rC(z) shown in Figure 6b. The corresponding best-fit eq 12 in this case resulted in a surface tension of around 65.97 mN m-1, comparable with that reported by Lin et al.5 for this set of data.
Langmuir, Vol. 21, No. 24, 2005 11249
Figure 7. Fourth derivative results for n-decane. (a) Curvatures. Dark continuous curves, back-calculated curvature from fourth derivative; light straight line, best-fit eq 11. (b) Radius ratio. Dark continuous curve, 2H × rC from fourth derivative; light continuous curve, best-fit eq 12.
From its physical meaning, the product 2H(z) × rC(z) is expected to increase smoothly from a small value close to the apex to a maximum of around 2 and then decreases smoothly with further increase in z. The high-frequency fluctuations observed in most of the 2H(z) × rC(z) versus z plots are not physically real. Suppressing them is therefore unlikely to lead to loss of physical information. Such fluctuations can be reduced by increasing the value of the λ, thereby giving a smoother f(z). It is suggested that instead of being guided by GCV, λ should be increased, by orders of magnitude if necessary, so that the resulting 2H(z) × rC(z) versus z plots exhibit a physically more realistic behavior. In addition, the selected λ should also satisfy the Morozov Principle,17 i.e., the average deviation between the back-calculated droplet profile based on f(z) and the experimental data should be comparable with the error bars of the measurements. Numerical experimentations show that these two physically based requirements often provide a computationally less expensive and yet reliable means of picking λ that in turn leads to a reliable β and hence surface tension. Other methods of overcoming the fluctuations in H(z) and in 2H(z) × rC(z) have also been explored. One possible method is to extend the two-term Taylor series expansion in eq 1 to three or four terms. The resulting analogue of eq 1 then becomes an integral equation for the third derivative or the fourth derivative to be solved by Tikhonov regularization. These higher derivatives can be then integrated and hopefully will give a smoother second and first derivatives to be used in the construction of the 2H(z) versus z plot and the 2H(z) × rC(z) versus z plot. Figure 7a and b shows examples of such plots constructed from the fourth derivative. These plots are for the n-decane profile shown in Figure 3a. They should be compared against the corresponding plots in Figure 3b and c based on the second derivative given by eq 6. Improvements are generally discernible and particularly close to the attachment point. The corresponding best-fit eqs 11 and 12 are shown as a lighter line and curve, respectively, on the corresponding plots. The surface tensions deduced from the two best-fit β are 21.60 mN m-1 from eq 11 and 22.70 mN m-1 from eq 12. Higher derivative computations with other profile data confirm that the improvement observed is general. However, the improvement observed should be weighed against the vastly increased computation cost of solving for the higher derivatives, particularly the effort needed in the construction the GCV plots for these higher derivatives. At present, the additional computation effort required is not considered as commensurate with the improvement achieved. Possible ways of reducing this computing cost is under investigation.
11250
Langmuir, Vol. 21, No. 24, 2005
5. Conclusions Tikhonov regularization in which GCV is used to guide the selection of the regularization parameter provides a reliable means of obtaining the radii of curvature of pendant droplets. For large droplets, least-squares regression based on the computed 2H(z) × rC(z) provides a reliable route to surface tension. For small droplets, it may be better to adjust the regularization parameter manually to get a well-behaved 2H(z) × rC(z) instead of relying on GCV. An alternative is to exclude a finite
Yeow et al.
number of data points in the neighborhood of the apex of a small droplet where the derivatives are close to singular prior to performing the Tikhonov regularization-GCV computation. Acknowledgment. The authors gratefully acknowledge the pendant droplet profile data provided by Prof M. T. Harris and by Prof S.-Y. Lin. LA051650P
