I n d . Eng. Chem. Res. 1988,27, 1313-1320
1313
From Molecular Models to the Solution of Flow Problems? Lewis E. Wedgewood* and R. Byron Bird Rheology Research Center and Department of Chemical Engineering, University of Wisconsin, Madison, Wisconsin 53706
One of the goals in polymer fluid dynamics is to use kinetic theory to derive a constitutive equation for a polymeric liquid starting from a molecular model, to use the constitutive equation to solve flow problems, and finally to use kinetic theory to describe molecular stretching and orientation. We illustrate this procedure by using a finitely extensible, nonlinear, elastic dumbbell with the Peterlin approximation (FENE-P dumbbell) as a crude model of a polymer molecule and then derive the constitutive equation. A number of new results for the FENE-P dumbbell are given which include an analytic expression for the elongational viscosity, numerical calculations of the stress growth after the start-up of elongational flow, and the third-order fluid constants. The technique of “model matching” is illustrated for squeezing flow between circular disks, the Weissenberg rod-climbing effect, and the torque on a rotating sphere. Conclusions can be drawn about the average extension and alignment of the dumbbell in these three flow fields. 1. Introduction One of the basic problems in polymer fluid dynamics is the development of a constitutive equation for the stress tensor, 7 , which is necessary for solving flow problems. Constitutive equations are developed in three ways: (1) The equation is chosen empirically to fit experimental data as well as possible. (2) Continuum mechanics can be used to generate “reduced equations” appropriate for specific classes of flows. (3) Molecular theory can lead to a constitutive equation in terms of the parameters of the mechanical model. In this last method, macroscopic fluid behavior is connected with the microscopic behavior of the constituent polymer molecules. Of course, the quality of the connection depends on the realism of the mechanical model. Unfortunately, many of the more realistic molecular models lead to such cumbersome expressions for the stress tensor that they are not used often for solving flow problems. In this paper we illustrate the process of developing a constitutive equation from kinetic theory and using the result to solve flow problems for a dilute solution of FENE dumbbells (Warner, 1972). We choose this simple model because it contains the basic characteristics of molecular stretching, orientation, and finite extensibility seen in polymer molecules, yet it is relatively simple mathematically. An approximate constitutive equation is derived which gives material functions that compare well with rheometric data. In addition, an approximate configurational distribution function is also derived, from which information about the orientation and extension of the dumbbells can be determined for a given flow field. Flow problems are solved by taking over published solutions developed from constitutive equations that are good for a wide class of fluids and a restricted flow regime and then matching their parameters with those of the FENE dumbbell model (Figures 6.2-2 of Bird et al. (1987a)).
2. The Diffusion Equation We model the polymer molecule by an elastic dumbbell (p 55 of Bird et al. (1987b, 1980)) consisting of two “beads” with equal friction coefficients, C. The beads are connected given by by a nonlinear spring with a connector force F(C) F(C)= f(Q2)Q (1) Here f is a scalar function of Q2 = Q-Q, where Q = r2A preliminary version of this material was presented by us in Wedgewood and Bird (1986). 0888-5885/88/2627-1313$01.50/0
r1 is the bead-to-bead vector, and r1and r2are the bead position vectors. The velocity field of the solvent is considered to be homogeneous: v = vo + [KT], where vo and K = (Vv)’ are in general functions of time, and r is the position vector; the superscript indicates the transpose of a tensor. Of interest is the configurational distribution function, $(Q,t), which is defined such that $(Q,t) dQ is the probability that an arbitrary dumbbell in the solution is in a configuration between Q and Q + dQ at time t. To find the configurational distribution function, $(Q,t), we use the equation of continuity in the configuration space
*
-= a+
-($.[[Bll$)
at where the symbol [ [ I] denotes a momentum-spaceaverage. The quantity [ [Q]] can be found by writing a force balance for each bead and then taking the difference of the two equations; when the mass of the inertia of the beads is neglected, this yields a balance of the forces for the internal motion
-C([[Qll
-
a aQ
[K-QI) - 2 k T - In $ - 2f(Q2)Q
= 0 (3)
where k is Boltzmann’s constant and T is the absolute temperature. The forces considered are hydrodynamic drag, Brownian, and intramolecular. Combining eq 2 and 3 by eliminating [ [Q]] produces the “diffusion equation” for $(Q,t):
Equation 4 has been solved for Hookean dumbbells (Lodge and Wu, 1971; van Wiechen and Booij, 1971), where f(Q2) = H a n d H i s a spring constant independent of the beadto-bead vector. Unfortunately, the Hookean dumbbell model is of limited applicability since it predicts a constant value for steady shear viscosity and an infinite elongational viscosity at a finite elongation rate. The Hookean dumbbell model can be improved upon by introducing a nonlinear connector force such as the FENE spring (Warner (1972) and p 16 of Bird et al. (198713)) which has a configuration-dependent “spring constant” f(Q2) = H/[1 -
(Q/QO)~I
(5)
where Qo is the maximum extension. Equation 4 cannot 0 1988 American Chemical Society
1314 Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988
be solved analytically with this configuration-dependent spring constant. Therefore, we replace the variable Q2 in the FENE "spring constant" eq 5 by its average ( Q2) which depends on time but not the configuration of the dumbbell. This substitution yields
f ( ( Q 2 ) ) = H/[1
-
((Q/Qd2)1 E HZ
(6)
where the angular brackets denote a configuration-space average. For an arbitrary second-order tensor B, the configuration-space average is defined as ( B ) = .fB+ dQ. We refer to F(C) = f((Q2))Qas the FENE-P connector force, and the elastic dumbbell model which employs it we refer to as the FENE-P dumbbell model. The "P" stands for Peterlin (1966), who was the first to suggest this kind of approximation. By combining eq 6 with eq 4, we get the diffusion equation for the FENE-P dumbbell model
which has the form of a multivariate linear Fokker-Planck equation where the coefficient tensors are in braces {I. Equation 7 has been rearranged, and the coefficient tensor containing f ( ( Q 2 ) ) can be brought in front of the spatial derivative since it is independent of the configuration of the dumbbell. Ottinger (1986) pointed out that eq 7 is solvable by a Gaussian ansatz (Kolmogoroff, 1931; van Kampen, 1981): +(Q,t) = 1 / ( ( 2 ~ k T / H (det ) ~ a)j'I2 exp(-(H/2kT)(a-':QQ){ ( 8 ) which leads to an ordinary differential equation for the tensor a: AHa(l)
= 6 - Za
a = (b/(b
+ 3))G
at t = -a
(9)
where AH = {/4H is a time constant. Subscript (1)on the tensor a indicates a "convected derivative"; a(1)= Da/Dt - (pa + The quantity 2 which is defined by eq 6 can be expressed in terms of the tensor a as .)'.sa
Z = (1- (l/b)tra)-'
(10)
where the tr denotes the trace of a tensor and b = HQ:/kT. The Hookean dumbbell configurational distribution function (p 72 of Bird et al. (1987b)) is recovered in the limit Qo m or b m. It should also be pointed out that, in the FENE dumbbell model of eq 5, the length of the dumbbell is restricted by Q2 < Qi2, whereas in the FENE-P dumbbell, as a consequence of the approximation in eq 6, it is the average length which is restricted by ( Q 2 )
- -
< Qo2.
Although no general solution for eq 9 is known, nonequilibrium solutions can be found, for example, for steady shear flow (u, = yy, uy = u, = 0); the tensor a in rectangular coordinates is given by a =
(;+2p2
(11)
where p = AH+/Z. This result along with eq 8 completely defines the probability, $(Q,t) dQ, of an orientation Q in steady shear flow and can be used in calculating other averages. For steady elongational flow (u, = -(1/2)2x, uy = -(1/2)ky, u, = gz), the tensor LY in rectangular coordinates is given by
which, along with eq 8, can be used to calculate configuration-space averages for elongational flow. 3. The Stress Tensor An expression for the polymer contribution to the stress tensor T~ can be derived in terms of averages of the momentum fluxes due to bead motion and spring tension (p 64 of Bird et al. (1987b)) T~
= -n(QF("))
+ nkTG
(13)
where n is the number density of polymer molecules. By substituting the FENE-P connector force into eq 13, we have
The polymer contribution to the stress tensor can be alternatively expressed by using a = ( Q Q ) ( H / k T ) (this relation can be easily shown by using eq 8 and performing the average) and eq 10 to obtain rp = -nkT(Za - 6) (15) By combining eq 15, the convected derivative of eq 15, eq 10, and eq 9, we can eliminate a and derive an expression for the constitutive equation for the FENE-P dumbbell: 7
=
-17sY(l)
+ 7p
Z = 1 + (3/b){1- [(trrP)/3a])
(164
(16~)
where y(lj = K * + K = (Vv) + (Vv)*. The constitutive equation contains four parameters: vS, the solvent viscosity; AH = c/4H, a characteristic relaxation time; b = HQo2/kT,a ratio of the elastic spring energy to the thermal energy; and a = nkT. Although these four quantities do have molecular interpretations, we treat all of them as adjustable constants to be determined from rheometric data. That is, the molecular theory is used to suggest a useful form for the Constitutive equation, and experiments are used to determine the values of the constants. As a result, the calculations of molecular stretching in this paper admit only a qualitative interpretation. It should be noted that the constitutive equation given by eq 16 may be considered as an approximate constitutive equation for the exact FENE dumbbell (i.e., without the Peterlin approximation). Indeed, eq 16 was first derived by Tanner (1975) by approximating the configurational distribution function by a delta function, +(Q,t) = S(Q S(t)), where S(t) is an arbitrary vector. 4. Material Functions In a steady shear flow (u, = qy,uy = u, = 0), it is possible to measure the shear stress and the first and second normal stress differences. The viscosity 17, the first-normal stress coefficient 'kl, and the second-normal stress coefficient 'k2 are defined by 7,y
=
-v(?)?
(17)
7xx
-
TJy
=
-'kA+)V
(18)
74y
-
72,
=
-*A?)+*
(19)
Ind. Eng. Chem. Res., Vol. 27, No. 7 , 1988 1315
ii
l
REDUCED SHEAR RATE P=[?l,?,My/RT
Figure 1. Intrinsic viscosity data versus reduced shear rate ( p = [ & o f l ~ / R Tfor ) monodisperse polystyrene in 6-solvents [(X) Nada et al. (1968) and (0,0, A, 0)Manke and Williams (198511 and calculated curves for the FENE-P dumbbell from eq 20 for b = 50, 100, and 1000.
'
-102
l
'
.
l
,
I
-IO -10-2 10-2 100 ELONGATION RATE X,i
,
I
102
. I 104
Figure 2. Elongational viscosity versus dimensionless elongation rate for FENE-P dumbbell (dotted curves) and FENE dumbbell (solid curves) for b = 10, 50, and 100.
= 2(v - 77J2/a
(21)
=0
(22)
7j
sinh [ (1/3) arcsinh (A'?) J (20)
q - qs = 3q0(A'?)-'
*2
in which qo = uAH[b/(b + 3)] and A' = (27/2)'/'[b/(b + 3)3i2]AH. A comparison of 7 and \kl for the FENE-P dumbbell, eq 20 and 21, with the numerical calculations 71 and *lfor the FENE dumbbell (i.e., without the Peterlin approximation introduced in eq 6) has been made by Fan (1985). Fan concluded that there is a semiquantitative agreement between the FENE-P dumbbell and the exact FENE dumbbell for b 3 10, and therefore, the use of the Peterlin approximation is justified in steady shear flow. In Figure 1, eq 20 is compared with the intrinsic viscosity data for monodisperse polystyrene in 8-solvents taken from N6da et al. (1968) and Manke and Williams (1985). The FENE-P dumbbell model does not contain excluded volume effects and one would expect it to have more success modeling dilute polymer solutions in @-solvents. Surprisingly, eq 20 has been more successful in comparisons with polymer solutions in good solvents (p 91 of Bird et al. (1987b)). For large shear rates, the material functions for the FENE-P dumbbell model can be approximated by simple power-law expressions. Taking the limit of eq 20 and 21 as 9 m, the material functions become
-
my-'
(23)
= m'W2
(24)
q - qs = *l
in which
m = aAH((b/2)'/2/AH)2/3 m' =
-
In a steady elongation flow (u, = -(1/2)ix, vy = -(1/2)iy, iz), it is possible to measure one normal stress difference. The elongational viscosity, i j , is defined by TZZ- 7,, = -q(i)i (25) where positive i corresponds to elongation and negative & corresponds to biaxial stretching. Again, there are three solutions because of the Peterlin approximation, which is used to derive the constitutive equation. Two are eliminated on physical grounds to give
where 9 = + [ ( l / 2 ) ( ~ ( l ~ : ~ ( l ~ ) J ' / 2 . For the FENE-P dumbbell model, the material functions can be calculated from the constitutive equation, eq 16. Both q and \kl are solutions to cubic equations with three roots, two of which can be eliminated on physical grounds, giving
*I
o
-104
n = 1/3
~ U A H ~ ( ( ~ / ~ ) ~ / ' / X , ) ~ / n'= ~
2/3
Equations 23 and 24 are appropriate for fast, steady shear flows. Mochimaru (1981) calculated numerically the shear stress and the first normal stress difference for the start-up of steady shear flow for the FENE-P dumbbell model. Most remarkable is the prediction of stress overshoot which has been observed experimentally for melts and concentrated polymer solutions (pp 118-123 of Bird et al. (1987a)).
u, =
= 377,
+ 6UA~[b/(b+ 3)]P'/' COS [(1/3) arCCOS (4P-3i2)+ CY]* W / ( W - 1)[(5A~i - W)/6x~i] (26)
in which w = ( b + 3)/b XHk
1
+ 4w)/AHi + ( 6 +~w2)/(AH$)2]
+ 12w)/XHi + (54 + 9W + 6W2)/(A~&)2 - (9W' + W 3 ) / ( A ~ k ) 3 )
4 = (1/216)( L 1 r- (w3 5
- (90
In Figure 2 the above expression for the FENE-P dumbbell is compared with numerical calculations of the elongational viscosity for the FENE dumbbell. The percent difference between for the FENE and the FENE-P dumbbell models is less than 1% for reasonable values of b (i.e., b > 50). For the start-up of steady elongational flow, the tensor a ( t ) was calculated from eq 8 using a fourth-order Runge-Kutta method. The time evolution of 7p was calculated by combining a ( t )from eq 8 with the expression for the stress tensor given by eq 15. The dimensionless elongational stress growth viscosity ($ - 3.rls)/(qm+ - 377,) for the FENE-P dumbbell is shown in Figure 3 for the dimensionless elongation rates of AHi = 0.1, 1.0, and 10.0. For the smallest elongation rate shown in Figure 3, the quantity (ij+ - 317s)/(5im+ - 3vS)increases monotonically with a monotonically decreasing slope, a result which is similar to that for the Hookean dumbbell in the same flow (p 75 of Bird et al. (1987b)). For the intermediate elongation rate of AH: = 1.0, the quantity ($ - 3q,)/(ijm+- 377,) for the FENE-P dumbbell is quantitatively similar to that of the Hookean dumbbell at short times (Le., XHt < 2.5); however, for longer times, the stress growth function for the Hookean dumbbell continues to grow without limit, whereas for the FENE-P dumbbell the same curve reaches
1316 Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 I
,
,
I
Note that b2 < bll and that b3 > b12 > blZll; these inequalities have been found for almost all molecular models studied so far (Bird, 1985). (3) The third example is the general linear viscoelastic model (p 255 of Bird et al. (198713))for flows with small displacement gradients
I
in which G(s)is the relaxation modulus and s = t - t'with t' < t . For the FENE-P dumbbell, G ( s ) is given by IO 20 30 40 DIMENSIONLESS TIME t / X ,
G(s) = -nkT expi-([b
50
Figure 3. Dimensionless elongational stress growth viscosity (q+ 3q8)/(qm+- 37,) versus dimensionless time (t/X,) for the FENE-P dumbbell with b = 50; dimensionless elongation rates hHi = 0.1, 1.0, and 10.0; and asymptotic (q-+ - 37,)/(3q0- 3qs) = 1.0, 16.7, and 31.7.
an asymptote. We also point out that the stress growth curve in Figure 3 for XHi = 1.0 lies below the curve for XHi = 0.1 because of the normalization (qm+ - 317J which is nearly 17 times larger at XHi = 1.0 than at X H i = 0.1. For the largest elongation rate shown in Figure 3, X H i = 10.0, the stress growth curve approaches a step function.
5. "Model Matching" the Constitutive Equation with Restricted Continuum-Mechanics Expressions Continuum mechanics provides several expressions for T that apply to a very wide class of viscoelastic fluids, but only within certain restricted flow regimes. Three examples are as follows: (1)The first example is the Criminale-Ericksen-Filbey (CEF) equation (p 503 of Bird et al. (1987b)) for steady shear flow: 7- = - d + ) Y ( l ) + (1/2)*1(+h(2) - ~ A + ~ ~ Y ( l ) ~ Y ( l ) I (27) For the FENE-P dumbbell, a(+), \k,(y) and \k2(+) are given by eq 20, 21, and 22, respectively, or eq 23 and 24 for large shear rates. If the terms containing \kl and \k2 are omitted from the CEF equation, then eq 27 reduces to the generalized Newtonian fluid equation. (2) The second example is the retarded-motion expansion (p 295 of Bird et al. (1987b) for slow flows, slowly varying in space and time: T = -[blY(1) + bzY(z) + bii(Y(if~(i)J+ b37(3) + b12(Y(I).Y(Z)+ Y(Z)'Y(l)I + bl:ll(Y~l):Y(l))Y(l~ + *-I (28) in which b,, b2, b,,, b3, etc., are constants and should not be confused with the model parameter b. The first term is linear in the velocity gradient, and bl is the zeroshear-rate viscosity a,,. The two terms containing b2 and bll contain terms that are quadratic in the velocity gradients, and the constants are related to the zero-shear-rate normal stress difference coefficients by \kl,o = -2bz and \kz,o = bll. Together these three terms form the second-order fluid. The terms shown in eq 28 containing b3, b12, and bl:ll are cubic in the velocity gradient and together with the previous terms form the third-order fluid. The constants corresponding to the FENE-P dumbbell are found by substituting eq 28 into eq 16b and matching like terms. This procedure leads to relations between the retarded motion constants and the FENE-P parameters given by bi = IS' + aXH[b/(b + 311 (29)
+ 3)2] bl1 = 0 b3 = a X ~ ~ [ b ' / ( b+ 3)3] = -(b + 3)bl.ll biz = 0 b2 = -aXH2[b2/(b
(30) (31)
+ 3]/b)(S/XH)j
(33)
This result is similar to that for the Hookean connector since the dumbbell is deformed only slightly by the flow. The three expressions for the stress tensor discussed in this section apply within the specified flow regimes. Flow problems that have been solved using one of these reduced expressions for the stress tensor can have the solution taken over by matching the parameters, as has been demonstrated above for the FENE-P dumbbell model. In the next section, direct solutions are considered. 6. Direct Solutions to Flow Problems Before employing the model-matching technique, we would like to point out two direct solutions by Mochimaru (1981, 1983) using the FENE-P model, eq 16. (1) Fast Squeezing Flow between Two Circular Disks. For the squeeze flow problem, Mochimaru made approximations consistent with the fast squeezing limit (i.e., the limit in which viscoelastic effects become pronounced), thereby simplifying the equations so that an analytical solution could be found (see eq 35 of Mochimaru (1981)). The results of this analysis compare well with the experimental data only in the fast squeezing regime. (2) Start-up of Couette Flow. Mochimaru gives an approximate solution for the problem of start-up of plane Couette flow using the FENE-P constitutive equation. The most interesting result of this analysis is the prediction of velocity overshoot. Newtonian fluids monotonically approach a linear velocity profile; however, the FENE-P model predicts that the transient velocity will exceed the steady-state velocity for certain parameter values and shear rates. Velocity overshoot has been recently observed for solutions of poly(isobuty1ene) solutions in decalin in the start-up of tangential flow between coaxial cylinders (pp 166-173 of Burdette (1987)).
7. Solutions to Flow Problems by the "Method of Model Matching" In section 5, it was shown how to relate reduced expressions of the stress tensor from continuum mechanics to the FENE-P constitutive equation. In this section we give three examples of this procedure for squeezing flow, rod climbing, and flow near a spinning sphere. (a) Fast Squeezing Flow between Circular Disks. The flow considered here is that of isothermal, time-dependent radial flow of a viscoelastic fluid between two disks after a constant force F has been applied normal to the disks of radius Rd and initial separation 2ho. To determine values for the model parameters, we fit experimental data for steady shear flow. Figure 4 shows Leider and Bird's (1973) data on the steady shear viscosity and first-normal stress difference for a solution of 5.2% PIB in 20.3% Vistanex LM-MS and 74.5% Oppanol B-1 along with the curves calculated from eq 23 and 24 where the parameters m and m'were found by using the least-squares method. If the solvent viscosity is neglected and a rea-
Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 1317
IO2 10-4
10-3
10-2
10-1
IO
I
SHEAR RATE 9 (s- I
io2
IO
Figure 4. FENE-P dumbbell material functions eq 23 and 24 rheometric data for PIB solution (Leider and Bird, 1973).
sonable value for b is chosen (Le., setting qs = 0 and b = 50), the remaining model parameters can be determined from the expressions for m and m‘in eq 23 and 24: XH = 8.4 s, a = 28 Pa. McClelland and Finlayson (1983) analyzed the squeeze flow problem based on the CEF equation although the authors did not state that they used the CEF equation. (In McClelland and Finlayson (1988), a paper which was published after the manuscript for this paper had been submitted, the authors reiterate the derivation of the squeeze flow problem and acknowledge their use of the CEF equation.) They found an expression for the dimensionless applied force F* = FA”/(xRd2m)using power law expressions for 17 and 9,and setting \k2 = 0. A more consistent analysis of the problem in which the effects of elasticity on the flow field are taken into consideration is given in example 10.4-1 of Bird et al. (1987a) where the expression for F * is given by
F* =
(y)n( 2) (n(-dh*/dt*)” + 3)h*(2n+1)+
yy(2)
” (-dh*
2(
/dt*)“’
[(n’+ 2) (n’ + 2)h*2n’
+ K ] (34)
where h* = h/ho, t* = t/X, and X = (m’/2m)1/(n’-n). Equation 34 differs from the result of McClelland and Finlayson by an additional factor of [(n’+ 2) + K ] in the second term where (3n + n’ + 1)- n(2n l)(n’ + 4) K= ( n + n’)(2n + n’+ 1)
+
and \k2 = 0. The first term in eq 34 corresponds to the Scott equation which is derived by using the generalized Newtonian fluid equation (p 190 of Bird et al. (1987a)). The power-law parameters in eq 34 were chosen to be those determined from experimental data for the FENE-P dumbbell model in eq 23 and 24. For the FENE-P dumbbell, we have 2n = n’so that the quantity -dh*/dt* is the solution of a quadratic equation. This result can then be used in the integral 1
tlj2/X = ij2(-dh*/dt*)-’ dh*
IO2
GEOMETRIC RATIO Rlh 0
)
(35)
which is the time required to squeeze out half of the liquid between parallel disks under a constant force F *. Equation 35 is evaluated by using Simpson’s rule of integration with variable quadrature. Figure 5 shows Leider and Birds data for the dimensionless time tl12/X versus the initial radius to height ratio Rd/ho. For values of tl12/X < 1, Figure 5 clearly demonstrates the inadequacy of the Scott equation to describe the data. This deficiency of the Scott equation is caused in part to the fact that it ignores the
Figure 5. Squeezing flow data (Leider and Bird, 1973) along with FENE-P dumbbell predictions based on rheometric data from Figure 4.
effects of elasticity. The solid curves in Figure 5 are generated from eq 34 for which all the parameters were determined from the steady shear flow data in Figure 4. By including elastic terms, there is a great improvement in describing the squeeze flow data. Finally, we write down an expression for the average stretching of the dumbbells by combining eq 10, 16c, 27, and 24
( Q 2 ) / Q o 2= 1 - (1+ (3/b)(l where E = (2/3)XH2{(b/2)1/2/XH)4/3. proximately
+ Ei,2/3))-1
(36)
The shear rate is ap-
where r and z are the radial and axial coordinates with the origin at the center of the disks. Since the parameters in eq 36 were determined by curve fitting with steady shear flow data, eq 36 is only a qualitative expression for molecular stretching. Armstrong et al. (1980) have used a similar relation for molecular stretching with qualitative success. From eq 36 and 37, it can be seen that the dumbbell end-bend distance increases with shear rate and is greatest near the surface of the disk, h = ho, and at the exit, r = Rd. The orientation of the dumbbells is described by the configurational distribution function eq 8. The covariance tensor is related to the stress by eq 15 and 16c where T~ = 7 for t s= 0. An expression for ‘I is given by the CEF equation, eq 27, where we again use the power law material functions, eq 23 and 24. The appropriate kinematic tensors are formed by using the velocity flow field (v = uJ,+ u,6J given by eq 10.4-48 and -49 of Bird et al. (1987a). When these results are combined, $(Q) for fast squeezing flow is $(Q) = l / ( ( 2 ~ k T / H (det )~ ( a , ; l ~ cos2 ~ 4 2a,;l~~ cos 4
+
exp(-(H/2kT) X + a 8 ; 1 ~sin2 2 4 CY,;~Z~)) (38)
+
-
where Q = iisr + SZand cos 4 is the dot product of the unit vectors 6, and 6,. In the limit r 0 3 , $(Q) = @(P)6(4)6(Z) where @(P) is Gaussian in P; in other words, the dumbbell is aligned with the direction of flow and has a Gaussian distribution of end-to-end distances. (b) The Weissenberg Rod-Climbing Experiment. When a vertical rod is made to rotate in a beaker of Newtonian fluid, the free surface of the fluid dips in the vicinity of the rod. Weissenberg was the first to explain the rise of the free surface characteristic of certain nonNewtonian fluids. Joseph and Fosdick (1973) and Joseph et al. (1984) have developed a theory for a second-order fluid describing to a good approximation the height of the
1318 Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 Table I. GrouDinP ( b , - 2 b l l ) from Molecular Theory and from Rod Climbing steady shear flow constants for eq 20 and 21 Dolvmer concn, ?& vo, Paas A”, 9 a , Pa 3.58 4.5 0.13 0.157 0.192 5.94 5.0 0.19 0.302 5.98 5.5 0.25 11.8 6.5 0.38 0.447 0.553 19.1 7.0 0.45
b2 - 2b1,, g/cm
eq 30” 0.76 1.80 4.31 17.8 38.1
b 47.9 38.3 30.6 25.0 16.0
eq 39” 0.13 0.76 1.71 14.7 36.2
“The constants in eq 30 were determined from the rheometeric data in Figures 6 and 7. The constants in eq 39 were determined from rod-climbing data (Joseph et al., 1984).
7
q,- 0.791 Pa-s X = 0.368
F
0.1
‘
1
1
I
100
10
I
SHEAR RATE y
0
1000
Figure 6. Viscosity of PMMA solutions (Joseph et al., 1984) along with FENE-P dumbbell curves based on eq 20.
I
UJ’ I
I
10
I
2
3 4 5 6 7 S SHEAR RATE y ( 8 - I )
9
10
Figure 8. Viscosity data for a poly(acry1amide) solution (Walters and Savins, 1965) along with FENE-P dumbbell curve. The zero shear rate viscosity vo is defined as 70 = aXHb/(b + 3).
1
1
1000
100
SHEAR RATE y (s- I
)
Figure 7. First-normal stress difference for PMMA solutions (Joseph e t al., 1984) along with FENE-P dumbbell curves based on eq 21.
free surface h as a function of the distance r away from the rod axis and the angular velocity W of the rod h(r,W) = ho(r)- (R,/2u(pg/~)”~)[(4(bz- 2 b d / ( 4 + A)} + (pR,2/(2 + A)}]W 2+ ... (39) in which u is the surface tension, A = R,(pg/u)’/2 is the dimensionless distance from the rod axis, ho(r)is the “static climb” (including meniscus effects), R, is the rod radius, and g is the gravitational acceleration. Joseph and coworkers measured h ( r = R,, w) for values of W approaching zero. The quantity (b, - 2bll) was determined from the measured value of the slope, dhldW2,at r = R, and as W 0 along with eq 39. Using the viscosity and the first-normal stress difference data on PMMA solutions of Joseph and co-workers (Figures 6 and 7 ) ,we determined the FENE-P dumbbell second-order fluid constants and calculated (b, - 2b11). The results are given in Table I. Surprisingly, the agreement between theory and experiment seems to improve for higher concentrations which is contrary to what is expected for a model for dilute solutions. Joseph and Fosdick’s (1973) and Giesekus’ (1963) analyses of the rod-climbing experiment yield an approximate expression to second order in W for the stress field from which the trace of the stress tensor can be found: trTp = -8(b11 - b2)(R,/r)4W2 (40)
-
0
(8-l)
n
Y
ANGULAR VELOCITY W
(9-’ )
Figure 9. Torque on a rotating sphere in poly(acry1amide) solution (Walters and Savins, 1965) of Figure 8 along with the FENE-P dumbbell prediction.
This result along with eq 10 and 16c shows that the dumbbells are most highly stretched at the rod surface, with the end-to-end distance decreasing to its equilibrium value as r increases. (e) The Torque on a Slowly Rotating Sphere. The flow around a rotating sphere of radius R, and constant angular velocity W, which is submerged in a vat containing a viscoelastic liquid of density p, has been studied both analytically and experimentally. Much of the analytical work has been focused on perturbation solutions using ordered fluids up to fourth order; for example, Thomas and Walters (1964) developed an expression for the torque, I , required to maintain steady rotation of a sphere in a second-order fluid. In a later paper by Walters and Waters (1964), the analysis was extended to third-order fluids, and the torque is given by
Using the viscosity data of Walters and Savins (1965) for a 1.5% aqueous solution of poly(acry1amide) (Figure 8),
Ind. Eng. Chem. Res., Vol. 27, No. 7, 1988 1319 we determined two groupings of molecular parameters, namely A’ = 0.36 s (see eq 20) and aAH[b/(b + 3)] = 0.791 P*s, taking qs to be zero. Since no \kl data were available, the values of qs, AH, a , and b could not be determined uniquely. From the groupings of constants above, we determined bl and bl:ll (bll and b12 are zero for all FENE-P dumbbells, and b2 was assumed to be negligible from information on streamlines). Then the torque was calculated from eq 41 as a function of W, the curve resulting from molecular theory is shown in Figure 9 along with experimental values. The agreement is regarded as satisfactory in view of the limited rheometric data. Finally, using the results of Fosdick and Kao (1980) for a fourth-order fluid, the trace of the stress tensor is trTp =
pR,2W2 5 - 1 ~
- (21 - 35) sin2 8 + 2(5 - 1) X _________________________
where 4 = r/R, is a dimensionless radial distance from the center of the sphere. The dotted underlined term is associated with fluid inertia, and the remaining terms arise from elasticity. By combining eq 10, 16c, and 42 and solving for ( Q 2 ) / Q 2 it , can be shown that the largest average end-to-end distance for the dumbbells takes place at 8 = 7r/2 and 5 = 1. Unfortunately, no data for molecular stretching are available for comparison. 8. Conclusions
-
The main result of this paper is the illustration of the complete sequence of activities: molecular model constitutive equation .+solution of flow problems .+ calculation of molecular stretching in the flow field. In this illustration we use the FENE-P dumbbell because it is capable of describing some rheometric properties and leads to results in complex flows which compare well with experimental data. Complete evaluations of the FENE-P dumbbell constitutive equation (or any other equation) are hampered by the lack of experimental rheometric data and experimental data from complex flow geometries, both for the same test fluid. It should be emphasized that the discussion given here is not intended to promote the FENE-P dumbbell model, since more realistic models do exist, but rather to demonstrate the connection between kinetic theory and the solving of flow problems. The FENE-P dumbbell model has been generalized to the FENE-P chaiqwith consistently averaged hydrodynamic interaction (Ottinger, 1987),and material functions for several shear flows have been calculated (Wedgewood and Ottinger, 1988).
Acknowledgment We thank the National Science Foundation (Grant CPE-8104705), the Vilas Trust Fund of the University of Wisconsin, the MacArthur Professorship, Dow Chemical, and the Universal Oil Products Fellowship for financial support. We also thank Dr. R. C. Armstrong of MIT for his helpful comments and discussions concerning the squeeze flow problem.
Nomenclature () =
phase-space average
[ [ ] ] = momentum-space average a = nkT b = dimensionless parameter ( = H Q o 2 / k T ) b,, b2, bll, etc. = retarded motion coefficients
F(C)= spring force Q = dumbbell end-to-end vector Qo = maximum extension of spring $ = configurational distribution function { = friction coefficient k = Boltzmann’s constant T = absolute temperature K = same as (Vv)* (that is, K~~ = au,/ax,) H = spring constant 8 = unit tensor cy = structure tensor ( = ( Q Q ) H / k T ) = ij component of the inverse of the cy tensor 2 = see eq 10 AH = time constant (={/4H) X = time constant (see eq 34) A’ = time constant (=(27/2)lI2[b/(b + 3)312]kH) T = stress tensor T~ = polymer contribution to the stress tensor i. = shear rate 7 = viscosity 7, = solvent viscosity \k, = first-normal stress coefficient \E2 = second-normal stress coefficient i j = elongational viscosity ij+ = elongational stress growth viscosity ijm+= value of q+ at t = m y(,),y(2j,etc. = rate of strain tensors Rd = disk radius R, = rod radius R, = sphere radius G = surface tension ho = half of initial plate separation F * = dimensionless applied force (see section 7a) t l j 2= squeezing time for h = ho/2 W = angular velocity I = torque
Literature Cited Armstrong, R. C.; Gupta, S. K.; Basaran, 0. Polym. Eng. Sci. 1980, 20, 466-472. Bird, R. B. Viscoelasticity and Rheology; Lodge, A. S., Renardy, M., Nohel, J. A,, Eds.; Academic: New York, 1985; pp 105-123. Bird, R. B.; Armstrong, R. C.; Hassager, 0. Dynamics of Polymeric Liquids, Volume 1, Fluid Mechanics; Wiley-Interscience: New York, 1987a. Bird, R. B.; Curtiss, C. F.; Armstrong, R. C.; Hassager, 0. Dynamics of Polymeric Liquids, Volume 2, Kinetic Theory; Wiley-Interscience: New York, 198713. Bird, R. B.; Dotaon, P. D.; Johnson, N. L. J. Non-Newtonian Fluid Mech. 1980, 7,213-235. [Because of an error in eq 58, pointed out by H. H. Saab, pp 224-234 are invalidated.] Burdette, S. R. “Experiments on the Flow of Viscoelastic Fluids between Eccentric Rotating Cylinders.” Ph.D. Dissertation, The University of Wisconsin a t Madison, 1987, pp 166-173. Fan, X. J. J. Non-Newtonian Fluid Mech. 1985, 17, 125-144. Fosdick, R. L.; Kao, B. G. Rheol. Acta 1980, 19, 675-697. Giesekus, H. Rheol. Acta 1963, 3, 59-71. Joseph, D. D.; Fosdick, R. L. Arch. Ration. Mech. Anal. 1973,49(5), 321-380. Joseph, D. D.; Beavers, G. S.; Cers, A.; Dewald, C.; Hoger, A.; Than, P. T. J . Rheol. 1984,28, 325-345. Kolmogoroff, A. Math. Ann. 1931, 104, 415-458. Leider, P. J.; Bird, R. B. Report 22, 1973; University of Wisconsin, Rheology Research Center, Madison. Lodge, A. S.; Wu, Y. Rheol. Acta 1971, 10, 539-553. Manke, C. W.; Williams, M. C. J. Non-Newtonian Fluid Mech. 1985, 19, 43-52. McClelland, M. A.; Finlayson, B. A. J. Non-Newtonian Fluid Mech. 1983, 13, 181-201. McClelland, M. A.; Finlayson, B. A. J. Rheol. 1988, 32, 101-133. Mochimaru, Y. J . Non-Newtonian Fluid Mech. 1981, 9, 157-178.
I n d . Eng: Chem. Res. 1988,27, 1320-1329
1320
Mochimaru, Y. J . Non-Newtonian Fluid Mech. 1983,12, 135-152. Nada, I.; Yamada, Y.; Nagasawa, M. J. Phys. Chem. 1968, 72, 2890-2898. Ottinger, H. C., private communication, 1986. Ottinger, H. C. J . Non-Newtonian Fluid Mech. 1987, 26, 201-246. Peterlin, A. J. Polym. Sci., Polym. Lett. 1966, 4B, 287-291. Tanner, R. I. Trans. SOC.Rheol. 1975, 19, 37-65. Thomas, R. H.; Walters, K. Q.J. Mech. Appl. Math. 1964,17,39-53. van Kampen, N. G. Stochastic Processes in Physics and Chemistry; North-Holland: Amsterdam, 1981; Section VIII.6. van Wiechen, P. H.; Booij, H.C. J . Eng. Math. 1971,5, 89-98.
Walters, K.; Savins, J. G. Trans. SOC.Rheol. 1965, 9, 407-416. Walters, K.; Waters, N. D. Rheol. Acta 1964, 3, 312-315. Warner, H. R., Jr. Ind. Eng. Chem. Fundam. 1972, 11, 379-387. Wedgewood, L. E.; Bird, R. B. In Integration of Fundamental Polymer Science and Technology;Kleintjens, L. A., Lemstra, P. J., Eds.; Elsevier,,Applied Science: London, 1986; pp 337-345. Wedgewood, L. E.; Ottinger, H. C. J . Non-Nontonian Fluid Mech. 1988,27, 245-264.
Received for review August 3, 1987 Accepted February 22, 1988
Computing All Real Solutions to Systems of Nonlinear Equations with a Global Fixed-point Homotopy Minoru Kuno* and J. D. Seader Department of Chemical Engineering, University of Utah, Salt Lake City, Utah 84112
A global fixed-point homotopy has not found wide application in chemical engineering for solving systems of nonlinear equations by continuation. Rather, the Newton and problem-dependent homotopies have been favored. However, it is conjectured here that the fixed-point homotopy would be expected to always place all real roots of the system on a global homotopy path because the path is forced to begin from a single starting point. If so, all roots could be computed by following the single path with a suitable continuation method. This would be particularly desirable when the number of real roots to a system cannot be predetermined and one wishes to compute all solutions. In this study, it was found that such a path does exist provided that a criterion is used for selecting a starting point which minimizes the number of real roots of the global fixed-point homotopy function a t an infinite value of the homotopy parameter. Several examples, including an adiabatic reactor with multiple solutions are presented to illustrate the application of the criterion. While methods have been devised recently to find all roots of polynomial equations, this is the first method for systematically locating all roots to general systems of nonlinear equations. The necessity of solving systems of nonlinear equations often arises in simulating and designing a chemical plant or optimizing a process. Newton's method or Powell's method are commonly applied to solve such systems. However, as described by Seader (1985), both methods have well-known disadvantages such as the requirement that the starting point (initial guess of the solution) must be in the neighborhood of a root. That is, these methods are only locally convergent. In addition, both methods are designed to locate, at best, just one root even though multiple solutions may exist. Continuation methods are rapidly being applied to numerous problems in engineering analysis as summarized in a recent review by Seydel and Hlavacek (1987). A solution method for nonlinear equations that is globally convergent from any starting point is the homotopy continuation method. This method, as described, e.g., by Garcia and Zangwill (1981)) does not solve the nonlinear equations directly, but gradually reaches a root by beginning from a starting point which satisfies a second, simpler system of equations. Both systems of equations are embedded into a so-called homotopy function. Thus, if f ( ~ ) is the system of nonlinear equations to be solved and g ( x ) is a second simpler system of the same number of equations, the homotopy function might be constructed as H(x,t) = tf(x)
+ (1- t ) g ( x ) = 0
(1)
*Author to whom a11 correspondence should be addressed.
Present address: Daitec Co., Ltd., 4-85 Chikara-machi Higashi-ku Nagoya-shi, Japan. Minoru Kuno was on leave from Daitec Co., Ltd., during the course of the study reported here. 0888-5885/88/2627-1320$01.50/0
where t is a scalar homotopy parameter which is gradually varied from 0 to 1as the path is tracked from the starting point to a solution. The homotopy function can be constructed in accordance with the characteristics of each system of nonlinear equations, but this effort can be time consuming. Alternatively, canonical homotopy functions can be used, two of which, as described by Garcia and Zangwill (1981) and Wayburn and Seader (1987), are the Newton homotopy and the fixed-point homotopy, where g ( x ) = f ( x ) - f ( x o ) and x - xo, respectively. Until now, most studies have focused on the use of homotopy continuation to obtain only one root even though other roots may exist. Studies to obtain all roots are recent and have been restricted largely to systems of polynomial equations, as discussed by Morgan (1987))where multiple starting points are used to obtain multiple solutions. Allgower and Georg (1980, 1983) discuss methods for approximating additional, but not necessarily all, roots of general systems of nonlinear equations. The aim of this paper is to present a method for determining a starting point from which one can find all the roots along one branch of the path when using a global fixed-point homotopy. The homotopy parameter is permitted to take on values greater than one, but the homotopy path remains bounded for all values of the homotopy parameter within the region where the real path connects all the zeros of f ( x ) . The real path may form an isola, which is a closed loop that does not contain the starting point. In general, the global Netwon homotopy and global fixed-point homotopy encounter difficulties when an at0 1988 American Chemical Society
