- Email: [email protected]
Three-dimensional menisci: Numerical simulation by finite elements
Three-Dimensional Menisci: Numerical Simulation by Finite Elements F. M. ORR, JR., R. A. B R O W N , AND L. E. S C R I V E N
Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, Minnesota 55455 Received May 22, 1976; accepted July 12, 1976 Menisci with complicated shapes occur in wetting and drying of powders and fabrics, coating of surfaces, fluid displacement in porous media, and other technologies. To predict meniscus shape and location requires the solution of the Laplace-Young equation, which is nonlinear and elliptic and depends on two space variables. The finite-element method has particular advantages that are further demonstrated here. Its recent application to menisci between cylinders is extended to three-dimensional menisci of unbounded extent or bounded by free contact lines (free boundary problems). Specimen solutions are displayed, including multiple solutions for the meniscus in a simple threedimensional array of spheres, an elementary type of porous medium. INTRODUCTION Menisci that possess neither translational nor rotational symmetry are three-dimensional. Although governed by the well-known LaplaceYoung equation, their full description requires the numerical solution of this nonlinear, elliptic equation. In the antecedent I (1), the finiteelement method was shown to be particularly effective because it can easily accommodate the irregular boundaries and nonlinear boundary conditions that generally accompany threedimensional menisci. Although this method, as developed in (1), can be used to solve the Laplace-Young equation for more complicated meniscus shapes than any other available method, there are important meniscus problems that cannot be solved by the method as there constituted. There are three sorts of difficulties that arise: (1) menisci of unbounded extent, i.e., domains of infinite expanse; (2) free contact lines, i.e., domain boundaries the location of which is part of the problem; and (3) menisci which do not have a single-valued representation in the customary Cartesian coordinate
system or an equivalent system. I n this paper we extend the finite-element method in order to deal with the first two difficulties and by way of illustration we solve several problems of scientific or practical interest. DIFFERENTIAL EQUATION AND BOUNDARY CONDITIONS In the notation used in I, the LaplaceYoung equation is 2(H -- t!o)g = g(Z -- Zo)Ap
[11
and the mean curvature of the surface is proportional to the surface divergence of the unit vector N, which is everywhere normal to the surface: 2H = - V s . N . When the surface can be represented by a single-valued elevation Z above a datum plane, it happens that V~-N = VII.N, where VII is the gradient operator in the datum plane. If that plane is parameterized with Cartesian coordinates (x, y), VII = i o / o X + jo/Oy; if polar coordinates (r, O) are used, VII = era~Or + eor-lO/O0. 137
C o p y r i g h t ~ 1977 by A c a d e m i c Press, Inc. All rights of reproduce.ion in a n y form reserved,
Journal oJ Colh)id and Interface Science, Vol. 60, No. l, J u n e 1, 1977 I S S N 0021-9797
138
ORR, BROWN AND SCRIVEN the expansion is [-cf. (1)]
In either case N = (VnZ)[-1 + (VliZ)2]-½.
[-2-]
3
h = E hi~, ee) (f, ,)
When the effect of gravity is significant, as in the cases considered here, it is convenient to scale all lengths with the capillary length a--(a/gAp)½. This is done in I and the Laplace-Young equation is shown to become VII" N + h = 0,
N = (k - ih~ - j h , ) / ( 1 + h~~ +
x ~-----,
h,~)½,
y Z -- Z0 •-------, h-~ +2H0a.[3]
Where the meniscus intersects a solid surface three types of boundary conditions are possible: a Dirichlet condition in which the location of the contact line is specified; a Neumann condition in which the contact angle between the fluid interface and the solid surface is specified; or a Robin condition in which location and angle are related. In the simplest instance the Neumann or contact angle boundary condition is N . n = cos 0o
[-4]
where n is the normal to the solid surface and 0¢ is a characteristic contact angle, sometimes held to be an equilibrium property, exhibited by the solid-liquid-fluid trio. Owing to the nonlinearity of N the contact angle condition is fairly difficult to handle. We adopt it for the sample problems solved here.
(~ and V refer to dimensionless Cartesian or cylindrical polar coordinates). We have experimented with linear trial functions (1) and with reduced Hermite cubic trial functions (3-5). We found that the increased accuracy of this cubic interpolation does not justify its added cost in computation time; that is, comparable accuracy can be obtained at lower cost with linear elements on a more refined mesh. Therefore, the results presented here were all obtained with linear trial functions. It is possible, however, that other higher order trial functions would have advantages. The Galerkin form of the finite-element method selects the approximate solution which makes the residual in the differential equation orthogonal to each of the trial functions
fa
4Ji(V.N+h)dA = 0 , i = 1 , 2 , . . . , Q , [6]
where A is the area of the domain and Q is the number of trial functions. With linear elements Q equals the number of nodes. An integration by parts followed by an application of the divergence theorem generates the nonlinear algebraic equations used to determine the expansion coefficients in Eq. [-5]: f
fo nb'dpiNds + I A
NUMERICAL METHOD The finite-element method used here follows closely that derived in I. Details are given in (2). The elevation of the surface is sought as an expansion in terms of locally defined polynomial trial functions. The problem domain is divided into triangular subdomains, or elements, and within each element a small number of nonzero trial functions is used to interpolate the elevations at the triangle vertices (nodes). Thus in a particular element
[-53
Jx
(4~ih - V$i'N)dA = 0, i = 1,2,...,Q,
[73
where nb is the normal to the domain boundary OA. The boundary conditions enter the approximate equations through the boundary integral. Equations [-7"]for the expansion coefficients hi are a nonlinear algebraic set which has the form
Journal of Colloid and Interface Science, Vol. 60, No. 1, June 1, 1977
F(h) = O.
E83
THREE-DIMENSIONAL MENISCI This set can be solved by successive approximations, which was done in I, or'by'means of a Newton-Raphson iteration which estimates the correction to the vector of expansion coefficients at the (k 4- 1)st iteration as
139
INTERIOR DOMAIN
DOMAIN y
h~k+l~ = h ( ~ -
F
OFi-1-t
L
OhsJ
l--/
F(h(k)) •
I-9-1
Whichever iterative scheme is used, the matrices which must be inverted are quite sparse ; they have a maximum of seven nonzero coefficients on any row. We found the sparse matrix techniques discussed by Gustavson (6) to facilitate greatly the required Gaussian eliminations. INTERFACES
OF UNBOUNDED
Fro. 1. Single cylinder of elliptical cross section in a pool of infinite extent. Sample curvilinear triangular
EXTENT
If a meniscus extends far enough outward from solid bodies it contacts, its equilibrium shape in a gravitational field becomes sensibly flat at some distance from its line of contact with the solid. From a mathematical point of view the interface may then be regarded as unbounded. But to impose a boundary condition at infinity is not possible in a finite set of finite elements. The simplest approximation is to impose the condition at some finite yet large distance, and then to increase the distance until it is sufficiently large that further increase leaves the problem solution unchanged at the desired level of accuracy. The disadvantage of this scheme is that a very large domain may be required, which in turn requires that a very large numerical problem be solved to the desired level of accuracy. The disadvantage can be circumvented by using knowledge of the asymptotic solution for the nearly flat meniscus far from solid bodies. Now lim Z (x, y) = Z0, OZ OZ lim-- = 0 = lim--, ~ Ox r ~ Oy
OA s
El0]
where r = (x 2 + y2)~ measures distance from solid bodies and Z0 is the elevation of a
finite e l e m e n t s d e r i v e d shown. , ~ A/B.
from
polar
coordinates are
perfectly flat meniscus. Thus at large distances r we have (OZ/Ox)2<< 1, (OZ/Oy)2<< 1, Ho = O, Z0 = O, and Eq. [1] becomes ~Z
=
gZAp.
Ell-]
This is Helrnholtz's equation, a well-known linear equation for the solution of which there are standard analytic methods. Its solution satisfying a boundary condition at infinity still contains unknown constants; these can be eliminated by matching the exterior asymptotic solution to an interior approximate solution of the full Eq. E1-]. To demonstrate this approach we use the meniscus around a cylindrical object standing vertically in a pool of unbounded extent. If the cylinder has a circular cross section, the meniscus is axisymmetric, and that problem was solved numerically by Huh and Scriven (7). Otherwise, the meniscus is three-dimensional, and no solution of such a problem on an infinite domain has appeared. A sample domain for a cylinder of elliptical cross section is shown in Fig. 1. Owing to its symmetry the problem can be solved in one quadrant. The domain is divided into an interior portion where the nonlinear, full problem is solved by the finiteelement method, and an exterior portion where
Journal of Colloid and Interface Science, Vol. 60, No. 1, June I, 1977
140
ORR, BROWN-AND-SCRIVEN TABLE I
Comparison of Results Computed for Axisymmetric Menisci around a Cylinder Case
R~
C o n d i t i o n a t r = R~
0 I 2 3 4
-a 3a 1.4a a
H u h and S c r ive n (7) 0 = 7r/2 0 = 7r/2 M a t c h i n g condition M a t c h i n g condition
Height at cylinder
Error (%) (vs Case 0)
Computed contact a ngle (degrees)
0.189 0.306 0.187 0.188 0.188
-100 1.09 0.68 0.53
60 60.8 62.6 61.1 60.7
the linear , asymptotic problem is solved. The two problems are coupled at the matching boundary, a circle of radius r = Rm where the matching condition is that the solution and its first derivatives are continuous. The exterior problem is V2h -- h = 0,
r >
Rln,
Oh
Oh O0
oa rhinb"Nds
lira h = l i r a - - = l i m - - = 0, r~Qo
r~:
Or
formly for r > Rm, yet its r derivative does not converge at r = Rm. The values of Oh/Or at r = Rm, which are needed for the matching condition, are therefore approximated with a forward finite difference just outside the matching boundary. When (OhiO0)2<<1, this series converges rapidly. The data from the asymptotic solution enter the interior problem through the boundary integral of [7-]. Integration is over the solid boundary for which the contact angle condition V4-] is given, over the matching boundary r = Rm on which the normal to the surface is given by the exterior solution, and over the two symmetry boundaries 0 = 0 and 0 = 7r/2. The contributions along the symmetry boundaries vanish because nb.N is always zero there. Thus the boundary contribution to the weighted residual integrals is
,~o
=
fOAs rhinb"Nds
[12]
Oh O0
-0
h = P(O)
at
0 = 0
and
at
r = Rm,
where P(O) is the finite-element expansion for the interior solution evaluated at the matching boundary. Solving the exterior problem with a finite Fourier transform yields a series solution in terms of modified Bessel functions of the second kind:
K2.(r) h(r, O)= ~, a n p , , n=o Kz,(Ro)
+ ] 4~inm"Nds. [14] .IoAm
0 = 7r/2,
On the solid boundary OAs, nb = --ns and on the matching boundary OAm,nm = er. Inserting the contact angle condition and evaluating nn~"N yields
Y fo ~inb"Nds = -- I A
+ fo
~o
cos (2n0)
[-13]
t2/Trl,
n = 0 n =
COS
Oo4s
4)ihrds [15] A,,~[1 + h~ + (ho/r)q~"
Because the square of the slope is supposed to be small at r = Rm, this can be simplified to
where 2/Tr)b
,,tOA s
1, 2, . . . ,
foA 4)inb"Nds -
--
fOAs dpiCOSOcds
and -]- f O A m ~ i h , d s .
Tr[2
p,, =--a.
f
P(O) cos (2nO)dO.
dO
The series and its derivatives converge uniJournal of Colloid and Interface Science,
[-16]
Equations [-7] and [-16] make up the nonlinear algebraic set to be solved for the expan-
Vol. 60, No. 1, J u n e 1, 1977
THREE-DIMENSIONAL MENISCI sion coefficients hi. A three-point Gaussian quadrature formula (8) was used to evaluate area integrals in each element and boundary integrals were evaluated analytically. The Newton algorithm [14"] converged quickly for all cases tested. Typically, three to five iterations were required to reduce the maxim u m adjustment in the h~ to less than 10-4. SPECIMEN MENISCI AROUND CYLINDERS When the ellipse indicated in Fig. 1 has ellipticity ~ = 1 the solid is a circular cylinder of radius B. If the contact angle e~ is uniform all around the cylinder the meniscus extending outward is axisymmetric and the finite-element solution should agree with the results of H u h and Scriven (7); this provides one test of the method adopted here. In addition, the matching strategy used here can be tested by setting the slope equal to zero at the matching radius r = R .... This corresponds to solving the problem of a meniscus in the annulus between r = B and R = R .... with contact angles of 0~ on the inner cylinder and 7r/2 on the outer wall. Table I compares with the more exact solution a number of cases in which the major
semiaxis is B -- 0.2a, 0e = 60 °, and the finite element mesh has ten nodes in the r direction and ten in the 0 direction. The matching radius was so chosen in Cases 3 and 4 that the derivatives of the solution are less than 0.1, in keeping with the presumption of low meniscus slope at r - - R .... The computed contact angles in Table I differ from 60 ° because the finite-element solution is approximate; i.e., it satisfies the differential equation and boundary condition only in a certain average sense. The computed contact angles give an indication of the error in the slope of the finite element solution. Case (2) demonstrates that without the matching procedure, reasonable results can be obtained, though at the expense of increased element size when Rn, is large. Cases (3) and (4) show the additional accuracy gained by using the matching technique, which allows the use of a smaller domain (smaller R,,,). Clearly, the finite element solutions obtained here for the limiting case of ~ = 1 agree well with the more exact results of H u h and Scriven. When the cylinder is elliptical the meniscus loses its axial symmetry, as shown in Fig. 2.
h.s.
"
,_-+-
141
:
I
'k
\ ": ~"
I\
\
\
FIG. 2. Oblique view of a meniscus in an infinite pool about a cylinder of elliptical cross section. Elevation contours are shown. ~ ~- 0.6, B = 0.2a, 0c = 20 °. Journal of Colloid and Interface Science, Vol. 60, No. 1, J u n e 1, 1977
142
ORR, BROWN AND SCRIVEN
FIG. 3. Oblique view of half the meniscus between and around two cylinders in an infinite pool (the cut is along the plane of mirror symmetry). Elevation contours are shown. R = 0.2a, L = 0.3a, Rm = 0.8a, 0 = 40°. I t attains its maximum height where the minor axis of the elliptical cross section intersects the solid surface; its minimum rise occurs at the major axis of the ellipse, which is where the solid has maximum curvature. The variation in height of rise around the solid surface is induced by the coupling of the shape of the liquid surface to that of the solid surface in the contact angle boundary condition. The variation in surface elevation diminishes smoothly with increasing distance from the cylinder, in keeping with the elliptic nature of the problem. The numerical solution technique was further validated by solving for the shape of a meniscus elevated between two vertical circular cylinders, a problem considered by Princen
(9). The surface shape depends on the spacing between cylinders, L, as well as the cylinder radius WR = B and the contact angle 6¢. When the cylinders are far apart the meniscus around each differs little from that around a lone cylinder, but when the cylinders are close the meniscus climbs sharply to a saddle shape in the region where their surfaces are nearest together, as can be seen in the specimen solutions plotted in Figs. 3 and 4. The existence of such a configuration between two cylinders is the basis of a force balance used by Princen (9) to predict the capillary rise between an isolated pair of cylinders. Princen considered the case when the capillary rise is many times the cylinder radius and the meniscus can be approximated by a surface which turns upward with ever increasing steepness as it approaches the saddle region. In this case the height at the lowest point in the saddle can be estimated by balancing the surface tension force across the meniscus at that elevation against the weight of the column of liquid below. Again, the agreement between the numerical and approximate solutions is good. When B = 0.4a, L = 0.48a, RM = 1.6a and/9¢ = 15 °, the finite-element program calculated a meniscus height at 0 = 0, r = 0 which differs by no more than 4~/o from Princen's approximation. That the discrepancy is not less is attributable to the meniscus slope just beneath the saddle region being less than infinite in reality, whereas Princen's approximation takes it to be infinite. Compare Fig. 5. LATERALLY FREE CONTACT LINES If the solid cylinders of the previous section are replaced by objects having cross sections that vary with elevation, the horizontal location of the contact line at which the meniscus intersects the solid objects is not known a priori. The contact line is in this sense laterally free ; in other words the LaplaceYoung equation has to be solved with a free boundary. Free boundary problems are frequently difficult and in fact no three-dimen-
Journal of Colloid and InterJace Science, Vol. 60, No. I, J u n e 1, 1977
143
THREE-DIMENSIONAL MENISCI
FIG. 4. Plan view of contours of surface elevation. Note the saddle and ramp regions between cylinders. The contours grow more nearly elliptical with distance from the cylinders and approach circularity at large distance. R = 0.2a, L = 0.3a, Rm = 0.8a, 0 = 40 °.
sional solution of the Laplace-Young equation with a free boundary has heretofore been published (the closest attempt appears to be that of Hinata et al. (10), who converted a minimal surface problem into a free surface problem which they solved by a finite-element method). For certain simple solid shapes, free boundary problems can be handled in a straightforward way. An iterative scheme which generates successive approximations to the free boundary is (a) Guess a solution for the elevation of the surface. (b) Compute a location of the solid-liquidgas contact line, thereby setting a domain for the current iteration. Compute mesh point locations for that domain. (c) Solve the meniscus confguration problem on the current domain using the method
The successive approximation scheme of (1) can be used to circumvent the nonlinearity of Eq. [21-] below. In fact, the computer program used for the solutions presented here is a modification of that discussed in (1) Esee (2) for additional details-]. In order to illustrate the solution of free boundary problems we return to the arrangement of solid objects in an infinite square array studied in (1), but with the circular cylinders replaced by either cones or spheres. The differential equation E1-] is unchanged but the local radius of the solid is now a function of elevation. The finite-element analog of the LaplaceYoung equation is as above with one exception. The normal to the solid boundary n is no longer simply the negative of rib, the normal to the domain boundary, because n is now inclined out of the horizontal. Thus the boundary integral
of (1). (d) Use new elevations in step (b) and repeat to convergence.
Bi =- rod ND ¢inds
[173
Journal of Colloidand InterfaceScience,Vol. 60, No. 1, June 1, 1977
ORR, BROWN AND SCRIVEN
144
integrals are f D (~pih -- ~ 7 i i ~ i . n ) d a + B i = O,
i=l,
2,...,Q,
[21]
where the Bi are given by Eq. [20]. S P E C I M E N M E N I S C I I N ARRAYS OF CONES A N D S P H E R E S
At elevation h the horizontal radius of a vertical cone is
r(h) = R + h tan/3,
[-223
where R is the cone radius at zero elevation and 5 is the angle between the generator and axis of the cone. At elevation h the horizontal radius of a sphere is
r(h) = [R ~ - (h --h0)2-] ;,
[23-]
where R is the sphere radius and h0 the elevation of the sphere center. The normal to the surface of either solid is
e,n
FIO. 5. Oblique view of half the meniscus between and around two cylinders in an infinite pool (the cut is along the plane of mirror symmetry). Elevation contours are shown. R = 0.4a, L = 0 . 4 8 a , 0 = 15 °, Rm = 1.6a.
changes. Let t~ be the tangent to the contact line. The normal to the domain boundary nb is the normal to the contact line which is also parallel to the base plane : nb -- (k >< t e ) / l k X to[.
[18-]
The tangent to the contact line is given by t~ = N X n.
[19"]
Hence Eq. [17] becomes Bi=
~ k x f. DlkXt.I *'
( n X N).nds.
[20]
Thus for solids of general shape arranged in an infinite square array, the weighted residual
(dr/dZ)k
=
.
[24-]
E1 + (dr/dZ)2] ~ The analysis of this section can be applied only if the surface representation remains single-valued in the coordinate system used for the computations. This condition is satisfied as long as 0 c - fl >__ 0 in the case of cones, and when 7r/2 -- arccos [-(h0 -- he)/R] < 0c for spheres (he is the elevation of the contact line). A solution for an array of cones is shown in Fig. 6. The semivertex angle of the cone is /3 = 20°; the half-spacing between centerlines, L = a; the cone radius at zero elevation, R -- 0.04a; and the contact angle, 0c = 20 °. The surface displays the same geometric features as the meniscus in an array of circular cylinders (1). The maximum elevation occurs where the line of closest approach intersects the cone, and the minimum is at the center of symmetry between cones. However, the contact line is a free boundary and in the plan view is no longer a circle.
Journal o.[ Colloid and D terJace Science, Vol. 60, No. 1, June 1, 1977
THREE-DIMENSIONAL MENISCI
145
Fro. 6. Oblique view of a meniscus in an array of cones on square centers. Three adjacent primitive cells around one of the cones are shown. R = 0.1a, L = a, 0 = 50°, fl = 20 °. The same features characterize the solution for an array of spheres. Fig. 7 is the plan view of the meniscus between spheres of radius R = 0.9a centered at h0 = 1.8a a n d spaced at L = a. The contact angle is 0c = 60 °. The upper part of Fig. 8 is an oblique view of the same meniscus. As the elevation h0 of the sphere centers decreases, the meniscus becomes less highly curved. If the solution continues to exist for h 0 ' = h o ~ n2R, n a n integer, then we can regard the infinite square array as a layer in a three-dimensional array of spheres and we see that more than one meniscus configuration can exist for a given pressure at the d a t u m elevation h -- 0. Consider, for instance, a layer of spheres of radius R and center elevation h0 = 0 with another layer of spheres of the same radius just above the first, with h0' = 2R. One equilibrium solution can be found for the meniscus in contact with only the first level of spheres, and a second can be found when the meniscus contacts only the second level of spheres. Which solution is obtained depends FIG. 7. Plan view of the meniscus between spheres in a square array. Only one primitive cell is shown. Interior curves are elevation contours. R = 0.9a, L = a, 0c = 60°,h0 = 1.8a.
on the history of the meniscus: if it is percolating down through the sphere pack the upper solution is reached first; if it is migrating upward the lower solution is attained first. I n
Journal of Colloid and Interface Science, Vol. 60, No. 1, J u n e 1, 1977
146
ORR,~BROWNZAND S C R I V E N
FIG. 8. Oblique view of two alternative meniscus configurations in a three-dimensional array of spheres. The horizontal layers are arrays on square centers and they are in vertical registry. Elevation contours are shown. R = 0.9a, L = a, 0c = 60°; ho = 0 for the lower meniscus and ho' = 1.8a for the upper one.
the computer simulation, upper solutions were obtained by finding the lower equilibrium configuration at h0 = 0 and then solving a sequence of cases in which the elevation of the sphere center was steadily increased until h0r = 2R was reached. In Fig. 8 the upper solution, hop = 1.8a, and the lower solution, h0 = 0, are shown for spheres of the same radius, R = 0.9a. Note that the upper meniscus is more highly curved, as it must be in order to support the greater weight of liquid above the datum elevation. The computer simulation has shown no additional solutions in this case, although they may exist in others, depending on the Bond number ( R / a ) ~ and the contact angle 0o. The existence of multiple solutions of the
type discussed here was noticed by Haines (11). This was discussed and demonstrated experimentally by Smith et al. (12) ; see also (13, 14). Numerical simulation by finite elements opens the way to detailed investigations of menisci in sphere packs and less regular porous media. However, the iterative scheme employed here proved to be much more sensitive in calculations for spheres than in calculations for cones. The reason is that the normal to a cone surface does not vary with elevation while that to a sphere surface does, and so in the latter case the form of the contact angle boundary condition changes with elevation as well as azimuthal angle. Underrelaxation was often required in Step (d) above, in order to limit the extent of domain alteration at each
Journal o] Colloid and Interface Science, Vol. 60, No. 1, June 1, 1977
147
THREE-DIMENSIONAL MENISCI
iteration. Otherwise there were severe oscillations in the sequence of elevations produced by the successive approximation process. CONCLUSIONS
The solutions to the Laplace-Young equation which are presented here demonstrate the power of the finite-element method for predicting and analyzing menisci which are neither translationally nor rotationally symmetric. Not only can fully three-dimensional surface shapes be computed, but also the further complications brought by menisci of unbounded extent and by free contact lines can be dealt with effectively. Indeed, solutions for three-dimensional menisci of infinite extent and for three-dimensional menisci with a free boundary are reported here for the first time. To obtain them we have relied heavily on the computational convenience of expanding the surface elevation in trial functions on subdomains, on the ease with which irregular domains can be approximated, and on the substantial simplification of natural boundary conditions such as the contact angle condition--all of which are particular advantages of the finite-element approach. It seems reasonable to expect that this approach will bring within reach the simulation of still more complex meniscus configurations that are important in science and technology.
ACKNOWLEDGMENTS This work is supported in part by the National Science Foundation and also in part by the University of Minnesota Computer Center. The second author is grateful for an Eastman Kodak Fellowship awarded through the Department of Chemical Engineering and Materials Science. REFERENCES 1. ORR,F. M., JR., SCRIVEN,L. E., AND RIVAS,A. P., J. Colloid Interface Sci. 52, 602 (1975). 2. ORR, F. M., JR., Ph.D. Thesis, University of Minnesota, 1976. 3. ZL~.a~AL,M., Numer. Math. 14, 394 (1970). 4. ZLAMAL, M., Int. J. Numer. Meth. Eng. 5, 367 (1973). 5. ZLAMAL, M., S I A M J. Numer. Anal. 10, 229 (1973). 6. GUSTAVSON,F. G., in "Sparse Matrices and Their Applications" (D. J. Rose and R. A. Willoughby, Eds.), pp. 41-52. Plenum, New York, 1972. 7. Hun, C., ANDSCRIVEN,L. E., J. Colloid Interface Sci. 30, 323 (1969). 8. STROUD, A. H., "Approximate Calculation of Multiple Integrals," Prentice-Hall, Englewood Cliffs, N. J., 1971. 9. PRINCEN, H. M., J. Colloid Interface Sei. 30, 69 (1969). 10. HINATA, M., SHIMASAKI,M., AND KIYONO, T., Math. Comp. 28, 45 (1974). 11. HAINES, W. B., J. Agric. Sci. 20, 97 (1930). 12. SMITrL W. O., FOOTE, P. D., AND BUSANG,P. F., J. Appl. Phys. 1, 18 (1931). 13. SMITH,W. 0., J. Appl. Phys. 4, 184 (1933). 14. MORROW,N. R., Ind. Eng. Chem. 62, 32 (1970).
Journal of Colloid and InterJace Science, Vol. 60, No. 1, June 1, 1977
Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, Minnesota 55455 Received May 22, 1976; accepted July 12, 1976 Menisci with complicated shapes occur in wetting and drying of powders and fabrics, coating of surfaces, fluid displacement in porous media, and other technologies. To predict meniscus shape and location requires the solution of the Laplace-Young equation, which is nonlinear and elliptic and depends on two space variables. The finite-element method has particular advantages that are further demonstrated here. Its recent application to menisci between cylinders is extended to three-dimensional menisci of unbounded extent or bounded by free contact lines (free boundary problems). Specimen solutions are displayed, including multiple solutions for the meniscus in a simple threedimensional array of spheres, an elementary type of porous medium. INTRODUCTION Menisci that possess neither translational nor rotational symmetry are three-dimensional. Although governed by the well-known LaplaceYoung equation, their full description requires the numerical solution of this nonlinear, elliptic equation. In the antecedent I (1), the finiteelement method was shown to be particularly effective because it can easily accommodate the irregular boundaries and nonlinear boundary conditions that generally accompany threedimensional menisci. Although this method, as developed in (1), can be used to solve the Laplace-Young equation for more complicated meniscus shapes than any other available method, there are important meniscus problems that cannot be solved by the method as there constituted. There are three sorts of difficulties that arise: (1) menisci of unbounded extent, i.e., domains of infinite expanse; (2) free contact lines, i.e., domain boundaries the location of which is part of the problem; and (3) menisci which do not have a single-valued representation in the customary Cartesian coordinate
system or an equivalent system. I n this paper we extend the finite-element method in order to deal with the first two difficulties and by way of illustration we solve several problems of scientific or practical interest. DIFFERENTIAL EQUATION AND BOUNDARY CONDITIONS In the notation used in I, the LaplaceYoung equation is 2(H -- t!o)g = g(Z -- Zo)Ap
[11
and the mean curvature of the surface is proportional to the surface divergence of the unit vector N, which is everywhere normal to the surface: 2H = - V s . N . When the surface can be represented by a single-valued elevation Z above a datum plane, it happens that V~-N = VII.N, where VII is the gradient operator in the datum plane. If that plane is parameterized with Cartesian coordinates (x, y), VII = i o / o X + jo/Oy; if polar coordinates (r, O) are used, VII = era~Or + eor-lO/O0. 137
C o p y r i g h t ~ 1977 by A c a d e m i c Press, Inc. All rights of reproduce.ion in a n y form reserved,
Journal oJ Colh)id and Interface Science, Vol. 60, No. l, J u n e 1, 1977 I S S N 0021-9797
138
ORR, BROWN AND SCRIVEN the expansion is [-cf. (1)]
In either case N = (VnZ)[-1 + (VliZ)2]-½.
[-2-]
3
h = E hi~, ee) (f, ,)
When the effect of gravity is significant, as in the cases considered here, it is convenient to scale all lengths with the capillary length a--(a/gAp)½. This is done in I and the Laplace-Young equation is shown to become VII" N + h = 0,
N = (k - ih~ - j h , ) / ( 1 + h~~ +
x ~-----,
h,~)½,
y Z -- Z0 •-------, h-~ +2H0a.[3]
Where the meniscus intersects a solid surface three types of boundary conditions are possible: a Dirichlet condition in which the location of the contact line is specified; a Neumann condition in which the contact angle between the fluid interface and the solid surface is specified; or a Robin condition in which location and angle are related. In the simplest instance the Neumann or contact angle boundary condition is N . n = cos 0o
[-4]
where n is the normal to the solid surface and 0¢ is a characteristic contact angle, sometimes held to be an equilibrium property, exhibited by the solid-liquid-fluid trio. Owing to the nonlinearity of N the contact angle condition is fairly difficult to handle. We adopt it for the sample problems solved here.
(~ and V refer to dimensionless Cartesian or cylindrical polar coordinates). We have experimented with linear trial functions (1) and with reduced Hermite cubic trial functions (3-5). We found that the increased accuracy of this cubic interpolation does not justify its added cost in computation time; that is, comparable accuracy can be obtained at lower cost with linear elements on a more refined mesh. Therefore, the results presented here were all obtained with linear trial functions. It is possible, however, that other higher order trial functions would have advantages. The Galerkin form of the finite-element method selects the approximate solution which makes the residual in the differential equation orthogonal to each of the trial functions
fa
4Ji(V.N+h)dA = 0 , i = 1 , 2 , . . . , Q , [6]
where A is the area of the domain and Q is the number of trial functions. With linear elements Q equals the number of nodes. An integration by parts followed by an application of the divergence theorem generates the nonlinear algebraic equations used to determine the expansion coefficients in Eq. [-5]: f
fo nb'dpiNds + I A
NUMERICAL METHOD The finite-element method used here follows closely that derived in I. Details are given in (2). The elevation of the surface is sought as an expansion in terms of locally defined polynomial trial functions. The problem domain is divided into triangular subdomains, or elements, and within each element a small number of nonzero trial functions is used to interpolate the elevations at the triangle vertices (nodes). Thus in a particular element
[-53
Jx
(4~ih - V$i'N)dA = 0, i = 1,2,...,Q,
[73
where nb is the normal to the domain boundary OA. The boundary conditions enter the approximate equations through the boundary integral. Equations [-7"]for the expansion coefficients hi are a nonlinear algebraic set which has the form
Journal of Colloid and Interface Science, Vol. 60, No. 1, June 1, 1977
F(h) = O.
E83
THREE-DIMENSIONAL MENISCI This set can be solved by successive approximations, which was done in I, or'by'means of a Newton-Raphson iteration which estimates the correction to the vector of expansion coefficients at the (k 4- 1)st iteration as
139
INTERIOR DOMAIN
DOMAIN y
h~k+l~ = h ( ~ -
F
OFi-1-t
L
OhsJ
l--/
F(h(k)) •
I-9-1
Whichever iterative scheme is used, the matrices which must be inverted are quite sparse ; they have a maximum of seven nonzero coefficients on any row. We found the sparse matrix techniques discussed by Gustavson (6) to facilitate greatly the required Gaussian eliminations. INTERFACES
OF UNBOUNDED
Fro. 1. Single cylinder of elliptical cross section in a pool of infinite extent. Sample curvilinear triangular
EXTENT
If a meniscus extends far enough outward from solid bodies it contacts, its equilibrium shape in a gravitational field becomes sensibly flat at some distance from its line of contact with the solid. From a mathematical point of view the interface may then be regarded as unbounded. But to impose a boundary condition at infinity is not possible in a finite set of finite elements. The simplest approximation is to impose the condition at some finite yet large distance, and then to increase the distance until it is sufficiently large that further increase leaves the problem solution unchanged at the desired level of accuracy. The disadvantage of this scheme is that a very large domain may be required, which in turn requires that a very large numerical problem be solved to the desired level of accuracy. The disadvantage can be circumvented by using knowledge of the asymptotic solution for the nearly flat meniscus far from solid bodies. Now lim Z (x, y) = Z0, OZ OZ lim-- = 0 = lim--, ~ Ox r ~ Oy
OA s
El0]
where r = (x 2 + y2)~ measures distance from solid bodies and Z0 is the elevation of a
finite e l e m e n t s d e r i v e d shown. , ~ A/B.
from
polar
coordinates are
perfectly flat meniscus. Thus at large distances r we have (OZ/Ox)2<< 1, (OZ/Oy)2<< 1, Ho = O, Z0 = O, and Eq. [1] becomes ~Z
=
gZAp.
Ell-]
This is Helrnholtz's equation, a well-known linear equation for the solution of which there are standard analytic methods. Its solution satisfying a boundary condition at infinity still contains unknown constants; these can be eliminated by matching the exterior asymptotic solution to an interior approximate solution of the full Eq. E1-]. To demonstrate this approach we use the meniscus around a cylindrical object standing vertically in a pool of unbounded extent. If the cylinder has a circular cross section, the meniscus is axisymmetric, and that problem was solved numerically by Huh and Scriven (7). Otherwise, the meniscus is three-dimensional, and no solution of such a problem on an infinite domain has appeared. A sample domain for a cylinder of elliptical cross section is shown in Fig. 1. Owing to its symmetry the problem can be solved in one quadrant. The domain is divided into an interior portion where the nonlinear, full problem is solved by the finiteelement method, and an exterior portion where
Journal of Colloid and Interface Science, Vol. 60, No. 1, June I, 1977
140
ORR, BROWN-AND-SCRIVEN TABLE I
Comparison of Results Computed for Axisymmetric Menisci around a Cylinder Case
R~
C o n d i t i o n a t r = R~
0 I 2 3 4
-a 3a 1.4a a
H u h and S c r ive n (7) 0 = 7r/2 0 = 7r/2 M a t c h i n g condition M a t c h i n g condition
Height at cylinder
Error (%) (vs Case 0)
Computed contact a ngle (degrees)
0.189 0.306 0.187 0.188 0.188
-100 1.09 0.68 0.53
60 60.8 62.6 61.1 60.7
the linear , asymptotic problem is solved. The two problems are coupled at the matching boundary, a circle of radius r = Rm where the matching condition is that the solution and its first derivatives are continuous. The exterior problem is V2h -- h = 0,
r >
Rln,
Oh
Oh O0
oa rhinb"Nds
lira h = l i r a - - = l i m - - = 0, r~Qo
r~:
Or
formly for r > Rm, yet its r derivative does not converge at r = Rm. The values of Oh/Or at r = Rm, which are needed for the matching condition, are therefore approximated with a forward finite difference just outside the matching boundary. When (OhiO0)2<<1, this series converges rapidly. The data from the asymptotic solution enter the interior problem through the boundary integral of [7-]. Integration is over the solid boundary for which the contact angle condition V4-] is given, over the matching boundary r = Rm on which the normal to the surface is given by the exterior solution, and over the two symmetry boundaries 0 = 0 and 0 = 7r/2. The contributions along the symmetry boundaries vanish because nb.N is always zero there. Thus the boundary contribution to the weighted residual integrals is
,~o
=
fOAs rhinb"Nds
[12]
Oh O0
-0
h = P(O)
at
0 = 0
and
at
r = Rm,
where P(O) is the finite-element expansion for the interior solution evaluated at the matching boundary. Solving the exterior problem with a finite Fourier transform yields a series solution in terms of modified Bessel functions of the second kind:
K2.(r) h(r, O)= ~, a n p , , n=o Kz,(Ro)
+ ] 4~inm"Nds. [14] .IoAm
0 = 7r/2,
On the solid boundary OAs, nb = --ns and on the matching boundary OAm,nm = er. Inserting the contact angle condition and evaluating nn~"N yields
Y fo ~inb"Nds = -- I A
+ fo
~o
cos (2n0)
[-13]
t2/Trl,
n = 0 n =
COS
Oo4s
4)ihrds [15] A,,~[1 + h~ + (ho/r)q~"
Because the square of the slope is supposed to be small at r = Rm, this can be simplified to
where 2/Tr)b
,,tOA s
1, 2, . . . ,
foA 4)inb"Nds -
--
fOAs dpiCOSOcds
and -]- f O A m ~ i h , d s .
Tr[2
p,, =--a.
f
P(O) cos (2nO)dO.
dO
The series and its derivatives converge uniJournal of Colloid and Interface Science,
[-16]
Equations [-7] and [-16] make up the nonlinear algebraic set to be solved for the expan-
Vol. 60, No. 1, J u n e 1, 1977
THREE-DIMENSIONAL MENISCI sion coefficients hi. A three-point Gaussian quadrature formula (8) was used to evaluate area integrals in each element and boundary integrals were evaluated analytically. The Newton algorithm [14"] converged quickly for all cases tested. Typically, three to five iterations were required to reduce the maxim u m adjustment in the h~ to less than 10-4. SPECIMEN MENISCI AROUND CYLINDERS When the ellipse indicated in Fig. 1 has ellipticity ~ = 1 the solid is a circular cylinder of radius B. If the contact angle e~ is uniform all around the cylinder the meniscus extending outward is axisymmetric and the finite-element solution should agree with the results of H u h and Scriven (7); this provides one test of the method adopted here. In addition, the matching strategy used here can be tested by setting the slope equal to zero at the matching radius r = R .... This corresponds to solving the problem of a meniscus in the annulus between r = B and R = R .... with contact angles of 0~ on the inner cylinder and 7r/2 on the outer wall. Table I compares with the more exact solution a number of cases in which the major
semiaxis is B -- 0.2a, 0e = 60 °, and the finite element mesh has ten nodes in the r direction and ten in the 0 direction. The matching radius was so chosen in Cases 3 and 4 that the derivatives of the solution are less than 0.1, in keeping with the presumption of low meniscus slope at r - - R .... The computed contact angles in Table I differ from 60 ° because the finite-element solution is approximate; i.e., it satisfies the differential equation and boundary condition only in a certain average sense. The computed contact angles give an indication of the error in the slope of the finite element solution. Case (2) demonstrates that without the matching procedure, reasonable results can be obtained, though at the expense of increased element size when Rn, is large. Cases (3) and (4) show the additional accuracy gained by using the matching technique, which allows the use of a smaller domain (smaller R,,,). Clearly, the finite element solutions obtained here for the limiting case of ~ = 1 agree well with the more exact results of H u h and Scriven. When the cylinder is elliptical the meniscus loses its axial symmetry, as shown in Fig. 2.
h.s.
"
,_-+-
141
:
I
'k
\ ": ~"
I\
\
\
FIG. 2. Oblique view of a meniscus in an infinite pool about a cylinder of elliptical cross section. Elevation contours are shown. ~ ~- 0.6, B = 0.2a, 0c = 20 °. Journal of Colloid and Interface Science, Vol. 60, No. 1, J u n e 1, 1977
142
ORR, BROWN AND SCRIVEN
FIG. 3. Oblique view of half the meniscus between and around two cylinders in an infinite pool (the cut is along the plane of mirror symmetry). Elevation contours are shown. R = 0.2a, L = 0.3a, Rm = 0.8a, 0 = 40°. I t attains its maximum height where the minor axis of the elliptical cross section intersects the solid surface; its minimum rise occurs at the major axis of the ellipse, which is where the solid has maximum curvature. The variation in height of rise around the solid surface is induced by the coupling of the shape of the liquid surface to that of the solid surface in the contact angle boundary condition. The variation in surface elevation diminishes smoothly with increasing distance from the cylinder, in keeping with the elliptic nature of the problem. The numerical solution technique was further validated by solving for the shape of a meniscus elevated between two vertical circular cylinders, a problem considered by Princen
(9). The surface shape depends on the spacing between cylinders, L, as well as the cylinder radius WR = B and the contact angle 6¢. When the cylinders are far apart the meniscus around each differs little from that around a lone cylinder, but when the cylinders are close the meniscus climbs sharply to a saddle shape in the region where their surfaces are nearest together, as can be seen in the specimen solutions plotted in Figs. 3 and 4. The existence of such a configuration between two cylinders is the basis of a force balance used by Princen (9) to predict the capillary rise between an isolated pair of cylinders. Princen considered the case when the capillary rise is many times the cylinder radius and the meniscus can be approximated by a surface which turns upward with ever increasing steepness as it approaches the saddle region. In this case the height at the lowest point in the saddle can be estimated by balancing the surface tension force across the meniscus at that elevation against the weight of the column of liquid below. Again, the agreement between the numerical and approximate solutions is good. When B = 0.4a, L = 0.48a, RM = 1.6a and/9¢ = 15 °, the finite-element program calculated a meniscus height at 0 = 0, r = 0 which differs by no more than 4~/o from Princen's approximation. That the discrepancy is not less is attributable to the meniscus slope just beneath the saddle region being less than infinite in reality, whereas Princen's approximation takes it to be infinite. Compare Fig. 5. LATERALLY FREE CONTACT LINES If the solid cylinders of the previous section are replaced by objects having cross sections that vary with elevation, the horizontal location of the contact line at which the meniscus intersects the solid objects is not known a priori. The contact line is in this sense laterally free ; in other words the LaplaceYoung equation has to be solved with a free boundary. Free boundary problems are frequently difficult and in fact no three-dimen-
Journal of Colloid and InterJace Science, Vol. 60, No. I, J u n e 1, 1977
143
THREE-DIMENSIONAL MENISCI
FIG. 4. Plan view of contours of surface elevation. Note the saddle and ramp regions between cylinders. The contours grow more nearly elliptical with distance from the cylinders and approach circularity at large distance. R = 0.2a, L = 0.3a, Rm = 0.8a, 0 = 40 °.
sional solution of the Laplace-Young equation with a free boundary has heretofore been published (the closest attempt appears to be that of Hinata et al. (10), who converted a minimal surface problem into a free surface problem which they solved by a finite-element method). For certain simple solid shapes, free boundary problems can be handled in a straightforward way. An iterative scheme which generates successive approximations to the free boundary is (a) Guess a solution for the elevation of the surface. (b) Compute a location of the solid-liquidgas contact line, thereby setting a domain for the current iteration. Compute mesh point locations for that domain. (c) Solve the meniscus confguration problem on the current domain using the method
The successive approximation scheme of (1) can be used to circumvent the nonlinearity of Eq. [21-] below. In fact, the computer program used for the solutions presented here is a modification of that discussed in (1) Esee (2) for additional details-]. In order to illustrate the solution of free boundary problems we return to the arrangement of solid objects in an infinite square array studied in (1), but with the circular cylinders replaced by either cones or spheres. The differential equation E1-] is unchanged but the local radius of the solid is now a function of elevation. The finite-element analog of the LaplaceYoung equation is as above with one exception. The normal to the solid boundary n is no longer simply the negative of rib, the normal to the domain boundary, because n is now inclined out of the horizontal. Thus the boundary integral
of (1). (d) Use new elevations in step (b) and repeat to convergence.
Bi =- rod ND ¢inds
[173
Journal of Colloidand InterfaceScience,Vol. 60, No. 1, June 1, 1977
ORR, BROWN AND SCRIVEN
144
integrals are f D (~pih -- ~ 7 i i ~ i . n ) d a + B i = O,
i=l,
2,...,Q,
[21]
where the Bi are given by Eq. [20]. S P E C I M E N M E N I S C I I N ARRAYS OF CONES A N D S P H E R E S
At elevation h the horizontal radius of a vertical cone is
r(h) = R + h tan/3,
[-223
where R is the cone radius at zero elevation and 5 is the angle between the generator and axis of the cone. At elevation h the horizontal radius of a sphere is
r(h) = [R ~ - (h --h0)2-] ;,
[23-]
where R is the sphere radius and h0 the elevation of the sphere center. The normal to the surface of either solid is
e,n
FIO. 5. Oblique view of half the meniscus between and around two cylinders in an infinite pool (the cut is along the plane of mirror symmetry). Elevation contours are shown. R = 0.4a, L = 0 . 4 8 a , 0 = 15 °, Rm = 1.6a.
changes. Let t~ be the tangent to the contact line. The normal to the domain boundary nb is the normal to the contact line which is also parallel to the base plane : nb -- (k >< t e ) / l k X to[.
[18-]
The tangent to the contact line is given by t~ = N X n.
[19"]
Hence Eq. [17] becomes Bi=
~ k x f. DlkXt.I *'
( n X N).nds.
[20]
Thus for solids of general shape arranged in an infinite square array, the weighted residual
(dr/dZ)k
=
.
[24-]
E1 + (dr/dZ)2] ~ The analysis of this section can be applied only if the surface representation remains single-valued in the coordinate system used for the computations. This condition is satisfied as long as 0 c - fl >__ 0 in the case of cones, and when 7r/2 -- arccos [-(h0 -- he)/R] < 0c for spheres (he is the elevation of the contact line). A solution for an array of cones is shown in Fig. 6. The semivertex angle of the cone is /3 = 20°; the half-spacing between centerlines, L = a; the cone radius at zero elevation, R -- 0.04a; and the contact angle, 0c = 20 °. The surface displays the same geometric features as the meniscus in an array of circular cylinders (1). The maximum elevation occurs where the line of closest approach intersects the cone, and the minimum is at the center of symmetry between cones. However, the contact line is a free boundary and in the plan view is no longer a circle.
Journal o.[ Colloid and D terJace Science, Vol. 60, No. 1, June 1, 1977
THREE-DIMENSIONAL MENISCI
145
Fro. 6. Oblique view of a meniscus in an array of cones on square centers. Three adjacent primitive cells around one of the cones are shown. R = 0.1a, L = a, 0 = 50°, fl = 20 °. The same features characterize the solution for an array of spheres. Fig. 7 is the plan view of the meniscus between spheres of radius R = 0.9a centered at h0 = 1.8a a n d spaced at L = a. The contact angle is 0c = 60 °. The upper part of Fig. 8 is an oblique view of the same meniscus. As the elevation h0 of the sphere centers decreases, the meniscus becomes less highly curved. If the solution continues to exist for h 0 ' = h o ~ n2R, n a n integer, then we can regard the infinite square array as a layer in a three-dimensional array of spheres and we see that more than one meniscus configuration can exist for a given pressure at the d a t u m elevation h -- 0. Consider, for instance, a layer of spheres of radius R and center elevation h0 = 0 with another layer of spheres of the same radius just above the first, with h0' = 2R. One equilibrium solution can be found for the meniscus in contact with only the first level of spheres, and a second can be found when the meniscus contacts only the second level of spheres. Which solution is obtained depends FIG. 7. Plan view of the meniscus between spheres in a square array. Only one primitive cell is shown. Interior curves are elevation contours. R = 0.9a, L = a, 0c = 60°,h0 = 1.8a.
on the history of the meniscus: if it is percolating down through the sphere pack the upper solution is reached first; if it is migrating upward the lower solution is attained first. I n
Journal of Colloid and Interface Science, Vol. 60, No. 1, J u n e 1, 1977
146
ORR,~BROWNZAND S C R I V E N
FIG. 8. Oblique view of two alternative meniscus configurations in a three-dimensional array of spheres. The horizontal layers are arrays on square centers and they are in vertical registry. Elevation contours are shown. R = 0.9a, L = a, 0c = 60°; ho = 0 for the lower meniscus and ho' = 1.8a for the upper one.
the computer simulation, upper solutions were obtained by finding the lower equilibrium configuration at h0 = 0 and then solving a sequence of cases in which the elevation of the sphere center was steadily increased until h0r = 2R was reached. In Fig. 8 the upper solution, hop = 1.8a, and the lower solution, h0 = 0, are shown for spheres of the same radius, R = 0.9a. Note that the upper meniscus is more highly curved, as it must be in order to support the greater weight of liquid above the datum elevation. The computer simulation has shown no additional solutions in this case, although they may exist in others, depending on the Bond number ( R / a ) ~ and the contact angle 0o. The existence of multiple solutions of the
type discussed here was noticed by Haines (11). This was discussed and demonstrated experimentally by Smith et al. (12) ; see also (13, 14). Numerical simulation by finite elements opens the way to detailed investigations of menisci in sphere packs and less regular porous media. However, the iterative scheme employed here proved to be much more sensitive in calculations for spheres than in calculations for cones. The reason is that the normal to a cone surface does not vary with elevation while that to a sphere surface does, and so in the latter case the form of the contact angle boundary condition changes with elevation as well as azimuthal angle. Underrelaxation was often required in Step (d) above, in order to limit the extent of domain alteration at each
Journal o] Colloid and Interface Science, Vol. 60, No. 1, June 1, 1977
147
THREE-DIMENSIONAL MENISCI
iteration. Otherwise there were severe oscillations in the sequence of elevations produced by the successive approximation process. CONCLUSIONS
The solutions to the Laplace-Young equation which are presented here demonstrate the power of the finite-element method for predicting and analyzing menisci which are neither translationally nor rotationally symmetric. Not only can fully three-dimensional surface shapes be computed, but also the further complications brought by menisci of unbounded extent and by free contact lines can be dealt with effectively. Indeed, solutions for three-dimensional menisci of infinite extent and for three-dimensional menisci with a free boundary are reported here for the first time. To obtain them we have relied heavily on the computational convenience of expanding the surface elevation in trial functions on subdomains, on the ease with which irregular domains can be approximated, and on the substantial simplification of natural boundary conditions such as the contact angle condition--all of which are particular advantages of the finite-element approach. It seems reasonable to expect that this approach will bring within reach the simulation of still more complex meniscus configurations that are important in science and technology.
ACKNOWLEDGMENTS This work is supported in part by the National Science Foundation and also in part by the University of Minnesota Computer Center. The second author is grateful for an Eastman Kodak Fellowship awarded through the Department of Chemical Engineering and Materials Science. REFERENCES 1. ORR,F. M., JR., SCRIVEN,L. E., AND RIVAS,A. P., J. Colloid Interface Sci. 52, 602 (1975). 2. ORR, F. M., JR., Ph.D. Thesis, University of Minnesota, 1976. 3. ZL~.a~AL,M., Numer. Math. 14, 394 (1970). 4. ZLAMAL, M., Int. J. Numer. Meth. Eng. 5, 367 (1973). 5. ZLAMAL, M., S I A M J. Numer. Anal. 10, 229 (1973). 6. GUSTAVSON,F. G., in "Sparse Matrices and Their Applications" (D. J. Rose and R. A. Willoughby, Eds.), pp. 41-52. Plenum, New York, 1972. 7. Hun, C., ANDSCRIVEN,L. E., J. Colloid Interface Sci. 30, 323 (1969). 8. STROUD, A. H., "Approximate Calculation of Multiple Integrals," Prentice-Hall, Englewood Cliffs, N. J., 1971. 9. PRINCEN, H. M., J. Colloid Interface Sei. 30, 69 (1969). 10. HINATA, M., SHIMASAKI,M., AND KIYONO, T., Math. Comp. 28, 45 (1974). 11. HAINES, W. B., J. Agric. Sci. 20, 97 (1930). 12. SMITrL W. O., FOOTE, P. D., AND BUSANG,P. F., J. Appl. Phys. 1, 18 (1931). 13. SMITH,W. 0., J. Appl. Phys. 4, 184 (1933). 14. MORROW,N. R., Ind. Eng. Chem. 62, 32 (1970).
Journal of Colloid and InterJace Science, Vol. 60, No. 1, June 1, 1977