Theory of Meniscus Shape In Film Flows. A Synthesis Brian G. Higgins, William J. Silliman, Robert A. Brown, and L. E. Scriven' Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, Minnesota 55455
Free surface film flows over solids are dominated by the boundary conditions of capillarity in concert with effects described by the components of the Navier-Stokes equation parallel and perpendicular to the solid surface. Strategies for simplifying, transforming, and combining the basic equations and solving for meniscus shape are presented: integral balances, differential approximations, and perturbation techniques. Examples are drawn from the literature on surface leveling, dip coating, and rimming flow. Assumptions inherent in the various analyses are brought out, their interrelationships are established, and new directions are indicated.
Introduction Small-scale free surface flows in which surface tension and viscous effects are important are central to fluid contacting operations in chemical engineering. Film flows are common in the manufacture of photographic films and a myriad of other coatings. They are also employed in various heat and mass transfer operations and figured prominently in the pioneering studies of the fundamentals of convective transport processes without and with chemical reaction (Johnstone et al., 1941;Emmert and Pigford, 1954;Fulford, 1964;Sherwood et al., 1975). While relatively simple film flows have not been neglected in basic fluid mechanics, the past score of years has seen a n upsurge of research devoted to analyzing and predicting more complicated film flows, and viscocapillary motions more generally. These flows are of course governed by the continuity and momentum principles, which yield a set of partial differential equations on the flow domain and equations of one lower dimension on the domain boundaries. The boundary conditions are usually nonlinear owing to the shape of the meniscus which, being unknown a priori except in limiting cases, also adds the complication of a free-boundary problem. The momentum equations within the flow domain may of course be made nonlinear by inertial effects and by non-Newtonian rheological behavior. Complete solution of the equation set is not possible except in quite special limiting cases. A burgeoning number of approximations have appeared in the literature, many of them addressed to the engineering need for estimates of the location of the meniscus in space and time. The various approximation schemes have not been compared. What has been lacking is a basis for analyzing themand for generating still other schemes which might be superior, for one type of flow or another. The situation is particularly interesting because the number of ways of recasting and combining the equations of the basic set, for purposes of approximation or otherwise, is large. There is then a basic issue of strategy in attacking a given viscocapillary flow problem. We present three primary strategies for arriving at approximate solutions: integral balances, differential approximations, and perturbation analysis. We show how the approaches used in the literature are connected and we bring out the relation to numerical solution of the full equation set. For simplicity most of the discussion is restricted to two-dimensional flow of Newtonian fluid, which suffices to develop the main points. Non-Newtonian effects alter the equations but not most of the strategies. Two-dimensional flows are what have been considered in the literature. Three flows are drawn from the literature to illustrate the interconnections of the solution strategies. They are dip coating, surface leveling, and rimming flow.
Underlying Equations The governing equations within the liquid regarded as incompressible are
v*v=o
(2)
where p is the density, v is the velocity, g the body force, p the pressure (relative to that in the adjacent gas phase) and T the viscous stress tensor. At the free surface the appropriate form of the momentum principle is np - n
- + n2%u + V,u = 0 T
-
(3)
where n is the field of unit normals to the surface, Yf = 7 ' , n/2 is the mean curvature, u is the surface tension, and V, is the surface gradient operator. Equation 3 applies when the viscosity and density of the adjacent gas phase are insensible. For the two-dimensional flows considered here the location of the free surface can be represented as h = h ( x , t )as shown in Figure 1. Then
n = (-hxi
+ j ) / & ; t = (i + h,j)/&
(4)
m,
where t is the field of unit tangents, G 3 and h, ahlax. The position of the free surface must satisfy the kinematic condition ah -1 _ =n.v (5) G at On parts of the boundary of the flow domain which are rigid the adherence or no-slip condition and impenetrability condition applies v = vs
where v, is the velocity of the solid surface. Except in the rare closed system like rimming flow there are parts of the boundary across which there is inflow and outflow of liquid. In what is generally an approximation one supposes there are upstream and downstream locations at which the velocity field on the boundary of the flow domain is known, Le., v = vb. The surface traction condition (3) can be resolved into two components, respectively normal and tangential to the free boundary. The normal component is greatly simplified by substituting the tangential component ( 7 ) into it. Thus
The mean curvature is 7-f = ( ~ 3 ~ h / a x ~ ) /Equations 2 & ~ / ~ .6 and 7 together with 4 and 5 and the inflow and outflow conditions Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977 393
n NUMERICAL SOLUTION
OF
PERTURBAT] ON ANALYSIS
FULL SET
ilnrden 8 Harlow ( 1 9 7 0 a . b ) O P ~8 S c r i v e n (1976)
b t h e m n S Homiy (1976) Benne, (1966) K r d P t i S G w e n (1570)
X
BASIC EQUATION SETS W I E N T U H , C O N T I N U I T Y A\D K I N E M A T I C P R I N C I P L E S FOP, L l O U l O AN3 VEf1IS:LS equations ( I ) (7)
-
constitute a full set of boundary conditions for free-surface flow of a liquid governed by (1) and (2). This boundary value problem is especially difficult because the location of the meniscus and hence, in mathematical terms, the domain of the problem is not known a priori. The meniscus location and shape are coupled with the velocity and pressure fields through (5), (6), and (7): thus the meniscus shape reflects the interplay of capillary, pressure, viscous, inertial, and gravitational forces. We note that when the meniscus is of main interest, eq 5-7 which govern it directly can be regarded as the problem and eq 1 and 2 which describe liquid behavior are subsidiary. The roles of the higher-dimension domain and its lower-dimension boundary are reversed. Some solution strategies are best appreciated from this novel point of view, as we bring out below. Solution Strategies If details of the velocity field within the liquid film are not required, eq 1 and 2 can be integrated across the film, with the velocity and pressure fields approximated by some more or less reasonable functions, and the result used with the normal and tangential stress conditions (6) and (7) to generate a shape equation. The change of shape with time is described by the kinematic condition (5). Together these two evolution equations govern the location of the meniscus in space and time. Once the evolution equations are derived the higher-dimension flow problem becomes secondary. Owing especially to the nonlinearity of the curvature expression the evolution equations are generally nonlinear. Sometimes qualitative information about the nature of the velocity and pressure fields and about the shape of the meniscus is available, yet more details of the former are desired along with accurate estimates of meniscus shape. Then an appropriate strategy is to size the terms, remove those that seem small from the full equation set, and in some cases to hold constant a variable which is nearly so. This sort of simplification is most successful when it reduces (1)and (2) and, in the extreme, (5)-(7) to versions that can be solved in closed form, which is of course most likely when all nonlinear terms are either removed or linearized. The result is again meniscus evolution equations which are nonlinear or, in the extreme of small departures from a highly symmetric shape, linear. Sometimes perturbation techniques can be used to generate approximations more systematically in the form of a set of linear problems which can be solved sequentially for increasingly detailed information about departures from a tractable, and hence simple, base flow. Suitable small parameters have to be identified in which to expand the variables, velocity, pressure, and meniscus location, and the results are applicable only when the chosen parameters are small enough. As with less systematic approximations the range of validity usually can be established only by comparison with experiment or with suitable numerical solutions of the full equation set. On the other hand, a perturbation analysis can have great value by providing a rigorous analytical solution for a limiting case against which a numerical simulation scheme can be tested, as in an illustration below. 394
Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
I DIFFERENT I A t APPEOX I MAT IONS Linear
'ionlinear
Anrhur (1973)
Deiber B C e m
(1976)
Orchard (1962)
Van R o r i u m ( 1 9 5 8 )
D e i b e r 8 Cerro (1576) Landau & Lev'ch (1942)
I
I
INTEGRAL BALANCES
I
-
S p i e r s et dl (19741 Williamson (1972)
Wehausen
L
taitnnp iiqfin\
Figure 2. Interrelation of solution strategies. Large arrows point to alternative approximations; small arrows indicate means of testing approximations.
Figure 3. Rimming flow in a horizontal cylinder of radius R spinning about its axis at angular velocity w .
Numerical solution of many two-dimensional viscocapillary flows is feasible today. The greater the departure of the meniscus from a highly symmetric shape the more expensive are the calculations in time and memory, because of the nonlinearity of the equation set and the variable location of the boundary. Even with the finite element technique, which is well-suited to film-flow problems, the available high-speed memory of any but the largest computers today is easily filled and the calculations are not cheap. Nevertheless numerical simulation can be substantially less costly than experimentation, whether for understanding the mechanics of a flow, or gathering design and troubleshooting data, or validating approximate solutions intended for similar purposes. The interrelation of these four strategies is shown in Figure 2. The researches cited are the main contributions regarding the three examples we discuss here. The first is surface leveling, the damping of surface irregularities on a nearly flat film: see Figure 1. The second is rimming flow, the motion of a film carried on the inside surface of a horizontal cylindrical vessel rotating about its axis: see Figure 3. The third is dipcoating flow, the continuous deposition of a film on a solid substrate emerging vertically upward from a large bath: see Figure 4. Although Newtonian liquids are considered here, flows are of course non-Newtonian in many applications. Non-Newtonian effects were considered in surface leveling (see review article by Quach, 1973), in dip-coating (e.g., Gutfinger and Tallmadge, 1965; Groenveld, 1970; and Spiers et al., 1975) and in rimming flow by Rao and Throne (1972). Approximate Differential Equation Sets To estimate relative magnitudes of the terms it is desirable to scale the variables and put the equations into dimensionless
Table I. Physical Parameters and Dimensionless Groups Symbol
Figure 4. Dip-coatingflow: the coating of a liquid film onto a vertical surface continuously withdrawn from a bath. The dashed line represents the static meniscus profile (i.e., U = 0). form. Appropriate units and dimensionless groups are given in Table I. In terms of dimensionless variables, which are denoted by overbars, (1)and (2) for a Newtonian liquid become
Quantity Acceleration of gravity Datum thickness x -Length scale y-Length scale Pressure unit Volumetric flow rate Time scale of flow x -Velocity unit y-Velocity unit Deviation of thickness from datum Kinematic viscosity Density Surface tension = capillary number = Euler number = Froude number = Reynolds number = Strouhal number = Weber number = Aspect ratio = Velocity ratio = Characteristic slope
and, a t a datum location y = ho for the meniscus
ah -=u at au - + - =aou ay ax aii - + - - =p oali
ax
(10)
ajj
and (5)-(7), the equations a t the meniscus come
y
=
E(?,?),
be-
Equations 12 and 13 are written here for the case in which surface tension is uniform (excluded therefore are practically important methods of controlling surface leveling, for example, with surfactants or mixed solvent systems). When as in the leveling of a very nearly horizontal film the velocities relative to the substrate solid are known to be small and so are the meniscus slopes, the parameters N Rand ~ e are both of low order and the equations are well approximated by linear ones (provided a and p are of order one). These are written here in the dimensional form preferred in the literature (Lamb, 1945) (14)
Because there is no net convection overall the upstream and downstream conditions in x degenerate into simple boundary conditions. If the film is effectively unbounded these may be taken as periodic boundary conditions, consistent with Fourier analysis of the meniscus shape (cf. Lamb, 1945). The advantage of a linear system like this one is of course that it can be solved in terms of tabulated functions by which the dependence of film flow and meniscus shape on the dimensionless parameters can be examined relatively easily. Starting from (14)-(19), Wehausen and Laitone (1960) used a technique of Basset (1888)and Lamb (1945) to develop an expression for w ( h ) , the time scale of leveling, which in the linear approximation is naturally an exponential process; i.e., a meniscus variation of wavelength X decays according to h(X) = exp[w(X)t]. The dispersion equation governing w(X) is transcendental. Wehausen and Laitone (1960) found a limiting form for o(X) in the low viscosity limit: w is directly proportional to viscosity and is complex, which indicates oscillatory approach to flatness. Orchard (1962) and Anshus (1973) also attacked the linearized leveling problem but neglected the local acceleration terms &/at and dulat in (14) and (15). They found w to be inversely proportional to viscosity and real, which indicates monotonic approach to flatness. This is the high viscosity limit (and also the thin film limit), in which one expects Nst to grow large enough to justify the Orchard and Anshus analyses. This is borne out by an example discussed by Degani and Gutfinger (1974). The range of validity of the much cited surface leveling formulas of Orchard and Anshus is clarified by the more systematic linear apInd. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
395
proximation, (14)-(19). When the nonlinear effects excluded by that approximation become significant remains to be assessed. Van Rossum’s (1958) analysis of dip coating provides another example of a linearized analysis. Looking at the virtually flat film that may exist far above the coating bath, he in effect supposed that capillarity and inertia are negligible compared to viscosity, the flow is very nearly rectilinear, and therefore most of the terms in (8)-(13) can be dropped. In dimensional form the only survivors are pg and pa2uldy2 in (8); the boundary conditions are u = UOa t y = 0 and duldy = 0 at y = h; and the velocity profile in this flat-meniscus approximation is a familiar one
Pressure can then be eliminated from the normal stress condition by means of the component momentum eq 8 and 9 in the neighborhood of the meniscus. The result is simplified by using the continuity equation (10) and tangential stress condition (13). The technique of using the directional derivative to eliminate the pressure from the normal stress boundary condition is not novel. It has been used by Krantz and Goren (1970) and Atherton and Homsy (1976) for studying the stability of falling liquid films. The differential equation of meniscus shape is (Higgins and Scriven, 1977)
-
u = Uo - g(hy y 2 / 2 ) / ~ (20) Film thickness h is related to volumetric flow rate q upward past a fixed station by
To3- 3To t 3Q = 0; T o h(pgIpUo)1/2; Q ~(glvUo3)l’2 (21) Q is the last vestige of the upstream and downstream conditions. It is not enough, because according to (21) the maximum value of Q is 2/3, corresponding to To = 1,and for 0 < Q < 2/3 there are two values of To and thus of film thickness h. Van Rossum and others specuiated that the upstream and downstream conditions could be replaced by one of the elusive “maximum principles” but his experiments showed that Q does not get close to 213, and when the assumptions of Van Rossum’s analysis are relaxed, analytical results can no longer The companion equation of shape evolution is the kinematic be obtained and strategies other than linearization are needcondition ed. Two analyses of rimming flow in which Deiber and Cerro (1976) neglect capillarity afford further examples. Terms are eliminated from (8)-(13) by order-of-magnitude arguments These equations are exact. The dimensionless parameters are until equations of boundary-layer type survive for film flow inside the rotating cylinder. In the limit of low N R and ~ N F ~ defined in Table I and all quantities are evaluated at the meniscus = E(?,?). Approximations are created by inexact (for the rotating, nearly circular system but defined analoexpressions for 51 and E, by outright omission of terms, and by gously to the numbers in Table 1) these equations become errors incurred in constructing solutions of (23) and (24). linear and lead to a transcendental, solvable equation for the location of the free surface. In the limit of large N R , and N F ~ A number of authors have adopted ad hoc simplifications of (23) as equations of static meniscus shape. Examples are the boundary-layer equations again become linear and can be given in Table 11. The original derivations by Landau and solved for a velocity distribution which substituted in the Levich (1942), White and Tallmadge (1965), Williamson integrated continuity equation yields the shape of the free (1972), Spiers et al. (1974), and Lee and Tallmadge (1975, surface. However, Deiber and Cerro (1976) find that much of 1976) were circuitous and incorporated assumptions along the the parameter range of interest falls outside these two linway. In some cases the resulting formula is inconsistent beearizations. cause it retains terms of the same order as others that are Deiber and Cerro’s (1976) equations of boundary-layer type omitted. The inconsistency of Williamson (1972) and Spiers exemplify useful nonlinear approximations that have been et al. (1974) in replacing (9) by @/a7 = 0 while retaining dEld7 developed. These authors finesse the free boundary compliin (12) was first pointed out by Esmail and Hummel(1975a), cation by transforming to Protean coordinates (Duda and who attempted to rectify the matter by using an integral apVrentas, 1967) in which the stream function becomes an inproach (see next section) to include further terms. dependent variable. Since the free surface is itself a streamline The rigorously derived eq 23 accommodates flows much less the transformed domain is enclosed in a known, fixed nearly one-dimensional than treated heretofore, and in conboundary (provided the flow is closed, for otherwise the upcert with (24) covers transient film flows as well. This pair of stream and downstream conditions interfere where streamevolutionary equations requires information about the liquid wise viscous stress is appreciable). However, the transforflow only locally, immediately inside the meniscus. In this way mation to Protean coordinates can be made only for flows with the higher dimension flow field is replaced by a neighborhood no interior recirculation; i.e., all streamlines must encircle the of virtually the same dimension as the meniscus itself. axis of rotation. Moreover, the resulting equations are intractable analytically and so were solved by the finite differIntegral Balances ence method by Deiber and Cerro (1976). When approximation is necessary and accurate information Another strategy of approximation exists in the dip-coating is lacking about the flow field just inside the meniscus, there literature. It consists of assuming a u-field and, presumably, may be advantage in evaluating momentum exchange with a consistent v-field, using the x - and y-momentum equations the meniscus by means of a balance over a differential slice to evaluate the pressure exerted by the liquid on the free of the film (differential thickness in the x direction). This is boundary, and thereby deriving a shape equation for the particularly appropriate if the velocity field itself is required meniscus. The key is the directional derivative, along the to satisfy the momentum equation only on the average over meniscus, of the normal stress boundary condition (12), here the film depth, as in the integral method or moment method written 9 (x,y,t) = 0 396 Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
Table 11. ADDroximate Eauations of Meniscus ShaDea Coupling strategy
Source
Neighborhood Equation 23
Differential or integrodifferential equation d d2H m 2 -3/2 1 -1 71 + 32/3N~a2/3 dx dX 3H2 (E) aY2 Y = I
[
]
(E)
+
(z)y=l
+ 32/3N~a2/3
(g)(E) *[
+ 12Nca4/3
dH d2H 1 - 32'3N
2/3
Ca
-
v - 3)
__ a2
au2
ax2
Landau and Levich (1942) Spiers et al. (1975)
Equation 25 Integral-ymomentum
+1 ( E ) + & Nca2/3[ Y=I
d3H
dX3 3H2 a Y 2
[ [
d d2H dX dX
dH
- 7 1 + 32/3N~a2/3 (E) d2 +--Nca2/3 3lI3 dX2
+
dX3 3H2
(z)
To2= +L ( 3 )E) - H axay d x Y = I 3
] ] + 3 1~ (2) 2 -3/2
Y=o
+ [J (z
dY
32/3NCa2/3-
ax
J' UVdY] - -f-= T2
---[HZ 1 d2 Nwe dX2
Integral-xEsmail and momentum ~~~~~l (1975a)
2
Y=I
1
0
1(g) - (g) ] aY
Y=I
ay
Y=O
(197513)
of finding approximate solutions to the boundary layer equations. The exchange of either y momentum or x momentum can be found by integrating the respective equation over the film thickness. Either choice leads to an integrodifferential shape equation. In the first case the y-momentum equation (9) is integrated with respect to y and differentiated with respect to x . The term ap/ax at y = 0 is eliminated with the x-momentum equation (8) evaluated there; the term dplax at y = h is eliminated with the directional derivative of the normal stress condition (12) there. This brings in splay at y = h, which is eliminated with (9) evaluated there. The result is simplified
by means of the adherence condition, continuity equation (lo), kinematic and tangential stress conditions (11) and (13). Thus
sh + sh
pa (A ax at a2
udy
A
+ 1.1 1 ax
h
2 ax
o uudy]
au (ay + $) dy + ( $ ) y = o
+ pgz + pgy
(25)
As an equation for the meniscus profile, this is likely to be the Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
397
more accurate the more closely u and u come to satisfying eq 8-11 and 13, of which the continuity equation and x-momentum equation near y = 0 are probably the most important although this point has not been settled. The continuity equation can be guaranteed of course by approximating the stream function rather than the velocity, and all of the equations here are easily recast in terms of stream function. When the assumed velocity field is steady, eq 25 reduces to the nonlinear, ordinary differential equation for h ( x ) which is shown in dimensionless form in Table 11. Essentially this time-independent specialization was first derived by Gardner and Adebiyi (1974) for the film flow around a two-dimensional bubble migrating steadily along an inclined channel. In the second case the x-momentum equation (8) is integrated with respect t o y to give
1
p aa t J h0
udy
+a ax o Jh
u2dy]
(z
- - :)y=h
=-Jh$dy+k
-p
(” - ”) ay ax
y=o
t pg,h
This formula can in principle be differentiated with respect to x and the result inserted in (26). The outcome is cumbersome until a specific velocity distribution is substituted and each of the terms is evaluated or neglected. This was done by Esmail and Hummel(1975a) for a steady film flow with slight meniscus inclination, for which the equations are
while the normal stress boundary condition is approximated as au d2h p = 2r-U T (29) ay dx The result is shown in Table 11. Esmail and Hummel(1975b) have extended this method to include inertia and extra viscous effects by integrating the full x -momentum equation, while still approximating the y -momentum equation as (30) ay-’ay2 and the normal stress condition by (29). It is plain that integrating the x-momentum equation does not ordinarily lead to the same meniscus equation as integrating the y-momentum equation, nor to the locally based eq 23. There is a velocity field, however, in which the three 398
Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
u=
uo + 3
($-5)
2
- hY)
(31)
where q is the volumetric flow rate. The question of which integral scheme provides greater accuracy for comparable effort is not settled, although the y -momentum integral obviously has the advantage of relative simplicity. Boundary Conditions No matter what scheme is used to derive a differential equation of meniscus shape or shape evolution, integrating the equation requires upstream and downstream boundary conditions. In principle, these depend on any meniscus attachment to solid boundaries and on the upstream and downstream flow conditions. In dip coating it has heretofore been customary to restrict consideration to films that asymptotically approach uniform thickness downstream
(26)
Unlike the integration of the y-momentum equation, the pressure cannot be eliminated directly from (26) and in its present form eq 26 requires assuming a pressure field as well as a velocity field. However, the pressure field can be evaluated in terms of the velocity field by integrating the y-momentum equation from the meniscus a t y = h(x,t) to the variable location y . This relates p(x,y) to the pressure just inside the meniscus and then the normal stress boundary condition connects the latter to the ambient pressure outside the meniscus. Thus with due precaution in interchanging orders of differentiation and integration the pressure field is
9 - a2u
meniscus equations are actually the same. This is so when the flow is assumed to be one-dimensional (Landau and Levich, 1942; White and Tallmadge, 1965) and the velocity profile is of the form
The problem is to predict that downstream thickness, ho. If uniform thickness is indeed approached, it is controlled, through mass conservation, by the upstream region where liquid is entrained into the forming film. This region includes the boundary layer adjacent to the rising solid substrate, some neighboring portion of the coating bath, and the outward extending meniscus (Figure 4). At low substrate velocities this meniscus and the liquid beneath it are almost motionless. Landau and Levich (1942) introduced a novel subterfuge of matching the dynamic meniscus downstream where the difference between viscous and gravitational forces is of the same order as the capillary pressure caused by meniscus curvature, with a static meniscus upstream that reaches far away from the moving substrate and is shaped by gravity and capillary pressure only. This would be successful if there were not an interval in which the meniscus responds to all three forces, viz., viscous, gravitational, and capillary. Perturbation analysis alone supplies a rigorous basis for matching the two meniscus regions, as mentioned in the next section. Recent authors (White and Tallmadge, 1965, Spiers et al., 1975) first linearize the meniscus shape equation and then apply the Landau and Levich matching condition, while Esmail and Hummel (1975b) use a different approach. Assuming the meniscus slope is everywhere small, they integrate a corresponding shape equation upstream to the presumed location of a stagnation line on the meniscus, beyond which they take the shape to be that of a static meniscus. When the meniscus is pinned to a solid surface, as in pre-metered slot coating (Ruschak, 1976),the contact line may or may not be the only stagnation line on the upstream part of the meniscus. If it is, Esmail and Hummel’s approach should be useful provided the meniscus is not highly curved and a reasonable approximation to the underlying flow is available. Perturbation Analysis If a film flow has a limit in which the governing equations (8)-(13) simplify enough to be solvable, departures from the limiting base case can be studied by perturbation analysis, for which mathematical techniques are well developed (Van Dyke, 1975; Cole, 1968). Suitable small parameters have to be found in which to expand the perturbation solution. For examples, in surface leveling the base case may be the rest state and the parameters the ratios ( h - ho)/ho and Xlho; in dip coating the base case may be a state of vanishing film thickness and the small parameter is the capillary number ~ U o / u ;
in rimming flow the base case is solid body rotation and the parameters the Froude number g / w 2 R and a Weber number u/po2R3.In free surface flows the unknown meniscus location as well as velocity and pressure have to be expanded in appropriate powers of the small parameters. When the strategy is most successful the expansions substituted in the governing equations generate a set of linear problems for successively smaller adjustments to the base case, and these problems can be solved sequentially. For rimming flow at high Weber number (comparatively large surface tension) and small Froude numbers, Le., centripetal acceleration much larger than gravity, Ruschak and Scriven (1976) showed that perturbation in Froude number about solid body rotation is mathematically regular, and they were able to solve the first-order problem in closed form, a result useful in understanding the fluid mechanics of the flow and in validating a numerical solution strategy developed for more general cases (Orr, 1976; Orr and Scriven, 1977). The perturbation strategy could be carried to higher orders, except that direct numerical solution is more fruitful. For dip coating at small capillary numbers, i.e., viscous stresses much smaller than capillary pressure (cf. eq 6 and 121, Ruschak (1974) showed that perturbation in capillary number is mathematically singular. The zero-order solution, or base case, proves to be a liquid film of zero thickness which upstream joins a static meniscus along a contact line at which the contact angle is zero (measured through the liquid). The first-order solution turns out to be in essence the solution found by Landau and Levich (1942);unlike their strategy the perturbation analysis reveals precisely where the dynamic meniscus must be matched upstream with a static meniscus that extends outward over the coating bath (Ruschak, 1974). The higher-order problems are also linear but require solving quite complicated partial differential equations for further adjustments to the velocity and pressure fields. Again, direct numerical solution of the full problem is much more likely to be fruitful. Another type of analysis which applies to unsteady perturbations is important for studying the stability of a viscocapillary flow found by any strategy. That is the method of small disturbances, or infinitesimal disturbances if the resulting equations are to be made linear (Joseph, 1976). Examples are not yet available, except for the much-studied stability of uniform, falling liquid films on stationary substrates. A particularly interesting technique was initiated by Benney (1966), developed by Krantz and Goren (1970) and applied by Atherton and Homsy (1976) for the evolution of finite-amplitude disturbances of long wavelength. Pressure and velocity fields are expanded in terms of the appropriate small parameter and the expansions are substituted in (8)(lo), (12)-( 13)-excluding the kinematic condition-to produce a set of differential equations which can be solved sequentially in terms of the unknown meniscus profile. One or more of these solutions can then be used to estimate the velocity field a t the meniscus for substitution in the kinematic condition, which is a partial differential evolution equation for meniscus shape. The resulting equation in practice has to be solved numerically.
Numerical Analysis Rarely does the meniscus of a film flow coincide with any standard coordinate surface. Consequently, a convenient set of basis functions for the flow domain is seldom available. Moreover, the domain itself is indistinct if the upstream and downstream conditions are known only asymptotically, as is often true. Often, too, different forces and accelerations are important in different regions, and gradients of meniscus location, pressure, and velocity components vary widely over the flow domain. When the need is for the meniscus shape and
the convective action close beneath the interface, a less accurate approximation throughout the rest of the film may be acceptable, as in the neighborhood and integral balance strategies described above. Among weighted residual methods for calculating approximate solutions to boundary value problems (Finlayson and Scriven, 1966; Finlayson, 1972), the one best suited to these circumstances is the subdomain method with a separate basis for each subdomain, and either the Galerkin or collocation procedure for determining the coefficients of the basis functions. Moreover, the complications in such an approach are reduced by dividing the flow domain into enough microdomains, usually triangular or quadrilateral in two-dimensional problems, that a single, simple polynomial basis function suffices on each microdomain, or element. This is the finite element method (Strang and Fix, 1973; Oden, 1972). In it the solution of the problem is expanded in terms of basis or trial functions each of which vanishes everywhere but one group of microdomains (the basis functions are said to have compact support). Thus the finite element method generates an analytic expansion approximation akin to Fourier-Bessel and more general eigenfunction expansions on separable macrodomains. Free boundary and nonlinear problems require iterative procedures, however. The finite element method is commonly regarded as a means of generating a numerical approximation and is compared with the older finite difference approach. Standard finite difference schemes are unsuitable for most film flow problems, for the same reasons that make the finite element method attractive. The disadvantages can be overcome by regional mesh refinement and marker-and-cell devices (Harlow and Welch, 1965; Amsden and Harlow, 1970a,b) although MAC techniques have been applied to time-dependent capillary free-surface flows with success (Nichols and Hirt, 1971;Daly, 1969). Generally, however, the finite element method appears to be superior, partly by virtue of its much closer relation to the integral conservation equations from which the differential equations of film flow are derived. This opinion is supported by calculations of static meniscus shapes (Orr et al., 1975) and a recent numerical simulation of rimming flow by means of the finite element method (Orr, 1976; Orr and Scriven, 1977). Meniscus location and the velocity and pressure fields are calculated for various combinations of Froude number, Reynolds number, Weber number, and liquid loading of the rotating cylinder. In many cases inertial, gravitational, pressure, viscous, and capillary forces are all important, and in some cases a recirculation zone exists which could not be predicted by the Protean coordinate strategy. Thus the finite element simulation extends the work of Deiber and Cerro (1976). It also extends the perturbation analysis by Ruschak and Scriven (1976) and the extension makes clear the limits of applicability of that analysis. The finite element method can be applied to dip-coating and other steady film flows. It can be extended to handle time-dependent film flows such as surface leveling. Numerical simulation, though often more expensive than a group of approximations based on simplified equations, does permit exploration of wide ranges of physical variables. It should be possible to probe in a parameter space to find conditions for the onset of numerical instabilities which may correspond to hydrodynamic instabilities of importance in the real system.
Conclusions There is an important question of the optimal strategy for a given viscocapillary flow. The six basic scalar equations of two-dimensional film flow, (8)-(13), can be transformed and combined into mathematicany equivalent sets in a large number of ways, some of which are indicated above. It is plain Ind. Eng. Chem., Fundam., Vol. 16, No. 4 , 1977
399
that different strategies lend themselves to different types of simplifications and approximation methods, the results of which have differing ranges of validity. At issue is the allocation of error over the equations of the set and, so far as the momentum and continuity equations for the liquid are concerned, over the neighborhood of the meniscus as against the deeper part of the film, and over various regions along a flowing film. Likewise in numerical simulation there are choices to be made of the distributions of element sizes, basis or trial functions, and convergence criteria. In open systems there is likely to be a question about where upstream and downstream to switch from error management to demanding that preconceived boundary conditions be met. When the meniscus is elevated to the status of system and the underlying liquid film, though of higher dimension, is demoted to the role of boundary, a choice has to be made among the meniscus equations, (11)-(13) and all their possible combinations, for the one or two to serve as the meniscus shape and evolution equations. Of particular interest is the newly derived shape equation (23). The parallel choice in numerical simulation is which equation to use as the criterion for relocating the meniscus in successive iterations (Orr and Scriven, 1977). The normal stress equation (12) may be preferable because it alone depends on the highest derivatives (second) of meniscus location. For a time-dependent flow, the kinematic condition (11)has also been used (Benney, 1966). Apparently the shear stress equation (13) has not been tried. The theory of meniscus shape in film flows which is synthesized here is still an open subject. It is becoming possible to proceed more systematically than intuitively in attacking a specific problem, in marshalling the available strategies and devising new ones. The theory will be advanced by accumulating experience with overlapping combinations of the more promising strategies. Acknowledgment This research stemmed from the Minnesota coating flows seminar. We are grateful to W. T. Donahue for helpful comments and suggestions. Partial support was provided by the National Science Foundation. B. G. Higgins and W. J. Silliman received partial support from the Council for Scientific and Industrial Research, South Africa, and the Bush Foundation, respectively. Nomenclature 3 = normal stress boundary condition, N/m2 g = acceleration of gravity, m/s2 g = gravity vector, m/s2 Ex = acceleration of gravity in x direction, dimensionless = acceleration of gravity in y direction, dimensionless h = height of meniscus, m h = height of meniscus, dimensionless ho = datum thickness, m h, = (ah/&),dimensionless H = h/ho, dimensionless H = mean curvature, m-1 L, = x-length scale, m L, = y-length scale, m n = normal to surface, dimensionless N = dimensionless group (see Table I) p = pressure,N/m2 = pressure, dimensionless P = pressure unit, N/m2 q = volumetric flow rate, m3/s Q = volumetric flow rate, dimensionless 4 = time, s t = time, dimensionless To = h(pg/MUO)1/2, dimensiofiless T = stress tensor, N/m2 400
Ind. Eng. Chem., Fundin., Vol. 16, No. 4, 1977
~i, =
components of stress tensor, N/m2
u = x velocity, m/s u = x velocity, dimensionless U = u/Uo, dimensionless UO = velocity of web, m/s U, = n velocity unit, m/s = y velocity, m/s = y velocity, dimensionless v = velocity vector, m/s v, = velocity of solid surface, m/s V = u / ~ ~ / ~ N c ~ ~ dimensionless /~HUO, V, = y velocity unit, m/s x = distance,m F = distance, dimensionless X = 31/3N~a1/3x/ho, dimensionless y = distance,m 7 = distance, dimensionless Y = y/h, dimensionless
-uu
Greek Letters a = L / L x ,dimensionless ti = (I?+ h, 2)1/2, dimensionless /3 = V,/U,, dimensionless V, = surface gradient operator, m-l t = p&,, dimensionless 70 = deviation of thickness from datum, m 0 = time unit, s X = wavelength,m p = viscosity, kg/m-s v = kinematic viscosity, m2/s p = density, kg/m3 u = surface tension, N/m w = angular velocity, s-1
Literature Cited Amsden, A. A., Harlow, F. H., J. Comput. Phys., 6, 322-325 (1970a). Amsden, A. A., Harlow, F. H., Los Alamos Scientific Laboratory Report LA-4370, (1970b). Anshus. B. E.. Am. Chem. Soc., Div. Ora. Coat. Plast. Chem. PaD., 33, No. 2, 493-502 (1973). Atherton, R. W.. Homsy, G. W., Chem. Eng. Commun., 2,57-77 (1976). Basset, A. B., "A Treatise on Hydrodynamics," George Bell and Sons, London, 1888. Benney, D. J., J. Math. Phys., 45, 150-155 (1966). Cole. J. D.,"Perturbation Methods in Applied Mathematics." Blaisdell, London, 1968. Daly, B. J.. J. Comput. Phys., 4, 97-117 (1969). Degani, D., Gutfinger, C., lsr. J. Tech., 12, 191-197 (1974). Deiber, J. A., Cerro, R. L., lnd. Eng. Chem., Fundam., 10, 102-110 (1976). Duda, J. L., Vrentas, J. S., Chem. Eng. Sci., 22, 855-869 (1967). Emmert, R. E., Pigford, R. L., Chem. Eng. Prog., 50, 87-93 (1954). Esmail, M. N., Hummel, R. L., Chem. Eng. Sci., 30, 1197-1198 (1975a). Esmail, M. N., Hummel, R. L., AlChEJ., 21, 958-965 (1975b). Finlayson, B. A., "The Method of Weighted Residuals and Variational Principles," Academic Press, New York, N.Y., 1972. Finlayson, B. A., Scriven, L. E., Appl. Mech. Rev., 19, 735-748 (1966). Fulford, G. D., Adv. Chem. Eng., 5, 151-229 (1964). Gardner, G. C.. Adebiyi, G. A., Chem. Eng. Sci., 29, 461-473 (1974). Groenveld, P., Ph.D. Thesis, Delft, 1970. Gutfinger. C.. Tallmadge, J. A,, AlChEJ., 11, 403-413(1965). Harlow, F. H., Welch, J. E., Phys. Fluids, 8, 2182-2189 (1965). Higgins, B. G., Scriven, L. E., submitted to lnd. Eng. Chem., Fundam., 1977. Johnstone, H. F., Pigford, R. L., Chapin, J. H., Trans. AlChE, 37, 95-133 (1941). Joseph, D. D.. "Stability of Fluid Motion," Springer, New York, N.Y., 1976. Krantz. W. B., Goren, S. L., lnd. Eng. Chem., Fundam., 9, 107-113 (1970). Lamb, S. H., "Hydrodynamics," Dover Publications, New York, N.Y., 1945. Landau. L., Levich, E., Acta Phys. U.R.S.S., 17, 42-54 (1942). Lee, C. Y., Tallmadge, J. A., lnd. Eng. Chem., Fundam., 14, 120-126 (1975). Lee, C. Y., Tallmadge, J. A,, Ind. Eng. Chem., Fundam., 15, 258-266 (1976). Nichols, B. D., Hirt, C. W., J. Comput. Phys., 8, 434-448 (1971). Oden, J. T., "Finite Elements of Nonlinear Continua," McGraw-Hill, New York, N.Y., 1972. Orchard, S. E., Appl. Sci. Rev., 11, 451-464 (1962). Orr, F. M., Jr., Ph.D. Thesis, University of Minnesota, Minneapolis, 1976. Orr, F. M., Jr., Scriven, L. E., J. Fluid Mech., in press, 1977. On,F. M.,Jr., Scriven, L. E., Rivas, A. P., J. Colloidlnterface Sci., 52, 602-610 (1975). Quach, A., lnd. Eng. Chem., Prcd. Res. Dev., 12, 110-116 (1973). Rao. M. A., Throne, J. L., Polym. Eng. Sci., 12, 237-264(1974). Ruschak, K. J., Ph.D. Thesis, University of Minnesota, Minneapolis, 1974. Ruschak, K. J., Chem. Eng. Sci., 31, 1057-1060 (1976). Ruschak, K. J.. Scriven, L. E., J. FluidMech., 76, 113-125(1976). Shewood, T. K., Pigford, R. L., Wilke, C. R., "MassTransfer," McGraw-Hill, New York. N.Y., 1975.
Spiers, R. P., Subbaraman, C. V., Wilkinson, W. L., Chem. Eng. Sci., 29,389-396 (1974). Spiers, R. P., Subbaraman, C. V., Wilkinson, W. L., Chem. Eng. Sci., 30,379-395 (1975). Strang, G.. Fix, G. J., "An Analysis of the Finite Element Method," Prentice-Hall, Enalewood Cliffs, N.J.. 1973. Van dyke, M., "Perturbation Methods in Fluid Mechanics," Parabolic Press, Stanford, Calif., 1975. Van Rossum, J. J., Appl. Sci. Res. A, 7, 121-144 (1958).
Wehausen, J. V., Laitone, E. V., "Handbuch der Physik," Vol. 9, pp 466-779 Springer-Verlag. New York, N.Y., 1960. White, D. A., Tallmadge, J. A., Chem. Eng. Sci., 20, 33-37 (1965). Williamson, A. S.,J. FluidMech., 52, 639-659 (1972).
Received for reuiew February 16,1977 Accepted June 24,1977
Some Theoretical and Experimental Observations of the Wave Structure of Falling Liquid Films F. Wayne Plerson and Stephen Whitaker' Department of Chemical Engineering, University of California, Davis, California 956 16
The linear stability problem associated with a vertical liquid film flowing under the action of gravity has been solved in terms of a numerical solution of the Orr-Sommerfeld equation. Results have been obtained for both the temporal and spatial representations of growing waves and for a wide range of the physical parameters characterizing this type of flow. Experimental values of the wavelength and wave velocity have been determined for water films and the results are in reasonably good agreement with the theory. Calculated values for the wavelength, wave velocity, and growth rate of the most unstable wave indicate only a small difference between the temporal and spatial formulations for water films.
1. Introduction
The wave characteristics of thin, falling liquid films have long been of interest to engineers because of the prevalence of thin liquid films in industrial processes. In addition, there is a general interest in the phenomenon as it provides an attractive challenge for the comparison of stability theory with experiment. Our own interest in this problem stems from the sometimes dramatic effect of surface active agents on the wave characteristics of thin liquid films. Under certain conditions small amounts of surfactant can completely inhibit wave formation to the extent that the liquid film appears as motionless as a plate of glass. Davies and Rideal (1963, p 266) have commented on this phenomenon, and have provided some photographic observations. As a prelude to an investigation of some of the anomolous experimental results for surfactant solutions (Strobe1 and Whitaker, 1969; Whitaker, 19711, a detailed survey of linear stability theory for pure liquids is in order. Linear stability theory has severe restrictions regarding the range of application for falling liquid films. Such films may be formed by fluid issuing from a channel or flowing over a weir as illustrated in Figure 1.Both types of flow have been analyzed theoretically by Cerro and Whitaker (1971a, 1974) and the calculated values for the free surface velocity and the film depth are in good agreement with the experimental data of Lynn (1960) and of Cerro and Whitaker (1971b). When the original channel depth in Figure l a is equal to the final film depth, h a , the distance required for the surface velocity to reach 99% of the final value is 0 . 5 h - N ~ This ~ . is an order of magnitude larger than the entrance lengths for pipes or rectangular channels; however, the transition here is from a confined flow to a free surface flow and we do not expect the entrance lengths to be the same. For falling liquid films the entrance length is somewhat sensitive to the way in which the film is formed, but a reasonable estimate of the length is given by
Le
N
hmNRe
(1.1)
where N R is~the Reynolds number based on the free surface velocity of the uniform flow and h mis the film thickness for the given uniform flow. These quantities are defined by ha = (3q~/g)l/~ u, = hm2g/2v
(1.2)
(1.3)
where q is the volumetric flow rate per unit width, g is the gravitational constant, and v is the kinematic viscosity. The film thickness can be written in terms of the Reynolds number h, = ( ~ N R~,' / g ) ' / ~ so that the entrance length can be expressed as
Le
-
( 2v2/g)1/3N~e4/3
-
(1.4) (1.5)
For water we have ( 2 ~ ~ / g ) l 6/ ~X cm, and an estimate of the entrance length as a function of the Reynolds number is given by Le
-
[6 X 10-3N~e4/3]cm
(for water)
(1.6)
This provides us with an indication of the region in which the flow is nonuniform and the usual small disturbance analysis leading to the Orr-Sommerfeld equation may not apply. For small Reynolds numbers, say N R ~ 1,this region is negligible; however, for N R =~ 100 the entrance length is estimated to be about 3 cm and a t a Reynolds number of 500 we find Le to be on the order of 24 cm. Certainly the entrance region cannot be ignored in any comprehensive study of the stability characteristics of falling liquid films. The distance down the film at which waves are first visible is referred to as the wave inception line, and the location of this line has been tabulated by several investigators (Pierson, 1974; Cerro and Whitaker, 1971b; Strobe1 and Whitaker, 1969; N
Ind. Eng. Chem., Fundam., Vol. 16, No. 4, 1977
401
