- Email: [email protected]
Analysis of surface tension driven leveling in viscoelastic films
Journal of Non-Newtonian Fluid Mechanics, 22 (1987) 219-233 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
219
ANALYSIS OF SURFACE TENSION DRIVEN LEVELING IN VISCOELASTIC FILMS
R. KEUNINGS
and D.W. BOUSFIELD
Center for Advanced Materials, Lawrence Berkeley Laboratory, Berkeley, CA 94720 (USA)
University of California,
(Received June 3, 1986)
We study the surface tension driven leveling of surface irregularities in viscoelastic films deposited on a flat surface. Analytical results are presented for the generalized Maxwell model in the limit of infinitesimal surface disturbances. The results show that elasticity retards the leveling process. In some cases, the retardation effect is dramatic. Keunings’ finite element algorithm for solving viscoelastic moving boundary problems is used to analyze the case of finite-amplitude surface disturbances. Results for an Oldroyd B fluid display significant elastic effects which are consistent with the asymptotic predictions.
1. Introduction Coating operations involve small-scale free surface flows by which a solid substrate is covered with a thin liquid film. The appearance and performance of the finish greatly depend on the ability of the liquid film to level the surface irregularities generated during the coating process. In paint coating, for example, it is often observed that the film left behind the brush exhibits nearly sinusoidal disturbances (the so-called brush-marks). These disturbances are generated by the ribbing instability mechanism that takes place in the wake of coating devices (see e.g. [l-4]). Another source of surface irregularities is found in the powder coating process, whereby a multi-layer of polymer particles is deposited on a metallic substrate; upon heating, the polymer melts and covers the substrate by spreading and sintering. Surface corrugations result from local .non-uniformities of the powder deposition. The extent to which these corrugations will level during 0377-0257/87/$03.50
0 1987 Elsevier Science Publishers B.V.
220 the remainder of the bake cycle will affect the quality of the finish (see e.g. F-71). The physics and mechanics of leveling have been studied since the 1920’s (see the review by Quach [S]). The driving force is generally a combination of surface tension and gravity, although the latter can be safely neglected in many applications. Most of the work has focussed on the surface leveling of viscous, horizontal films of finite thickness and infinite extent. It is usually assumed that the amplitude of the corrugations is small compared to the mean film thickness, and that inertial forces can be neglected in comparison with viscous forces. In this context, Orchard [9] has proposed a linear theory of the surface tension driven leveling of two-dimensional disturbances in Newtonian films. The author expressed an arbitrary (but small) surface disturbance as a Fourier series, and obtained an expression relating the decay rate of each individual Fourier mode to the wave number, the equilibrium thickness of the film, the surface tension coefficient, and the fluid viscosity. Similar studies that include the effects of pseudoplasticity, thixotropy, and loss of volatiles are reviewed by Quach [8]. In view of its asymptotic character, Orchard’s analysis applies at the final stages of the leveling process. The case of finite-amplitude disturbances is difficult to analyse since it involves the solution of a complex transient free surface flow. Numerical procedures are clearly needed. Degani and Gutfinger [lo] used a version of the marker and cell method to simulate their laboratory experiments with glucose-water solutions. Good agreement with experimental data was reported. More recently, Kheshgi and Striven [ll] used a penalty finite element method for unsteady free surface flows to study the capillarity driven, viscosity dominated leveling of finite-amplitude sinusoidal corrugations (see also Kheshgi [12]). They found that the deeper region of the film levels faster than the shallower region at the early stages of the leveling process. This non-linear effect is observed in practice [13]. At later times, the leveling rate is accurately predicted by Orchard’s analysis. All these investigations have dealt with inelastic fluids. Viscoelasticity has long been recognized as an influential factor, however [8]. To our knowledge, the only attempt to determine viscoelastic effects on the leveling process dates back to 1968 [14]. The author extended Orchard’s linear analysis using the Rivlin-Ericksen second-order model. Since then, it has been established that such a model should not be used in the analysis of unsteady flows (see e.g. [15]). Re-examination of the problem is thus called for. In the present paper, we study the surface tension driven leveling of surface irregularities in viscoelastic films. First, we generalize Orchard’s asymptotic analysis in a way that is consistent with the linear viscoelastic limit of Simple Fluid Theory. To this end, we use the general linear
221 viscoelastic model with a relaxation function defined as a discrete sum of exponent&. The results show that elasticity retards the leveling process. The effect can be drastic. Numerical results for finite-amplitude corrugations are obtained with Keunings’ finite element algorithm for viscoelastic moving boundary problems [16]. We use the Oldroyd B constitutive equation, a useful model in the context of dilute polymeric solutions [17]. The results exhibit the drastic retardation effect of elasticity found with the linear theory. At the early stages of the leveling process, we predict non-linear effects similar to those observed in the Newtonian case by Kheshgi and Striven [ll]. At later times, the asymptotic theory is seen to provide an accurate measure of the leveling rate. 2. Formulation A schematic of the leveling problem is shown in Fig. 1. We assume that the flow is incompressible, isothermal, and two-dimensional. Inertial and body forces are neglected. The set of conservation laws reduces thus to -vp+v.T=O,
0)
vau=o,
(2) where p is the pressure, T is the extra-stress tensor, and 2) is the velocity vector. A constitutive equation is needed to relate the extra-stress tensor to the deformation experienced by the fluid. Our selection will be discussed later in the text. The unknown fields T, u, and p depend on time t and the space variables x and y. In the present application, the integration domain is also an unknown function of time. Its evolution is determined through the kinematic condition g
+ “$
= vu,
(3)
which relates the film thickness h to the velocity components free surface (Fig. 1).
ho
v,, vv at the
h
I //////////////////////////////// Substrate Fig. 1. Schematic of the leveling process. The liquid film extends from x = - COto x = + 00.
222 We consider the evolution of disturbances applied to the equilibrium film thickness. The surface disturbance is denoted by 6, i.e. h(x,
t) = h, + S(X, t),
(4
where h, is the thickness of the unperturbed boundary conditions apply: (i) stress balance at the free surface (y = 6):
(-pz+
T).n=a
a*h/ax*
[l + (i3h,ax)2]3'2n'
(flat) film. The following
(5)
where n is the normal to the free surface and u is the coefficient of surface tension [18]. The atmospheric pressure has been arbitrarily set to zero in (5). (ii) no-slip at the substrate (y = -h,): u, = uv = 0.
(6)
The selection of initial conditions for the film thickness and the velocity and stress fields is discussed hereafter. 3. Linear analysis We first consider a linearized version of the leveling problem valid for infinitesimal disturbances (S small). In this context, the choice of a specific rheological model is largely facilitated. We use in the present paper the general linear viscoelastic model with a relaxation function defined by a sum of N exponential terms [19]. The resulting constitutive equation is known as the generalized Maxwell model. It reads T=
&
(7)
i=l
8T. q. + r& = 2piD. The symbol D denotes the rate of strain tensor i( vo + VU=). The constants pi and 7i are viscosity coefficients and relaxation times, respectively. When all the oh’svanish, the model (7)-(8) defines a Newtonian fluid with viscosity p = C~v-l~i. Hereafter, we shall refer to this particular case as the Newtonian counterpart of the viscoelastic model (7)-(8). After linearization, the kinematic condition (3) reduces to
as/at = vy.
(9)
223 The free surface condition (5) becomes TX,,= 0, -P + T, = u
a%/ax2,
and is specified at the known location y = 0. 3.1 Analytical
solution
Following a standard Fourier procedure, we seek for each unknown TX,, TYY, TX,, ux, uv, and p, a solution of the linearized equations in the complex form
f (x, y, t) =f( y) eikxpat,
(12)
where k = 2n/A is the disturbance wave number and (Y is the leveling rate. The surface disturbance 6 is sought in the-form 6(x,
t) = S eikxpat,
(13)
where the coefficient 8 is a constant. In the context of this linearized analysis, the solution corresponding to an arbitrary (but small) thickness disturbance is obtained by linear superposition of Fourier modes (12)-(13). Inserting (12) in the momentum equation (1) and eliminating the pressure, we obtain ikcX + c;
+ k2E,
- ikqY = 0,
(14)
where the symbol ’ denotes the derivative with respect to y. On the other hand, the constitutive equations (7)-(8) yield TX = 2p.,[ik&],
05)
T,=21r&],
06)
<,,=py[C”+ikUy],
(17)
with cl,, defined as
Note that equations (15)-(17) have the form of a Newtonian constitutive model relating the Fourier coefficients of T and D, with the viscosity p replaced by IL,,_ The incompressibility constraint (2) gives
224
and with the help of (15)-(17) and (19), eqn. (14) reduces to an ordinary differential equation for C__:
- -kZ 12u,=o. d2
(20)
[ dY2 The general solution of this equation is Ei,(y)=a
cash ky+b
sinh ky+cycosh
ky+dy
sinh ky,
(21)
where a, b, c, d are constants determined by the boundary conditions. The linearized interface conditions (lo)-(11) yield ak=
b=--8,
-d,
ok
(22)
2l-G
and the no-slip condition (6) at the substrate gives c= -kb+(t). a = W(5), Here, 5 = 2?rh &I, and the functions 4, (p are defined by
(23)
4w = tanh 5 - 5 sech2t 1
HO=
+ E2 sech2t
’
1 + t2 1sech25 -
(24) (25)
Finally, use of the linearized kinematic condition (9) provides an implicit expression for the leveling rate CC
where I_L ,,( CX)is given by (18). Equation (26) is formally equivalent to the result of Orchard [9] for a Newtonian fluid, with p replaced by pV( a). 3.2 Newtonian fluid
The case of a Newtonian flu? has been analysed in detail by Orchard [9] and will serve as a reference SC ion. In order to investigate the effect of elasticity on the leveling rate, we freeze the values of the parameters (A, h,, u, p), and define the Newtonian leveling rate according to (26) (27)
For thin films (5 -C l), (27) becomes
(28)
225 TABLE 1 Leveling rate and half time as functions of the relaxation time; one-mode Maxwell model r(s) ct
(s-l)
h/Z
(4
0.
0.1
1.
10.
50.
0.088 8.
0.087 8.
0.081 9.
0.047 15.
0.016 43.
100. 0.009 77.
This result will be used later to define a time scale for the non-linear computations. 3.3 One-mode Maxwell fluid For a single Maxwell element, i.e.p = pr, 7 = 7r, and N = 1 in (7)-(8), obtain from (US), (26) and (27)
a= 1+
aNewt raNewt
-
we
(29)
Since the Newtonian leveling rate is always positive [9], the present analysis predicts that elasticity has a retardation effect on the leveling process. The following example shows that the leveling rate can be significantly reduced with highly elastic fluids. Consider a film of thickness h, = 0.01 cm and a disturbance of wave length X = 0.05 cm. Realistic numbers for a polymer melt are lo4 poise for the viscosity, and 40 dyne/cm for the surface tension coefficient. We show in Table 1 the leveling rates predicted on the basis of (27) and (29) for values of the relaxation time between 0 and 100 s. Over this range, the leveling rate experiences a ten-fold decrease from 0.088 s-r to 0.009 s-l. A useful time constant for the leveling process is the half time t,,, defined as the time required for the amplitude of the disturbance to decrease by a factor of two. It follows from (13) that $2 = a - ’ In 2. Values of this time constant are also shown in Table 1. 3.4 Jeffreys’ model If we particularize the constitutive model (7)-(8) 72 = 0, we obtain Jeffreys’ model:
with N = 2, 7 = 7i, and
where 7 * is the retardation time p2~/p and p= ~1~+ p2. Equation (30) is the linear viscoelastic limit of the dumbbell model for dilute polymeric
226 solutions [20]; ~1~and p are the solvent and solution viscosities, respectively. For this model, we find (y=
1 - rcy ~_T*aaNewt*
(31)
Equation (31) is quadratic in 0~.It can be shown that its roots are always real and positive. Clearly, the root of least magnitude will dominate the leveling process. Inspection of (31) reveals that at least one of the roots is smaller than aNewt, which implies that elasticity retards the leveling process. This result is obvious if one root satisfies 1 - TCX > 0. The case 1 - ~a -K0 and 1 - alar > 0 is impossible since (Yis positive. Now, assume that one root & is such that both 1 - T(Y(‘) and 1 - T*&) are negative. This implies that (Y(l)2 aNewt.Then, the second root (Y(*)must be smaller than aNewtsince the product of the two roots is equal to aNewJr*. For the sake of illustration, consider the polymeric solution (Boger fluid Bll) studied by Jackson et al. [21]. The material constants are in this case p = 1040 poise, T = 2.54 s, and T* = 1.97 s. We adopt the values of the previous example for the other parameters of the leveling problem. The Newtonian leveling rate is equal to 0.85 s-l, while the smallest root of (31) is 0.34 s-l. The corresponding values of the half time t,,, are 2 s for the viscoelastic fluid, and 0.8 s for its Newtonian counterpart. The retardation effect of elasticity is again significant. 3.5 N-mode Maxwell model For a discrete spectrum of N relaxation times, the implicit definition of the leveling rate (26) is a polynomial equation of degree N. As a result, the analysis becomes difficult when N >, 3 and direct computation of the roots of (26) for specific spectra { ri, pi} is called for. To do so, we have used the symbolic manipulation system MACSYMA [22] which implements an algorithm of Jenkins [23] for finding the roots of real polynomials. In ah our computations, the solutions of (26) were found to be real positive and the root of least magnitude was smaIler than the corresponding Newtonian leveling rate. We describe below two examples of such computations for values of the leveling parameters (X, h,, a) identical to those used in our previous illustrations. The results obtained with the discrete spectrum are compared to those obtained with either the largest relaxation time pi or the average relaxation time T,, defined by
(32) Our first example makes use of the relaxation data tabulated by Laun [24] for a low-density polyethylene similar to the IUPAC sample A (Meissner
227 TABLE 2 Relaxation data for LDPE Melt 1 at 150°C (Lauu [24]); c = 511000 Poise, 7av= 59 s
i
,a
pi (Poise)
1 2 3 4 5 6 7 8
1000. 100. 10 1. 0.1 0.01 0.001 0.0001
10000. 180000. 189000. 98000. 26700. 5860. 948. 129.
[25]). The 8-mode spectrum is given in Table 2. The Newtonian leveling rate is 0.0017 s-l, while the smallest root of (26) is 0.00097 s-l in the viscoelastic case. These values yield a half time t,,, of 404 s and 712 s, respectively. The one-mode analysis using the average relaxation time gives an overestimated value (0.0015 s-l) for the viscoelastic leveling rate. The largest relaxation time, on the other hand, yields a value of 0.00063 s-l which is closer to that obtained with the discrete spectrum. A more drastic retardation effect of elasticity is obtained with the data of Greener and Connelly [26] for an aqueous solution of polyacrylamide (Table 3). Here, the Newtonian leveling rate is 5 s-l and the smallest root of (26) is 0.01 s-l. The corresponding values of t,,, are 69 s for the polymeric solution and 0.14 s for its Newtonian counterpart. In this case, use of the average relaxation time yields a good approximation (0.06 s-‘) of the
TABLE 3 Relaxation data for au aqueous solution of polyacrylamide (Greener and Couuelly [26]); p = 168 Poise, T,, = 16 s
i
Ti
1 2 3 4 5 6 7 8 9 10 11
100. 32. 10. 3.2 1. 0.32 0.1 0.032 0.01 0.001 O.O@Ql
ls)
pi (Poise) 7.6 44. 41. 34.4 20.6 11.8 4.8 2.42 0.79 0.845 0.0254
228 leveling rate. The result obtained on the basis of the maximum relaxation time is even better (0.01 s-l). The leveling process appears to be dominated by the largest relaxation time, at least in the two examples discussed here. 4. Non-linear fiite
element analysis
In view of its asymptotic character, the above analysis applies at the final stages’of the leveling process where the fluid is close to the state of rest. It is quite general in the sense that the viscoelastic model (7)-(g) constitutes a good approximation of any Simple Fluid in the linear viscoelastic limit. The issue of finite amplitude disturbances of the film thickness is of course more difficult to address in a general way. Indeed, the results will be specific to a constitutive model (which must be objective), and a set of initialconditions. Moreover, analytical developments are impossible since the flow problem becomes highly non-linear. These difficulties are typical of the analysis of complex viscoelastic flows. We describe in this section the results of a non-linear finite element analysis of the leveling process. We have selected the Oldroyd B model to define the rheology of the film: T+G=2p
D+r*vD [
. I
(33)
Here, the symbol v denotes the upper-convected derivative. For example,
V T=g+oe~~-~.~~--vuT.~.
(34)
The material constants EL,7, and 7 * have the same meaning as in (30). The Oldroyd B model is a generalization of Jeffreys’ equation for arbitrary deformations and can be derived from the dumbbell kinetic theory [20]. It has proved quite useful in the modeling of dilute polymeric solutions (see e.g. [27,28]). We consider in the present section the evolution of a finite amplitude disturbance applied to the film kness. The initial disturbance is sinusoidal with a wave length A. We assume that the free surface remains spatially periodic as it levels (it will of course depart from a sinusoidal shape because of non-linear effects.) In addition, we restrict our attention to disturbances that are symmetric about the plane x = X/2 (Fig. 2). This assumption translates into the boundary conditions u
x
_au,=ah=, ax
ax
’
(35)
specified at the crest (x = 0) and trough (x = X/2) of the free surface. With
229
X 0.
x/2
Fig. 2. Typical geometryand finite element mesh for the finite-amplitudecomputations.
this simplication, the integration domain extends from x = 0 to x = X/2. The mathematical problem now consists of solving the governing equations (l)-(3) and (33) in terms of T, u, p, and h, with the boundary conditions (S)-(6) and (35) and proper initial conditions. Thus formulated, the leveling process is a two dimensional, transient, free surface viscoelastic flow driven by surface tension. The numerical results described hereafter have been obtained by means of the finite element algorithm developed by Keunings [16] for the solution of viscoelastic moving boundary problems. Briefly, the method consists of a Galerkin principle invoked on deforming finite elements. The motion of the nodes is anchored to the displacement of the free surface in a way that prevents gross distortion of the finite element layout. The Gale&in formulation takes this nodal motion into account and contains the standard Eulerian and Lagrangian approaches as particular cases. A semi-implicit predictorcorrector scheme is used for the temporal integration. This technique has been applied successfully in the analysis of ‘the surface tension driven breakup of viscoelastic jets [16,28]. The flow domain is discretized by means of nine-node Lagrangian elements. A typical finite element mesh is shown in Fig. 2. The free surface representation h is approximated by Co - P* one-dimensional elements; we use Co - P’ interpolations for the pressure and the extra-stress tensor, and
230 Co-P2 basis functions for the velocity. The initial free surface shape is set to h(x, 0) =h,[l
+r cos(2Izx/A)],
(36)
where c is the dimensionless amplitude of the disturbance. The other three dimensionless groups that arise in the present application are the Deborah number De, the dimensionless wave number 5, and the ratio /3of characteristic times:
De=7
16m4ahfj 3jLx4
’
t=,,
2nh,
(37)
In our definition of the Deborah number, we have used the result (28) to obtain a typical time scale of the flow. 4. I Newtonian fluid
The evolution of finite-amplitude disturbance in Newtonian films has been studied by Kheshgi and Striven [ll] (see also [12]). The authors used a penalty finite element method for solving the moving boundary problem. The results for creeping flow indicated that the thick region of a sinusoidally disturbed film levels faster than the thin region. This non-linear effect
1
.c
1:
0.9
%
%
2 y 1.2
0.8 E 2 l-
z
I-
E
.EI1.3
0.7
4
1.4
0.6
1.5
0.5
E i ‘E .I
Time
Fig. 3. Maximum and minimum film thicknesses as a function of time (Newtonian fluid); finite-amplitude computations: -, linear theory: - - - - - -, Kheshgi’s results [12]: l ; time ant thicknesses are made dimensionless with 3~A4/16a4uh i and h,,, respectively.
231 deforms the initially sinusoidal shape of the free surface. In the particular case of creeping Newtonian flow, Keunings’ algorithm contains the classical displacement (u-u-p) and mixed Gale&in formulations of the Stokes problem [29]. We have repeated with both methods a simulation of Kheshgi [12] for 5 = r/10 and c = 0.5. Figure 3 shows the evolution of the maximum and minimum film thicknesses as they approach the steady state h = h,. The results obtained with the displacement method are indistinguishable from those of the mixed technique at the scale of the drawing. The agreement with Kheshgi’s computations is excellent. Also shown for comparison is the linear prediction derived from (13) and (27). It is inadequate at small times, but gains in accuracy as the film reaches the final flat state. 4.2 Oldroyd B fluid We have investigated the non-linear response of a viscoelastic film for several values of the wave number E and the amplitude e. Typical results are shown in Fig. 4 for E = 1, c = 0.8, De = 15, and B = 0.25. For the sake of comparison, we also present the Newtonian solution and the linear predictions derived from (27) and (31). We have specified the Newtonian velocity and stress fields as initial conditions for the viscoelastic computations. As a result, the viscoelastic solution follows its Newtonian counterpart at the very
1.01
I
I
1.1-
[ = 1
, E = 0.8
1
2
0
I
3
I
I
I
4
5
6
Il.0
7
Time Fig. 4. Maximum and minimum film thicknesses as a function of time for Newtonian and Oldroyd B (De =15, /? = 0.25) fluids. Same conventions as in Fig. 3.
232
- 0.6
0
50
100
150
200
Time
Fig. 5. Maximum and minimum film thicknesses as a function of time for Newtonian and Oldroyd B (De = 195, j? = 0.25) fluids. Same convections as in Fig. 3.
beginning of the leveling process. At later times, however, elastic effects strongly retards the flow. The actual evolution of the free surface depends on the selection of initial conditions for the stress and velocity fields. It is thus irrelevant to compare at small times non-linear and linear predictions in the viscoelastic case (in the Newtonian case, there is no need for initial conditions for the velocity field since we have assumed creeping flow:) At later times, however, the non-linear results should exhibit the leveling rate predicted by the linear theory. This does indeed happen. As in the Newtonian case, the non-linear computations indicate that the thin region of the viscoelastic film levels more slowly than the thick region. The results for a disturbance of short wave length (5 = r) are shown in Fig. 5. In this case, c = 0.5, De = 195, and p = 0.25. Again, significant viscoelastic effects are predicted with the finite element analysis, and the linear theory is seen to provide an accurate estimate of the leveling rate after the initial portion of the transient. Acknowledgments The authors gratefully acknowledge stimulating discussions with Professor Morton M. Denn. This work was supported by the Director, Office of Energy Research, Office of Basic Energy Sciences, Materials Science Divi-
233 sion of the U.S. Department of Energy under Contract No. DE-AC0376SFOOO98. The numerical simulations described in this paper have been conducted on a Cray X-MP supercomputer of the National Magnetic Fusion Energy Computer Center, Lawrence Livermore National Laboratory. References 1 2 3 4 5 6 7 8 9 10 11
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
D. Degani and C. Gutfinger, Israel J. Technol., 12 (1974) 191. D.J. Coyle, Ph.D. Thesis, Univ. of Minnesota, 1984. J. Greener, T. Sullivan, B. Turner, and S. Middleman, Chem. Eng. Commun., 5 (1980) 73. T. Bauman, T. Sullivan, and S. Middleman, Chem. Fng. Commun., 14 (1982) 35. B.E. Anshus, ,Amer. Chem. Sot. (Org. Coating and Plastics Chem.), 33(2) (1973) 493. G.T. Spitz, Amer. Chem. Sot. (Org. Coating and Plastics Chem.), 33(2) (1973) 502. H. Vanoene, J. Adhesion, 4 (1972) 247. A. Quach, Ind. Eng. Chem. Prod. Res. Develop., 12.2 (1973) 110. S.E. Orchard, Appl. Sci. Res., 11 (1962) 451. D. Degani and C. Gutfinger, Comput. Fluids, 4 (1976) 149. H.S. Kheshgi and L.E. Striven, Penalty finite element analysis of unsteady free surface flows. In R.H. Gallagher, J.T. Oden, O.C. Zienkiewia, T. Kawai, and M. Kawahara (Eds.), Finite Elements in Fluids, Vol. 5, Wiley, 1984, p. 393. H.S. Kheshgi, Ph.D. Thesis, University of Minnesota, 1983. H.S. Kheshgi and L.E. Striven, Chem. Eng. Sci., 38 (1983) 525. V.M. Biermann, Rheol. Acta, 7(2) (1968) 138. R.I. Tanner, Engineering Rheology, Clarendon Press, Oxford, 1985. R. Keunings, J. Comput. Phys., 62 (1986) 199. J.G. Oldroyd, Proc. R. Sot. A., 200 (1950) 523. V.G. Levich, Physicochemical Hydrodynamics, Prentice-Hall, Englewood Cliffs, N.J., 1962. R.B. Bird, R.C. Armstrong, and 0. Hassager, Dynamics of Polymeric Liquids, Vol. 1, Fluid Mechanics, Wiley, New York, 1977. R.B. Bird, 0. Hassager, R.C. Armstrong, and C.F. Curtiss, Dynamics of Polymeric Liquids, Vol. 2, Kinetic Theory, Wiley, New York, 1977. K.P. Jackson, K. Walters, and R.W. Williams, J. Non-Newtonian FluidMech., 14 (1984) 173. MACSYMA Reference manual (version lo), The Mathlab Group, Laboratory for Computer Science, M.I.T., 1983. Jenkins, Algorithm 493, TOMS, l(1975) 178. H.M. Laun, Rheol. Acta, 17. 1 (1978) 1. J. Meissner, Pure and Appl. Chem., 42 (1975) 553. JGreener and R.W. Connelly, J. Rheol., 30 (1986) 285. R.J. Binnington and D.V. Boger, J. Rheol., 29 (1986) 887. D.W. Bousfield, R. Keunings, G. Marrucci, and M.M. Dem, J. Non-Newtonian Fluid Mech., 21 (1986) 79. M.J. Crochet, A.R. Davies, and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, 1984.
219
ANALYSIS OF SURFACE TENSION DRIVEN LEVELING IN VISCOELASTIC FILMS
R. KEUNINGS
and D.W. BOUSFIELD
Center for Advanced Materials, Lawrence Berkeley Laboratory, Berkeley, CA 94720 (USA)
University of California,
(Received June 3, 1986)
We study the surface tension driven leveling of surface irregularities in viscoelastic films deposited on a flat surface. Analytical results are presented for the generalized Maxwell model in the limit of infinitesimal surface disturbances. The results show that elasticity retards the leveling process. In some cases, the retardation effect is dramatic. Keunings’ finite element algorithm for solving viscoelastic moving boundary problems is used to analyze the case of finite-amplitude surface disturbances. Results for an Oldroyd B fluid display significant elastic effects which are consistent with the asymptotic predictions.
1. Introduction Coating operations involve small-scale free surface flows by which a solid substrate is covered with a thin liquid film. The appearance and performance of the finish greatly depend on the ability of the liquid film to level the surface irregularities generated during the coating process. In paint coating, for example, it is often observed that the film left behind the brush exhibits nearly sinusoidal disturbances (the so-called brush-marks). These disturbances are generated by the ribbing instability mechanism that takes place in the wake of coating devices (see e.g. [l-4]). Another source of surface irregularities is found in the powder coating process, whereby a multi-layer of polymer particles is deposited on a metallic substrate; upon heating, the polymer melts and covers the substrate by spreading and sintering. Surface corrugations result from local .non-uniformities of the powder deposition. The extent to which these corrugations will level during 0377-0257/87/$03.50
0 1987 Elsevier Science Publishers B.V.
220 the remainder of the bake cycle will affect the quality of the finish (see e.g. F-71). The physics and mechanics of leveling have been studied since the 1920’s (see the review by Quach [S]). The driving force is generally a combination of surface tension and gravity, although the latter can be safely neglected in many applications. Most of the work has focussed on the surface leveling of viscous, horizontal films of finite thickness and infinite extent. It is usually assumed that the amplitude of the corrugations is small compared to the mean film thickness, and that inertial forces can be neglected in comparison with viscous forces. In this context, Orchard [9] has proposed a linear theory of the surface tension driven leveling of two-dimensional disturbances in Newtonian films. The author expressed an arbitrary (but small) surface disturbance as a Fourier series, and obtained an expression relating the decay rate of each individual Fourier mode to the wave number, the equilibrium thickness of the film, the surface tension coefficient, and the fluid viscosity. Similar studies that include the effects of pseudoplasticity, thixotropy, and loss of volatiles are reviewed by Quach [8]. In view of its asymptotic character, Orchard’s analysis applies at the final stages of the leveling process. The case of finite-amplitude disturbances is difficult to analyse since it involves the solution of a complex transient free surface flow. Numerical procedures are clearly needed. Degani and Gutfinger [lo] used a version of the marker and cell method to simulate their laboratory experiments with glucose-water solutions. Good agreement with experimental data was reported. More recently, Kheshgi and Striven [ll] used a penalty finite element method for unsteady free surface flows to study the capillarity driven, viscosity dominated leveling of finite-amplitude sinusoidal corrugations (see also Kheshgi [12]). They found that the deeper region of the film levels faster than the shallower region at the early stages of the leveling process. This non-linear effect is observed in practice [13]. At later times, the leveling rate is accurately predicted by Orchard’s analysis. All these investigations have dealt with inelastic fluids. Viscoelasticity has long been recognized as an influential factor, however [8]. To our knowledge, the only attempt to determine viscoelastic effects on the leveling process dates back to 1968 [14]. The author extended Orchard’s linear analysis using the Rivlin-Ericksen second-order model. Since then, it has been established that such a model should not be used in the analysis of unsteady flows (see e.g. [15]). Re-examination of the problem is thus called for. In the present paper, we study the surface tension driven leveling of surface irregularities in viscoelastic films. First, we generalize Orchard’s asymptotic analysis in a way that is consistent with the linear viscoelastic limit of Simple Fluid Theory. To this end, we use the general linear
221 viscoelastic model with a relaxation function defined as a discrete sum of exponent&. The results show that elasticity retards the leveling process. The effect can be drastic. Numerical results for finite-amplitude corrugations are obtained with Keunings’ finite element algorithm for viscoelastic moving boundary problems [16]. We use the Oldroyd B constitutive equation, a useful model in the context of dilute polymeric solutions [17]. The results exhibit the drastic retardation effect of elasticity found with the linear theory. At the early stages of the leveling process, we predict non-linear effects similar to those observed in the Newtonian case by Kheshgi and Striven [ll]. At later times, the asymptotic theory is seen to provide an accurate measure of the leveling rate. 2. Formulation A schematic of the leveling problem is shown in Fig. 1. We assume that the flow is incompressible, isothermal, and two-dimensional. Inertial and body forces are neglected. The set of conservation laws reduces thus to -vp+v.T=O,
0)
vau=o,
(2) where p is the pressure, T is the extra-stress tensor, and 2) is the velocity vector. A constitutive equation is needed to relate the extra-stress tensor to the deformation experienced by the fluid. Our selection will be discussed later in the text. The unknown fields T, u, and p depend on time t and the space variables x and y. In the present application, the integration domain is also an unknown function of time. Its evolution is determined through the kinematic condition g
+ “$
= vu,
(3)
which relates the film thickness h to the velocity components free surface (Fig. 1).
ho
v,, vv at the
h
I //////////////////////////////// Substrate Fig. 1. Schematic of the leveling process. The liquid film extends from x = - COto x = + 00.
222 We consider the evolution of disturbances applied to the equilibrium film thickness. The surface disturbance is denoted by 6, i.e. h(x,
t) = h, + S(X, t),
(4
where h, is the thickness of the unperturbed boundary conditions apply: (i) stress balance at the free surface (y = 6):
(-pz+
T).n=a
a*h/ax*
[l + (i3h,ax)2]3'2n'
(flat) film. The following
(5)
where n is the normal to the free surface and u is the coefficient of surface tension [18]. The atmospheric pressure has been arbitrarily set to zero in (5). (ii) no-slip at the substrate (y = -h,): u, = uv = 0.
(6)
The selection of initial conditions for the film thickness and the velocity and stress fields is discussed hereafter. 3. Linear analysis We first consider a linearized version of the leveling problem valid for infinitesimal disturbances (S small). In this context, the choice of a specific rheological model is largely facilitated. We use in the present paper the general linear viscoelastic model with a relaxation function defined by a sum of N exponential terms [19]. The resulting constitutive equation is known as the generalized Maxwell model. It reads T=
&
(7)
i=l
8T. q. + r& = 2piD. The symbol D denotes the rate of strain tensor i( vo + VU=). The constants pi and 7i are viscosity coefficients and relaxation times, respectively. When all the oh’svanish, the model (7)-(8) defines a Newtonian fluid with viscosity p = C~v-l~i. Hereafter, we shall refer to this particular case as the Newtonian counterpart of the viscoelastic model (7)-(8). After linearization, the kinematic condition (3) reduces to
as/at = vy.
(9)
223 The free surface condition (5) becomes TX,,= 0, -P + T, = u
a%/ax2,
and is specified at the known location y = 0. 3.1 Analytical
solution
Following a standard Fourier procedure, we seek for each unknown TX,, TYY, TX,, ux, uv, and p, a solution of the linearized equations in the complex form
f (x, y, t) =f( y) eikxpat,
(12)
where k = 2n/A is the disturbance wave number and (Y is the leveling rate. The surface disturbance 6 is sought in the-form 6(x,
t) = S eikxpat,
(13)
where the coefficient 8 is a constant. In the context of this linearized analysis, the solution corresponding to an arbitrary (but small) thickness disturbance is obtained by linear superposition of Fourier modes (12)-(13). Inserting (12) in the momentum equation (1) and eliminating the pressure, we obtain ikcX + c;
+ k2E,
- ikqY = 0,
(14)
where the symbol ’ denotes the derivative with respect to y. On the other hand, the constitutive equations (7)-(8) yield TX = 2p.,[ik&],
05)
T,=21r&],
06)
<,,=py[C”+ikUy],
(17)
with cl,, defined as
Note that equations (15)-(17) have the form of a Newtonian constitutive model relating the Fourier coefficients of T and D, with the viscosity p replaced by IL,,_ The incompressibility constraint (2) gives
224
and with the help of (15)-(17) and (19), eqn. (14) reduces to an ordinary differential equation for C__:
- -kZ 12u,=o. d2
(20)
[ dY2 The general solution of this equation is Ei,(y)=a
cash ky+b
sinh ky+cycosh
ky+dy
sinh ky,
(21)
where a, b, c, d are constants determined by the boundary conditions. The linearized interface conditions (lo)-(11) yield ak=
b=--8,
-d,
ok
(22)
2l-G
and the no-slip condition (6) at the substrate gives c= -kb+(t). a = W(5), Here, 5 = 2?rh &I, and the functions 4, (p are defined by
(23)
4w = tanh 5 - 5 sech2t 1
HO=
+ E2 sech2t
’
1 + t2 1sech25 -
(24) (25)
Finally, use of the linearized kinematic condition (9) provides an implicit expression for the leveling rate CC
where I_L ,,( CX)is given by (18). Equation (26) is formally equivalent to the result of Orchard [9] for a Newtonian fluid, with p replaced by pV( a). 3.2 Newtonian fluid
The case of a Newtonian flu? has been analysed in detail by Orchard [9] and will serve as a reference SC ion. In order to investigate the effect of elasticity on the leveling rate, we freeze the values of the parameters (A, h,, u, p), and define the Newtonian leveling rate according to (26) (27)
For thin films (5 -C l), (27) becomes
(28)
225 TABLE 1 Leveling rate and half time as functions of the relaxation time; one-mode Maxwell model r(s) ct
(s-l)
h/Z
(4
0.
0.1
1.
10.
50.
0.088 8.
0.087 8.
0.081 9.
0.047 15.
0.016 43.
100. 0.009 77.
This result will be used later to define a time scale for the non-linear computations. 3.3 One-mode Maxwell fluid For a single Maxwell element, i.e.p = pr, 7 = 7r, and N = 1 in (7)-(8), obtain from (US), (26) and (27)
a= 1+
aNewt raNewt
-
we
(29)
Since the Newtonian leveling rate is always positive [9], the present analysis predicts that elasticity has a retardation effect on the leveling process. The following example shows that the leveling rate can be significantly reduced with highly elastic fluids. Consider a film of thickness h, = 0.01 cm and a disturbance of wave length X = 0.05 cm. Realistic numbers for a polymer melt are lo4 poise for the viscosity, and 40 dyne/cm for the surface tension coefficient. We show in Table 1 the leveling rates predicted on the basis of (27) and (29) for values of the relaxation time between 0 and 100 s. Over this range, the leveling rate experiences a ten-fold decrease from 0.088 s-r to 0.009 s-l. A useful time constant for the leveling process is the half time t,,, defined as the time required for the amplitude of the disturbance to decrease by a factor of two. It follows from (13) that $2 = a - ’ In 2. Values of this time constant are also shown in Table 1. 3.4 Jeffreys’ model If we particularize the constitutive model (7)-(8) 72 = 0, we obtain Jeffreys’ model:
with N = 2, 7 = 7i, and
where 7 * is the retardation time p2~/p and p= ~1~+ p2. Equation (30) is the linear viscoelastic limit of the dumbbell model for dilute polymeric
226 solutions [20]; ~1~and p are the solvent and solution viscosities, respectively. For this model, we find (y=
1 - rcy ~_T*aaNewt*
(31)
Equation (31) is quadratic in 0~.It can be shown that its roots are always real and positive. Clearly, the root of least magnitude will dominate the leveling process. Inspection of (31) reveals that at least one of the roots is smaller than aNewt, which implies that elasticity retards the leveling process. This result is obvious if one root satisfies 1 - TCX > 0. The case 1 - ~a -K0 and 1 - alar > 0 is impossible since (Yis positive. Now, assume that one root & is such that both 1 - T(Y(‘) and 1 - T*&) are negative. This implies that (Y(l)2 aNewt.Then, the second root (Y(*)must be smaller than aNewtsince the product of the two roots is equal to aNewJr*. For the sake of illustration, consider the polymeric solution (Boger fluid Bll) studied by Jackson et al. [21]. The material constants are in this case p = 1040 poise, T = 2.54 s, and T* = 1.97 s. We adopt the values of the previous example for the other parameters of the leveling problem. The Newtonian leveling rate is equal to 0.85 s-l, while the smallest root of (31) is 0.34 s-l. The corresponding values of the half time t,,, are 2 s for the viscoelastic fluid, and 0.8 s for its Newtonian counterpart. The retardation effect of elasticity is again significant. 3.5 N-mode Maxwell model For a discrete spectrum of N relaxation times, the implicit definition of the leveling rate (26) is a polynomial equation of degree N. As a result, the analysis becomes difficult when N >, 3 and direct computation of the roots of (26) for specific spectra { ri, pi} is called for. To do so, we have used the symbolic manipulation system MACSYMA [22] which implements an algorithm of Jenkins [23] for finding the roots of real polynomials. In ah our computations, the solutions of (26) were found to be real positive and the root of least magnitude was smaIler than the corresponding Newtonian leveling rate. We describe below two examples of such computations for values of the leveling parameters (X, h,, a) identical to those used in our previous illustrations. The results obtained with the discrete spectrum are compared to those obtained with either the largest relaxation time pi or the average relaxation time T,, defined by
(32) Our first example makes use of the relaxation data tabulated by Laun [24] for a low-density polyethylene similar to the IUPAC sample A (Meissner
227 TABLE 2 Relaxation data for LDPE Melt 1 at 150°C (Lauu [24]); c = 511000 Poise, 7av= 59 s
i
,a
pi (Poise)
1 2 3 4 5 6 7 8
1000. 100. 10 1. 0.1 0.01 0.001 0.0001
10000. 180000. 189000. 98000. 26700. 5860. 948. 129.
[25]). The 8-mode spectrum is given in Table 2. The Newtonian leveling rate is 0.0017 s-l, while the smallest root of (26) is 0.00097 s-l in the viscoelastic case. These values yield a half time t,,, of 404 s and 712 s, respectively. The one-mode analysis using the average relaxation time gives an overestimated value (0.0015 s-l) for the viscoelastic leveling rate. The largest relaxation time, on the other hand, yields a value of 0.00063 s-l which is closer to that obtained with the discrete spectrum. A more drastic retardation effect of elasticity is obtained with the data of Greener and Connelly [26] for an aqueous solution of polyacrylamide (Table 3). Here, the Newtonian leveling rate is 5 s-l and the smallest root of (26) is 0.01 s-l. The corresponding values of t,,, are 69 s for the polymeric solution and 0.14 s for its Newtonian counterpart. In this case, use of the average relaxation time yields a good approximation (0.06 s-‘) of the
TABLE 3 Relaxation data for au aqueous solution of polyacrylamide (Greener and Couuelly [26]); p = 168 Poise, T,, = 16 s
i
Ti
1 2 3 4 5 6 7 8 9 10 11
100. 32. 10. 3.2 1. 0.32 0.1 0.032 0.01 0.001 O.O@Ql
ls)
pi (Poise) 7.6 44. 41. 34.4 20.6 11.8 4.8 2.42 0.79 0.845 0.0254
228 leveling rate. The result obtained on the basis of the maximum relaxation time is even better (0.01 s-l). The leveling process appears to be dominated by the largest relaxation time, at least in the two examples discussed here. 4. Non-linear fiite
element analysis
In view of its asymptotic character, the above analysis applies at the final stages’of the leveling process where the fluid is close to the state of rest. It is quite general in the sense that the viscoelastic model (7)-(g) constitutes a good approximation of any Simple Fluid in the linear viscoelastic limit. The issue of finite amplitude disturbances of the film thickness is of course more difficult to address in a general way. Indeed, the results will be specific to a constitutive model (which must be objective), and a set of initialconditions. Moreover, analytical developments are impossible since the flow problem becomes highly non-linear. These difficulties are typical of the analysis of complex viscoelastic flows. We describe in this section the results of a non-linear finite element analysis of the leveling process. We have selected the Oldroyd B model to define the rheology of the film: T+G=2p
D+r*vD [
. I
(33)
Here, the symbol v denotes the upper-convected derivative. For example,
V T=g+oe~~-~.~~--vuT.~.
(34)
The material constants EL,7, and 7 * have the same meaning as in (30). The Oldroyd B model is a generalization of Jeffreys’ equation for arbitrary deformations and can be derived from the dumbbell kinetic theory [20]. It has proved quite useful in the modeling of dilute polymeric solutions (see e.g. [27,28]). We consider in the present section the evolution of a finite amplitude disturbance applied to the film kness. The initial disturbance is sinusoidal with a wave length A. We assume that the free surface remains spatially periodic as it levels (it will of course depart from a sinusoidal shape because of non-linear effects.) In addition, we restrict our attention to disturbances that are symmetric about the plane x = X/2 (Fig. 2). This assumption translates into the boundary conditions u
x
_au,=ah=, ax
ax
’
(35)
specified at the crest (x = 0) and trough (x = X/2) of the free surface. With
229
X 0.
x/2
Fig. 2. Typical geometryand finite element mesh for the finite-amplitudecomputations.
this simplication, the integration domain extends from x = 0 to x = X/2. The mathematical problem now consists of solving the governing equations (l)-(3) and (33) in terms of T, u, p, and h, with the boundary conditions (S)-(6) and (35) and proper initial conditions. Thus formulated, the leveling process is a two dimensional, transient, free surface viscoelastic flow driven by surface tension. The numerical results described hereafter have been obtained by means of the finite element algorithm developed by Keunings [16] for the solution of viscoelastic moving boundary problems. Briefly, the method consists of a Galerkin principle invoked on deforming finite elements. The motion of the nodes is anchored to the displacement of the free surface in a way that prevents gross distortion of the finite element layout. The Gale&in formulation takes this nodal motion into account and contains the standard Eulerian and Lagrangian approaches as particular cases. A semi-implicit predictorcorrector scheme is used for the temporal integration. This technique has been applied successfully in the analysis of ‘the surface tension driven breakup of viscoelastic jets [16,28]. The flow domain is discretized by means of nine-node Lagrangian elements. A typical finite element mesh is shown in Fig. 2. The free surface representation h is approximated by Co - P* one-dimensional elements; we use Co - P’ interpolations for the pressure and the extra-stress tensor, and
230 Co-P2 basis functions for the velocity. The initial free surface shape is set to h(x, 0) =h,[l
+r cos(2Izx/A)],
(36)
where c is the dimensionless amplitude of the disturbance. The other three dimensionless groups that arise in the present application are the Deborah number De, the dimensionless wave number 5, and the ratio /3of characteristic times:
De=7
16m4ahfj 3jLx4
’
t=,,
2nh,
(37)
In our definition of the Deborah number, we have used the result (28) to obtain a typical time scale of the flow. 4. I Newtonian fluid
The evolution of finite-amplitude disturbance in Newtonian films has been studied by Kheshgi and Striven [ll] (see also [12]). The authors used a penalty finite element method for solving the moving boundary problem. The results for creeping flow indicated that the thick region of a sinusoidally disturbed film levels faster than the thin region. This non-linear effect
1
.c
1:
0.9
%
%
2 y 1.2
0.8 E 2 l-
z
I-
E
.EI1.3
0.7
4
1.4
0.6
1.5
0.5
E i ‘E .I
Time
Fig. 3. Maximum and minimum film thicknesses as a function of time (Newtonian fluid); finite-amplitude computations: -, linear theory: - - - - - -, Kheshgi’s results [12]: l ; time ant thicknesses are made dimensionless with 3~A4/16a4uh i and h,,, respectively.
231 deforms the initially sinusoidal shape of the free surface. In the particular case of creeping Newtonian flow, Keunings’ algorithm contains the classical displacement (u-u-p) and mixed Gale&in formulations of the Stokes problem [29]. We have repeated with both methods a simulation of Kheshgi [12] for 5 = r/10 and c = 0.5. Figure 3 shows the evolution of the maximum and minimum film thicknesses as they approach the steady state h = h,. The results obtained with the displacement method are indistinguishable from those of the mixed technique at the scale of the drawing. The agreement with Kheshgi’s computations is excellent. Also shown for comparison is the linear prediction derived from (13) and (27). It is inadequate at small times, but gains in accuracy as the film reaches the final flat state. 4.2 Oldroyd B fluid We have investigated the non-linear response of a viscoelastic film for several values of the wave number E and the amplitude e. Typical results are shown in Fig. 4 for E = 1, c = 0.8, De = 15, and B = 0.25. For the sake of comparison, we also present the Newtonian solution and the linear predictions derived from (27) and (31). We have specified the Newtonian velocity and stress fields as initial conditions for the viscoelastic computations. As a result, the viscoelastic solution follows its Newtonian counterpart at the very
1.01
I
I
1.1-
[ = 1
, E = 0.8
1
2
0
I
3
I
I
I
4
5
6
Il.0
7
Time Fig. 4. Maximum and minimum film thicknesses as a function of time for Newtonian and Oldroyd B (De =15, /? = 0.25) fluids. Same conventions as in Fig. 3.
232
- 0.6
0
50
100
150
200
Time
Fig. 5. Maximum and minimum film thicknesses as a function of time for Newtonian and Oldroyd B (De = 195, j? = 0.25) fluids. Same convections as in Fig. 3.
beginning of the leveling process. At later times, however, elastic effects strongly retards the flow. The actual evolution of the free surface depends on the selection of initial conditions for the stress and velocity fields. It is thus irrelevant to compare at small times non-linear and linear predictions in the viscoelastic case (in the Newtonian case, there is no need for initial conditions for the velocity field since we have assumed creeping flow:) At later times, however, the non-linear results should exhibit the leveling rate predicted by the linear theory. This does indeed happen. As in the Newtonian case, the non-linear computations indicate that the thin region of the viscoelastic film levels more slowly than the thick region. The results for a disturbance of short wave length (5 = r) are shown in Fig. 5. In this case, c = 0.5, De = 195, and p = 0.25. Again, significant viscoelastic effects are predicted with the finite element analysis, and the linear theory is seen to provide an accurate estimate of the leveling rate after the initial portion of the transient. Acknowledgments The authors gratefully acknowledge stimulating discussions with Professor Morton M. Denn. This work was supported by the Director, Office of Energy Research, Office of Basic Energy Sciences, Materials Science Divi-
233 sion of the U.S. Department of Energy under Contract No. DE-AC0376SFOOO98. The numerical simulations described in this paper have been conducted on a Cray X-MP supercomputer of the National Magnetic Fusion Energy Computer Center, Lawrence Livermore National Laboratory. References 1 2 3 4 5 6 7 8 9 10 11
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
D. Degani and C. Gutfinger, Israel J. Technol., 12 (1974) 191. D.J. Coyle, Ph.D. Thesis, Univ. of Minnesota, 1984. J. Greener, T. Sullivan, B. Turner, and S. Middleman, Chem. Eng. Commun., 5 (1980) 73. T. Bauman, T. Sullivan, and S. Middleman, Chem. Fng. Commun., 14 (1982) 35. B.E. Anshus, ,Amer. Chem. Sot. (Org. Coating and Plastics Chem.), 33(2) (1973) 493. G.T. Spitz, Amer. Chem. Sot. (Org. Coating and Plastics Chem.), 33(2) (1973) 502. H. Vanoene, J. Adhesion, 4 (1972) 247. A. Quach, Ind. Eng. Chem. Prod. Res. Develop., 12.2 (1973) 110. S.E. Orchard, Appl. Sci. Res., 11 (1962) 451. D. Degani and C. Gutfinger, Comput. Fluids, 4 (1976) 149. H.S. Kheshgi and L.E. Striven, Penalty finite element analysis of unsteady free surface flows. In R.H. Gallagher, J.T. Oden, O.C. Zienkiewia, T. Kawai, and M. Kawahara (Eds.), Finite Elements in Fluids, Vol. 5, Wiley, 1984, p. 393. H.S. Kheshgi, Ph.D. Thesis, University of Minnesota, 1983. H.S. Kheshgi and L.E. Striven, Chem. Eng. Sci., 38 (1983) 525. V.M. Biermann, Rheol. Acta, 7(2) (1968) 138. R.I. Tanner, Engineering Rheology, Clarendon Press, Oxford, 1985. R. Keunings, J. Comput. Phys., 62 (1986) 199. J.G. Oldroyd, Proc. R. Sot. A., 200 (1950) 523. V.G. Levich, Physicochemical Hydrodynamics, Prentice-Hall, Englewood Cliffs, N.J., 1962. R.B. Bird, R.C. Armstrong, and 0. Hassager, Dynamics of Polymeric Liquids, Vol. 1, Fluid Mechanics, Wiley, New York, 1977. R.B. Bird, 0. Hassager, R.C. Armstrong, and C.F. Curtiss, Dynamics of Polymeric Liquids, Vol. 2, Kinetic Theory, Wiley, New York, 1977. K.P. Jackson, K. Walters, and R.W. Williams, J. Non-Newtonian FluidMech., 14 (1984) 173. MACSYMA Reference manual (version lo), The Mathlab Group, Laboratory for Computer Science, M.I.T., 1983. Jenkins, Algorithm 493, TOMS, l(1975) 178. H.M. Laun, Rheol. Acta, 17. 1 (1978) 1. J. Meissner, Pure and Appl. Chem., 42 (1975) 553. JGreener and R.W. Connelly, J. Rheol., 30 (1986) 285. R.J. Binnington and D.V. Boger, J. Rheol., 29 (1986) 887. D.W. Bousfield, R. Keunings, G. Marrucci, and M.M. Dem, J. Non-Newtonian Fluid Mech., 21 (1986) 79. M.J. Crochet, A.R. Davies, and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, 1984.