ATOMIZAÇÃO_phd4

7/30/2019 ATOMIZAO_phd4

1/15

Chapter 4

STATIC INTERFACE SHAPES

WITH VOLUME CONSTRAINT

I shall now proceed to lay before mathematicians another case, whichis still more remarkable, from the variety and singularity of thephenomena depending upon it, and from its being susceptible ofan equally accurate analysis; the case alluded to is that of capillary

attraction.

P. S. Laplace, Traite de Mecanique Celeste (1805)

As was discussed in chapter 3, test problems are needed to verify the VOF-

CSF program. One of the ones chosen is the problem of determining the shape of

the meniscus in a vertical cylindrical capillary tube. This tests the surface force

algorithm as well as the contact angle boundary condition. This problem is an

old and important one, predating the monumental work of Laplace (1805) and

Young (1805). It has applications in many fields and has become a prototype

problem explored by many workers, a recent historical perspective on which can

be found in Finn (1986). The capillary tube meniscus fits into the general class

of axisymmetric menisci that meet the axis of revolution (Princen, 1969). This

history reflects the fact that numerical analysis can provide the solution to the

problem to any desired accuracy, and by a variety of approaches. After some early

work of Bashforth and Adams (1883), Concus (1968) was the first to investigate

the problem systematically. More recently, Cuvelier and Schulkes (1990) solved the

problem including the volume constraint explicitly using a finite element method.

Results pertinent to this problem have also been obtained for the related

system of pendent drops and bubbles attached to horizontal supports (Finn,

56


2/15

57

1986). Lord Kelvin (1886) proposed a mechanical geometric integration and found

solutions to the pendent drop problem that had multiple bulges. Concus and Finn

(1979) concluded that not all surfaces have a simple single-valued projection onto

the reference plane, and parameterized the problem in terms of the arc length along

a vertical section of the interfacial surface. They then integrated the pendent drop

solution as an initial value problem but without inclusion of a volume constraint.

In this chapter we combine several concepts in producing a formulation

that allows the efficient, systematic evaluation of the solution to the capillary

problem corresponding to a given volume constraint. When implemented compu-

tationally, the method achieves similar accuracy with significantly less storage and

computational work for intermediate results compared to the finite-element meth-ods (Cuvelier and Schulkes, 1990) recently proposed for this problem, although

the latter are generalizable to two dimensions. The present method also permits

the bifurcation properties of the solution to be calculated in a natural way. The

numerical solution of the equilibrium meniscus shapes in a vertical cylinder with a

given volume of fluid is illustrated for all interesting values of the two parameters:

Bond number and contact angle.

4.1 Mathematical Description

The geometry of the problem is illustrated in Figure 4.1. A fluid of density

1 rests below a fluid of density 2, with the reference level at the bottom of

the tube. The surface is described either simply as a function z = u(r), or

parametrically for the more general case when z is not a single-valued function of

r, as will be presented in the next section. The governing equation is obtained

from a differential force balance on an element of axisymmetric interfacial area,

with a necessary condition for equilibrium being that the sum of the projection

of the forces in the z direction must vanish. For a derivation of this equation,

see appendix B. With gravitational acceleration, g, directed in the negative z

direction and the cylinder radius R used as the characteristic length, we have in


3/15

58

dimensionless form (Cuvelier and Schulkes, 1990):

1

rd

dr

rdu

dr1 +

du

dr

2

= Bo u + (4.1)

where r r/R is the dimensionless radial distance, u u/R is the dimensionless

meniscus height, Bo (1 2) g R2/ is the Bond number, expressing a measure

of the ratio of gravitational to surface tension forces, is the interfacial surface

tension, and is an integration constant yet to be determined. This equation (4.1)

is to be solved with symmetry and contact angle boundary conditions:

du

dr= 0 at r = 0 (4.2)

du

dr= cot at r = 1 (4.3)

where is the contact angle with the wall measured within fluid 2. In addition,

the volume of fluid 2 is to be constrained with the reference level at z = 0:

Vc =10 2r

u

dr

(4.4)

The fluid 2 is wetting or non-wetting depending on whether the contact

angle is in the range 0 < 2

or 2

< respectively. The case of a non-

wetting fluid can be obtained from that of a wetting fluid by replacing by

and u by 2Vc u. Thus, only a wetting fluid will be considered. The trivial

solution for = 2

is the flat interface, u Vc .

Typically, equation (4.1) can be solved numerically as a boundary value

problem, subject to the conditions (4.2) and (4.3) with the value of determined

by imposing (4.4) as an additional constraint (Cuvelier and Schulkes, 1990).

However, this approach requires a good initial guess for the interface profile,

due to the nonlinearity of (4.1). This is typically provided by an asymptotic


4/15

59

z

r

g

z = u(s)

r = r(s)

1

2

Fluid 2

Fluid 1

Figure 4.1: Schematic of the meniscus in a vertical cylindrical tube.

solution valid for some range of parameters and continuous in the parameter space.

What we propose in the present chapter is a much simpler numerical procedure

following a rearrangement of the original equation (4.1) dictated from an a priori

implementation of the constraint (4.4) and boundary conditions (4.2) and (4.3).

As can be seen by inspection of equation (4.1), is a constant equal to the

dimensionless pressure difference po (p

2o p

1o) = 2/R

o across the interface

minus the hydrostatic head, Bo uo, where u

o is the meniscus height and R

o is

the dimensionless radius of curvature at the centerline. It was Laplace (1805)


5/15

60

who found a first integral of equation (4.1) and related the unknown, , to the

boundary conditions. If we multiply (4.1) by r and integrate from zero to one we

obtain an expression for the unknown (Cuvelier and Schulkes, 1990):

= 2 cos Bo Vc

(4.5)

This incorporation of (4.3) and (4.4) directly into equation (4.1) eliminates

either one of (4.3) or (4.4) from further consideration in the solution and forms

the basis for our computational approach. This integration is possible only for

Neumann boundary conditions, such as (4.2) and (4.3), and cannot be carried out

for Dirichlet conditions, such as a given contact line.

The much simplified computational strategy, then, involves simply solution

of (4.1), using (4.5) for , by a forward integration using (4.2) and u (0) = u0

as

initial conditions. The proper value ofu0

is found iteratively by imposing either of

(4.3) or (4.4). The key advantage as compared to the previous scheme is the need,

for a given value of the parameters (Vc, , B o), of the approximate value of only

one constant, u0

, as opposed to a full interface profile, for the iterative scheme to

converge. In addition, powerful adaptive ODE solvers can be used to implement

the numerical method effortlessly.It is of interest to calculate the total energy of the system to decide which

states are thermodynamically favored. Given multiple solutions at the same value

of (Vc, , B o), the configuration with the lowest total energy will be the most

favored. Referring to Figure 4.1, if we choose z = 0 as the reference level, the

total energy is given in dimensionless form by the sum of the potential energy,

surface, and contact Helmholtz free energy:

Et =10

Bou2r dr +10

2r sec dr cos [1 + 2u (1)] (4.6)

where Et is scaled by R2, and is the angle of the tangent to the interface

with the horizontal. Use has been made of Youngs (1805) equation to derive the


6/15

61

last term in (4.6):

cos = 1s 2s (4.7)

where 1s and

2s are the interfacial surface tensions between the solid wall and

fluid 1 and fluid 2, respectively. It is interesting that, alternatively, the Young-

Laplace equation (4.1) can be derived from the energy equation (4.6) and the

volume constraint (4.4) using variational techniques (Finn, 1986).

4.2 Numerical Solution

The problem then is to solve (4.1) subject to (4.2)-(4.4). Before doing so,

however, we first reformulate the equations to allow for the possibility that the

meniscus height, u, is not a simple function of radius, r, which occurs when

the interface folds back on itself (Concus and Finn, 1979). If we introduce the

dimensionless arc length s and use the chain rule, we obtain a system of first-order

coupled nonlinear differential equations, with independent variable s. Equation

(4.1) becomes:

du

ds= sin , u(0) = uo (4.8)

with the variation of r and given by:

dr

ds= cos , r(0) = 0, r(sf) = 1 (4.9)

d

ds= Bo u + 2 cos

Bo Vc

sin

r, (0) = 0 (4.10)

where sf is the total dimensionless meniscus arc length from the centerline to the

wall. The integral (4.4) for the volume constraint becomes:

dV

ds= 2ru cos , V (0) = 0, V

sf

= Vc (4.11)

The equations (4.1)-(4.4) thus reduce to a set of four ordinary differential

equations (4.8)-(4.11) with four initial conditions and two boundary conditions.


7/15

62

The (decoupled) equation for total energy (4.6) becomes, in terms of arc length:

d

ds= Bo ur cos + 2r, (0) = 0 (4.12)

Et =

sf cos

1 + 2u

sf

(4.13)

The numerical integration of these equations is limited by the fact that the

centerline height, uo, and the total arc length, s

f, are not known a priori. These

unknowns are found by using Newtons method to satisfy the volume constraint

(4.11) and the wall constraint (4.9), allowing the problem solution to be obtained

for a given (Vc, , B o). The volume and wall constraint equations are, with Vc and

fixed:

f1

uo, s

f; Bo

= V

uo, s

f; Bo Vc = 0 (4.14)

f2

uo, s

f; Bo

= r

uo, s

f; Bo 1 = 0 (4.15)

where uo, s

f are used as primary variables for which equations (4.14)-(4.15) are

solved, with Bo considered as a continuation parameter. In addition, because the

nonlinearity of the equations can lead to turning points and solution bifurcations,

a parametric arc length continuation method (Keller, 1977) was used to followsolution branches through turning points. For continuation through turning points

with respect to Bo, the system equations (4.14)-(4.15) are augmented by:

f3

uo, s

f; Bo,

= 2

uo uRo

2

sf sRf

2

Bo BoR2

= 0 (4.16)

where is the parametric arc length distance between a reference point xR =

uRo , s

Rf ; Bo

R

on the branch and another point x = u

o, s

f; Bo also belongingto the same solution family. The problem then generalizes into that of solving

the system of equations f (x) = 0 with Newtons method starting with a known

solution xo, considering f (f1, f2, f3) a function of an expanded set of variables

x

uo, sf; Bo

and the parametric arclength, , as the new parameter.


8/15

63

The expanded Jacobian matrix in the Newtons method has the advantage

that it is nonsingular at turning points wheredBo

d= 0. Critical points can be

identified by the vanishing of the determinant of the Jacobian J det

(V, r)

uo, s

f

.To use Newtons method the nine partial derivatives

f

xwere calculated by

differentiating equations (4.8) to (4.11) analytically (see appendix C), resulting in

eight additional ordinary differential equations. These equations contain some

terms that are indeterminate as r, s 0, and so lHopitals rule was used

on the indeterminate terms in a manner similar to that used by Siekmann et

al. (1981). The equations (4.8)-(4.12), and the ordinary differential equations for

the partial derivatives, were integrated with a commercially available, variable

step size, fifth- and sixth-order Runge-Kutta-Verner technique (IMSL, 1987).

The continuation method was started from the exact solution at Bo = 0 (see

below), with the parametric step size adaptively determined by the speed of

convergence of the Newtons method. The reference point, xR, was updated every

five steps to maintain a stable calculation. Convergence was based on the criterion

f 107.

At the zero gravity limit, Bo = 0, inspection of equation (4.1) reveals that

the total curvature is constant, therefore the solution is the part of the surface

of a sphere that satisfies the symmetry and contact angle boundary conditions

(4.2)-(4.3). This case serves as a check on the numerical solution, and as a known

starting point for the continuation method described above. The solution can be

expressed in closed form in terms of the arc length for the range 0 s sf by

extending the results of Concus (1968):

= s cos (4.17)

u = uo +(1 cos )

cos (4.18)


9/15

64

r =sin

cos (4.19)

where the final arc length is given by:

sf =2 cos

(4.20)

and (4.11), (4.18), and (4.19) are used to obtain the centerline height:

uo =Vc

1

cos +

2

3

1 sin 3

cos 3(4.21)

The total energy from (4.12), (4.13) for this case becomes:

E

t =2 (1 sin )

cos 2 cos

1 + 2

u

o +1 sin

cos

(4.22)

4.3 Results

A range of results obtained using the above method are shown in Figures 4.2

to 4.4. Figure 4.2 depicts meniscus centerline height, uo, meniscus arc length, s

f,

the determinant of the Jacobian, J, and total energy, Et , vs. positive Bond

number, Bo, for volume Vc = , and for contact angles = 0, 30, and 60.

Figure 4.2 for positive Bond numbers shows the same monotonic trend in all

variables from the spherical solution at Bo = 0 to the asymptotically flat solution

u Vc

as Bo . The determinant of the Jacobian, J, remains positive, except

for = 0 for which J = 0 for all Bond numbers. This can be seen by inspection of

equations (4.9) and (4.11) evaluated at the wall where = 2

. In the absence of

turning points for this singular case, a theoretical infinity of solutions satisfying the

equation and boundary conditions exists. We interpret these solutions physically

to be very thin films of varying heights wetting the wall of the tube.

Figure 4.3 shows that for the same variables as Figure 4.2 , but for negative

Bo, turning points are present. Curves in this study were terminated when the

meniscus just touched the wall of the tube and further changes in the Bond number

in the same direction would cause it to bulge out of the tube. The Figure 4.3,


10/15

65

(A) (B)

* *

(C) (D)

*

Figure 4.2: Meniscus centerline height, uo, meniscus arc length, s

f, determinantof Jacobian, J, and total energy, Et , vs. positive Bond number,Bo, for volume, Vc = . Dashed curve for = 0, dotted curve for = 30, and solid curve for = 60.


11/15

66

= 0, curves terminate at a lowest Bond number of about Bo = 0.842. As

can be seen from Figure 4.4, which depicts the meniscus profile as height, u,

vs. radius, r, for contact angle = 0, 30, and 60 for various Bond numbers,

Bo = 0.842 is where the meniscus starts to develop a bulge outside the tube,

and this is also a turning point. This limiting Bond number is precisely the same

as that at which the meridional curvature of u,d

ds, in equation (4.10), becomes

zero at the wall, in agreement with Concus (1968).

The curves in Figure 4.3 for = 30 show a turning point at Bo = 2.506,

with the solution branch terminated beyond the turning point at Bo = 0.579,

where again the meniscus bulge attempts to leave the tube (Figure 4.4b). However,

the meridional curvature becomes zero at the wall earlier, on the main solutionbranch at Bo = 2.03, in agreement with Concus (1968). He concluded that a

stable equilibrium surface is no longer possible at lower Bond number in agreement

with the previous variational analysis by Reynolds and Satterlee (1966). A more

rigorous approach to static stability analysis as was used by Wente (1980) for the

pendent drop could be applied to the cylinder. It can be seen also that the total

energy reaches a minimum value at the turning point, at which it then begins to

increase. This would indicate that for a given Bond number, the configurations in

the main family are favored thermodynamically relative to those solutions past the

turning point. However, no claim is being made here about dynamic stability, as

this would require inclusion of velocity variations. For example, Maxwell (1890),

performing an inviscid dynamic analysis, found that a plane horizontal interface =

2

on a circular boundary cannot be formed in stable equilibrium when,

in our nomenclature, Bo < 3.39 . . ., this being the first root squared of the

derivative of the J1 Bessel function.

The curves in Figure 4.3 for = 60 reveal similar behavior to the = 30

curve except that no terminal Bond number is reached and the solution branch was

arbitrarily terminated at uo = 3. Even though u

o is allowed to go below zero, the

bottom of the container, this offers no conceptual difficulty since any other volume


12/15

67

(A) (B)

*

*

(C) (D)

*

Figure 4.3: Meniscus centerline height, uo, meniscus arc length, s

f, determinantof Jacobian, J, and total energy, Et , vs. negative Bond number,Bo, for volume, Vc = . Dashed curve for = 0, dotted curve for = 30, and solid curve for = 60.


13/15

68

could be specified which would simply translate the meniscus vertically (Finn,

1986; Concus, 1968). The shapes with bulges in Figure 4.4b-4.4c are known as

unduloids. This unduloid behavior is identical to the initial shapes of the pendent

drops presented by Lord Kelvin (1886), and Concus and Finn (1979).

(A) (B) (C)

*

*

*

*

*

*

Figure 4.4: Meniscus height, u, vs. radius, r, for volume Vc = , and forcontact angle = 0, 30, and 60 deg. Bond numbers for contactangle = 0 in order from the flattest to the most arched areBo = 100, 10, 0,0.842. Bond numbers for contact angle = 30

are Bo = 100, 10, 0,2.506,0.579. Bond numbers for contact angle = 60 are Bo = 100, 10, 0, -5.886, -1.643, -2.038, -2.214, -2.319.

No solution bifurcations have been found in all cases investigated (with the

exception of the singular case, = 0). Each time the determinant of the Jacobian,

J, passed through zero, a turning point was found. The total energy, Et , again

shows the configurations to be increasingly unfavored thermodynamically as each


14/15

69

turning point is passed and the interface develops increasingly inflection points

and folds in agreement with intuition. However, other solutions were found off

the main physical solution branch. These exhibit looping self intersections that

are known as nodoids and are illustrated in Figure 4.17 of Finn (1986) for the

pendent drop, but they have been eliminated from further consideration as being

not physically realizable (Lord Kelvin, 1886).

4.4 Conclusions

The solution of the equilibrium shape of the meniscus in a cylindrical

tube of a given fluid volume has been evaluated for various values of Bond

number and contact angle. A formulation has been developed that explicitly

incorporates the volume constraint and that allows accurate and efficient numerical

implementation. The results have been extended to negative Bond numbers not

previously reported for the capillary problem. No solution bifurcations have been

found for the parameters studied. Results confirm the existence of turning points

and thermodynamically unfavored menisci with multiple bulges (Finn, 1986; Lord

Kelvin, 1886; Concus and Finn, 1979) as with the pendent drop.

The advantages of the solution method presented in this chapter vs. pre-

vious formulations are several. First, because of the multivalued shape of the

menisci, the arc length parameterization in s allows these complex shapes to be

calculated with a standard ordinary differential equation solver in an accurate,

efficient manner. Second, the closed-form first integration of the Young-Laplace

equation (4.1) eliminates the unknown and the contact angle boundary condition

(4.3) in favor of the more physically intuitive unknown centerline height. Third,

for a given Vc, , and Bo, the entire solution can be reconstructed knowing only

the values ofu

o, and s

f, presented in Figures 4.2 and 4.3, eliminating the need forstorage of the entire meniscus profile. Fourth, by expanding the system with the

parametric continuation parameter, , the solution can be continued systemat-

ically through turning points in the expanded x

uo, s

f; Bo

space and while


15/15

70

still explicitly incorporating the volume constraint. Finally, the method presented

in this chapter can be generalized to other axisymmetric systems such as spheres,

cones, and ellipsoids with or without axial rotation.

Documents

ATOMIZAÇÃO_phd4