Bifurcation of axially symmetric capillary surfaces

Bifurcation of Axially Symmetric Capillary Surfaces ULRICH H O R N U N G AND HANS DETLEF MITTELMANN 1 SCHI, P.O. Box 1222, W-8014 Neubiberg, Germany Received October 1, 1990; accepted February 26, 1991 Axially symmetric drops that are entrapped between two parallel planes and influenced by gravity and rotation are studied. Their shapes and the resulting bifurcation and stability diagrams are calculated numerically. The parameters involved are the volume and the gravitational and centrifugal forces. © 1991 Academic Press, Inc.

The paper is organized as follows: In Section 2 we state the variational problem and its Euler-Lagrange equation. The boundary value problem to be solved is derived and written in terms of an ODE. In Section 3 a method that uses a shooting technique for the solution of this problem is explained in detail. In Section 4 the numerical method used in this paper is described and the results are presented. This method is based on a finite difference approximation of the variational problem; it uses path-following techniques developed recently for bifurcation problems; see (9, 21 ). For an analysis of the bifurcation behavior of planar drops using singularity theory; see (19).

1. INTRODUCTION

It is known that certain capillary surfaces must be axially symmetric. For sessile, pendent, and rotating drops this was proved 10 years ago by Wente; see (29, 31 ). Arguments similar to those in these papers were used by Vogel for the case of drops entrapped between two parallel planes; see (26, 27). In these papers Vogel studied the stability under the assumption of zero gravity and without rotation. The present paper is devoted to the numerical calculation of bifurcation and stability diagrams for drops that are entrapped between two parallel planes, including the influence of gravity and rotation. It turns out that the diagrams change their type as soon as gravity is taken into consideration whereas rotation modifies the figures only quantitatively. For earlier results using perturbation methods and some computer simulation, see (24). Similar studies were carried out several years ago for drops rotating in free space; see (4, 6, 22, 28). Liquid bridges were investigated in (3, 7, 20, 23, 25). Also sessile and pendent drops have been analyzed with respect to their stability; see (30, 1, 2, 18). The behavior of liquids in axially symmetric containers was considered in ( 12, 14, 10).

2. THE BOUNDARY VALUE PROBLEM

We assume that the domain in question can be described as =

{ (X, y, Z)"

/9 =

VX 2 +

y2

< r ( z ) , 0 < z < 1}

with a function r: [0, 1] -+ R (see Fig. 1). The boundary of this domain is 0~2 = I? tO Z with r = {(x,y,z):p

= r ( z ) , O < ~ z < ~ 1}

and ~To whom correspondence should be addressed at Department of Mathematics, Arizona State University, Tempe, Arizona 85287-1804.

Z = {(x,y,z):p
= 1}.

219 0021-9797/91 $3.00 Journal of Colloid and Interface Science, Vol. 146, No. 1, October 1, 1991

Copyright © 1991 by Academic Press, Inc, All fights of reproduction in any form reserved.

220

HORNUNG AND MITTELMANN Z

E

a' = 1 R r 2 _ B z + p

sin a

1

r

The boundary conditions are cos a = - cos 0 for z = 0 and /

/ / / / /

/

i / / / / / /

Z

P

The dimensionless energy can be written as (Bz - ~Rp2)d~

This energy, expressed by integrals with respect to the variable z and divided by 2a-, is 2--g =

8zr

- g R r 4 + !/1 +

1 -~cos0(r2(0)

1.

Altogether, we arrive at the following boundary value problem:

FIG. 1. Notations.

E =

cosa=cos0forz=

I

r'2rldz

+ r2(1)).

a'= ½Rr:-Bz+p

sin a

1 r'

0
r' = cot a,

0
a=Tr-0,

z=0,

a=O,

z=

1, [3]

1.

In the special case B = 0 and R = 0 the problem reduces to finding rotationally symmetric surfaces of constant m e a n curvature, a problem which was studied by Delaunay in (8); its solutions can be given in closed form (see, e.g., (12)).

[ll 3. THE SHOOTING METHOD

The volume is V 21r

fo 1 ~1 r 2 d z .

[2]

Using standard arguments from the calculus of variations, one obtains for the variational problem with a volume constraint the differential equation d - - r' +

dz l 77ar,

1 rV1 + / 2 = 1Rr2 - Bz +p. 2

The left-hand side of the last equation m a y be written as d dzC°Sa+

sina r

(a, sina

We assume that the numbers 0, B, and R are given constants. We consider the initial value problem that is derived from problem [ 3 ] by dropping the condition at z = 1 and instead imposing another initial condition r = p for z = 0 with some initial value p. In general, the solution of this problem will not satisfy the boundary condition a = 0 for z = 1. In order to accomplish this, we modify the parameter p appropriately; i.e., we determine the n u m b e r p such that the condition at z = 1 is fulfilled. This can easily be done using a Newton iteration method. For this purpose we calculate the derivatives a p and rp of the functions a and r with respect to the parameter p. We get

1) +

,

which is identical to 2H, where H is the mean curvature of the surface P in 3D-space. Therefore, the resulting differential equation is Journal o f Colloid a n d Interface Sciencg, Vol. 146, No. 1, October 1, 1991

, 1 ap - sin a

(

Rrrp+ I-

(12 Rr

-Bz

o
BIFURCATION OF CAPILLARY SURFACES a, sin2a '

rp

a,=0,

O
4. NUMERICAL SOLUTION A N D STABILITY

1,

z=O,

rp=0,

[41

z=O.

The boundary value problem [31 is now equivalent to the equation G(p) = a( 1 ) - 0 = 0. The derivative of the function G is (d/ dp)G = ap(1). The next step is to modify the solution strategy such that both the numbers p and p are considered as parameters. This suggests the use of a path-following procedure. We assume that a pair (o, P) = (p0, Po) has been found for which the boundary value problem [ 3 ] has a solution. We want to find another pair (o, p) which approximately has the distance h from (p0, P0) in the O-p-plane. This is done by solving the pair of equations

G(p, P) = 0 c(p-po)+s(p-Po)-h

=0,

[5]

where C ~

OpG(po, Po) g

and

s

OpG(po, Po) g

with g = ]/(OeG(po, po)) z + (OpG(po, po)) z. The Jacobian of the system [ 5 ] is

There are at least two methods to numerically determine the shape of the capillary surfaces defined above. Both were utilized in the computations reported below, but the second approach outlined below was used for the determination of stability of solutions. One method is to realize the shooting method of Section 3 for which all the details including derivatives are given and for which standard ODE solvers can be called. The other way is to discretize the energy functional E in [ 1] and to minimize it subject to the constraint H = V - V0 = 0, with V as in [ 2 ]. If the computed stationary points are constrained minimizers of the energy, their stability can be determined. Assuming that asymmetric perturbations may be eliminated from consideration, a solution is called stable if it represents a local m i n i m u m of the functional [1 ] subject to the constraint [ 2 ]. This is analogous to the analysis given in (25) with the straightforward generalization of L e m m a 2.1 in (25) to the energy functional [ 1] containing the additional terms due to gravity and rotation. On a uniform grid zi = ih for i = 0 , . . . , n with h = 1 / n we introduce the grid function rh = (ro . . . . . rn), ri = rh( zi ), suppressing the h in the following. A straightforward difference discretization yields

E(r)

, =

1

sin oz

=

h

-

~Rr2 - B z

Rrrp-

-2 (zir i +

Z i _ 1 r i2- 1 )

i=l

The derivative of G with respect to p is found by solving the initial value problem ae

221

X

-

R (r 4 + 8

rL,)

+

ri

+ (ri + ri-1) --

:-,),

]2

-- h cos 0(r 2 + r z) and ,

re =

_ __ap

sin2a,

a e = O,

z = O,

re=l,

z=0.

O
1,

H ( r ) = h ~ r2i + r2i-1 2 i=1 [61

Now we obtain OeG = ae( 1 ) from system [ 6 ] and OpG = ap( 1 ) from system [ 4].

Vo.

The K u h n - T u c k e r conditions

VE +pVH= H=0

O,

[7]

Journal of Colloid and Interface Science, Vol, 146, No. 1, October 1, 1991

222

HORNUNG

AND

MITTELMANN

are used to find stationary points. This system may be solved by Newton's method. For this the Jacobian o

VH T

L = V2E + pV2H,

o f the system is needed. A solution r*, p * of [ 7 ] is a minimizer if S r L * S is positive definite, where L * = L(r*, p * ) and S i s an (n + 1 ) X n matrix whose columns form an orthonormal basis of the subspace M = { y E R n+l : V H ( r * ) T y = 0 }. Well-known methods from numerical linear algebra may be used to solve these problems. Here, the modified Gram-Schmidt process was used to compute S, and the Cholesky decomposition yielded the positive definiteness of S T L * S . For this, first an (n + 1 ) X (n + 1 ) matrix S'is formed by augmenting VII(r* ) by an (n + 1 ) × (n + I ) identity matrix o f which the, say, j t h column is omitted, wherej is such that t h e j th component of~TH(r* ) is not zero. Then the columns o f S are made mutually orthogonal. The well-known Gram-Schmidt process produces a factorization S = QR, where Q is an orthogonal and R an upper tri-

o

eq

i e5 0.03 0.13

0.32

0.51

0.71

0.90

V(B = 0., R = 0., A N G L E = 90)

FIG. 2. r ( 0 ) versus V, B = 0, R = 0, 0 = 9 0 °.

Journal o f Colloid and Interface Science. Vol. 146, No. 1, October 1, 1991

!3 0113 013:

0151

0:71

0:90

V(B = O., R = O., A N G L E = 9 0 )

FIG. 3. p versus V, B = 0, R = 0, 0 = 9 0 °.

angular matrix. The columns of Q, from the second to the last, may then be taken to define the desired matrix S. During this process, however, a severe loss of orthogonality may occur and it is thus advantageous to rearrange the order of the calculation yielding the modified Gram-Schmidt process; see Section 5.2.8 in(ll). The Cholesky decomposition (11) computes the decomposition A = GGT, where here A = S T L * S. It is clear that such a computation can succeed, producing a regular G, only ifA is positive definite, since xT A x = IIGx II2 > 0 for x :b 0. Thus, a premature termination, due to the need to compute the square root of a negative real number, indicates a lack of positive definiteness. In order to compute bifurcation diagrams the continuation routine ALCON2 (9) was called and the stability computation was integrated into this routine. This computational setup would allow one to solve more general problems; e.g., the contact angles with the two planes may be different as in (26), or nonaxisymmetric surfaces could be considered as in (13, 15-17). We preferred to study the case of nonzero Bond and rotational Bond numbers and to consider various quantities as continuation parameters.

BIFURCATION

OF CAPILLARY

SURFACES

223

c5

d ¸

e5 ¸

0.03 0.13

0.32

0~52

0171

0190

0.03 0.14

V(B = 0., R = 0., ANGLE= 80)

0.33

0.52

0.71

0.90

V(B = 0., R = 0., A N G L E = 80)

FIG. 4. r ( 0 ) versus V, B = 0, R = 0, 0 = 80 °.

FIG. 6. r ( 0 ) versus V, B = 1, R = 0, 0 = 80 °.

First we consider the case B = 0, R = 0. The diagrams in the (r, V)- and (p, V)-planes for 0 = 90 ° are given in Figs. 2 and 3. The following figures are given only for the (r, V)plane. If0 is varied from 90 °, then the second bifurcation where solutions with two inflection points bifurcate from the cylindrical solutions is perturbed; cf. Fig. 4. On the other hand, if B > 0 is chosen, then all bifurcations are per-

turbed; cf. Figs. 5,6. The choice of a nonzero R only quantitatively modifies the bifurcation diagrams. The stability computations yielded that only the solutions with no inflection points are stable down to a certain volume. For B = 0, R = 0, and 0 = 90 ° this volume is

f

?~.

¢5

m

0.03 0.14

0.'33

0.52

0171

0.00

0190

V(B = 0., R = 0 , A N G L E = 90)

FIG. 5. r ( 0 ) versus V, B = 1, R = 0, 0 = 90 °.

0.60

* 1.20 1.80 B (R = 0 , N = 30)

, 2.40

3.00

FIG. 7. A n g l e for smallest stable bridge as a function o f =0.

B,R

Journal of Colloid and Interface Science, V o l . 146, No. 1, October 1, 1991

224

HORNUNG AND MITTELMANN TABLEI

5. CONCLUSION

Smallest Stable Volume and Type of Point at Loss of Stability

Numerical methods that were designed for bifurcation problems with one or more parameters are used to study axially symmetric capillary surfaces. Diagrams are obtained for various sets of parameters, such as volume, Bond number, and rotational Bond number.

B

0

Vn,i.

Type

0 0 1

90 ° 80 ° 90 °

0.3184 0.3213 0.5608

Bifurcation point Bifurcation point Turning point

ACKNOWLEDGMENTS

known to be Vmin = 1/~r; see (25). In order to get an overall picture of the stability behavior, the stability limit Vmin w a s computed as a function of B, R, and 0. For these computations, the newest version of the program PITCON (20) was used since it allows better control of the continuation step length. A onesided bisection algorithm was used to find Vmin for decreasing values of the volume. The loss of stability occurred at the first bifurcation point (Figs. 2-4) and at the first turning point (Figs. 5,6) on the branch with solutions exhibiting at most one inflection point. On the bifurcating branches the number of inflection points was higher by one. As was found in (5) for the case of B = 0, R = 0, the m i n i m u m value o f V0 with respect to 0 is attained at 0 = 90 °. This is not true for B > 0. In fact, the minimizing 0 moves to considerably larger values while the effect of a nonzero rotational Bond number is a further shift, here also to larger angles. In order to see the change in the angle for the minimal stable volume with B, finally Fig. 7 was plotted showing the nonlinear growth o f 0 in the range 0 ~ B ~< 3. In order to compare our results with those known from the literature, we computed the volume at the point of loss of stability using a finer discretization. As quoted above from (25, Corollary 4.4), the m i n i m u m possible volume is 1/~r = 0.3183 (to four digits). Table I lists the stability limits obtained for R = 0 and a discretization with n = 50 and the type of singular point.

JournalofColloidandInterfaceScience,Vol. 146,No. 1, October1, 1991

The work of the second author was supported by the U.S. Air Force Office of Scientific Research under Grant AFOSR 90-0080. The authors thank Hua Chen for her help with the computations. REFERENCES 1. Babu, S. R., J. Colloid Interface Sci. 116, 350 (1987). 2. Bashforth, F., and Adams, J. G., "An Attempt to Test the Theories of Capillary Action by Comparing the Theoretical and Measured Forms of Drops of Fluid." Cambridge University Press, London/New York, 1883. 3. Bayramli, E., Abou-Obeid, A., and Van De Ven, T. G. M., J. Colloid lnterface Sci. 116, 490 (1987). 4. Brown, R. A., and Scriven, L. E., Proc. R. Soc. London, Set. A 371, 331 (1980). 5. Carter, W. C., Acta Metall. 36, 2283 (1988). 6. Chandrasekhar, S., Proc. R. Soc. London Ser. A 286, I (1965). 7. Chen, L.-H., and Chang, H.-C., J. Colloid Interface Sci. 120, 377 (1987). 8. Delaunay, C., J. Math. Pures Appl. 6, 309 (1841). 9. Deuflhard, P., Fiedler, B., and Kunkel, P., S l A M J. Nurner. Anal 24, 912 (1980). I0. Finn, R., Manuscr. Math. 61, 347 (1988). 11. Gohib, G. H., and Van Loan, Ch. F., "Matrix Computations," 2nd ed. Johns Hopkins Univ. Press, Baltimore/London, 1989. 12. Hornung, U., Technical Report 105, Department of Mathematics, Arizona State University, Tempe, 1988. 13. Hornung, U., GAMM-Mitteilungen 2, 33 (1989). 14. Homung, U. in "Free Boundary Problems. Theory and Applications" (K.-H. Hoffmann, Ed.), V o l . 5, p. 367. Pitman, Boston, 1990. 15. Hornung, U., and Langbein, D., in "Geometric Analysis and Computer Graphics" (P. Concus and R. Finn, Eds.), p. 1090. Mathematical Sciences Research Institute Series. Springer-Verlag, New York/ Berlin, 1991.

BIFURCATION OF CAPILLARY SURFACES 16. Hornung, U., and Mittelmann, H. D., J. Colloidlnterface Sci. 133, 409 (1989). 17. Hornung, U., and Mittelmann, H. D., J. Comput. Phys. 87, 126 (1990). 18. Kumar, A., and Hartland, S., J. ColloM Interface Sci. 124, 67 (1988). 19. Lewis, D., and Marsden, J., J. Math. Phys. 28, 2508 (1987). 20. Mazzone, D. N., Tardos, G. 1., and Pfeffer, R., J. Colloid Interface ScL 113, 544 (1986). 21. RheinboldL W. C., and Burkhardt, J. V.,ACM Trans. Math. Software 9, 215 ( 1983 ). 22. Slobozhanin, L. A., FluidMech. 15, 1 (1986). Soy. Res.

225

23. Tsaopoulos, J. A., Poslinski, A. J., and Ryan, M. E., J. FluidMech. 197, 523 (1988). 24. Ungar, L. H., and Brown, R. A., Philos. Trans. R. Soc. London Ser. A 306, 347 (1982). 25. Vogel, T., Pac. J. Math. 103, 205 (1982). 26. Vogel, T., SIAMJ. AppL Math. 47, 516 (1987). 27. Vogel, T., SIAMJ. Appl. Math. ,19, 1009 (1989). 28. Wang, T. G., Trinh, E. H., Croonquist, A. P., and Elleman, D. D., Phys. Rev. Lett. 56, 452 (1986). 29. Wente, H. C., Pac. J. Math. 88, 387 (1980). 30. Wente, H. C., Pac. J. Math. 88, 421 (1980). 31. Wente, H. C., Manuscr. Math. 39, 287 (1982).

Journal of Colloid and Interface Science, Vol. 146, No. 1, October 1, 1991