Himmelblau, D. M., Chem. Rev. 64,527 (1964). Ibrahim, S. H., Kuloor, N. R., Brit. Chem. Eng. 5, 795 (1960). Johnson, P. A., Babb, A. L., Chem. Rev. 56,387 (1956). Krieger, I. M., Mulholland, G. W., Dickey, C. S., J . Phys. Chem. 71, 1123 (1967). Lange, N. A., Ed., “Handbook of Chemistry,” 10th ed, hlcGrawHill, New York, N. Y., 1961 Malik, V. K., Hayduk, W., Can. J . Chem. Eng. 46, 462 (1965). Olander, D. R., A.I.Ch.E. J. 9, 207 (1963). Reddy, K. A,, Doraiswamy, L. K., IND.ENG.CHEM.,FUNDAM.
Reid, R. C., Sherwood, T. K., “The Properties of Gases and Liquids,” 2nd ed, McGraw-Hill, New York, N. Y., 1966. Ross, ill., Hildebrand, J. H., J . Chem. Phys. 40,2397 (1964). Tammann, G., Jessen, V., 2. Anorg. AZlg. Chem. 179, 125 (1929).
Unver, A. A., Himmelblau, D. M., J . Chem. Eng. Data 9 , 428 (1964).
Vivian, E. J., King, C. J., A.1.Ch.E. J . 10, 220 (1964). Wilke, C. R., Chang, P., A.I.Ch.E. J. 1 , 264 (1955). Wise, D. L., Houghton, G., Chem. Eng. Scz. 21, 999 (1966).
6, 77 (1967).
RECEIVED for review July 21, 1971 ACCEPTED April 10, 1972
Ree, F. H., Ree, T., Eyring, H., Ind. Eng. Chem. 50,1036 (1958). Ree, T. S., Ree, T., Eyring, H., J . Phys. Chem. 68, 3262 (1964).
Kinetic Theory and Rheology of Dilute Suspensions of Finitely Extendible Dumbbells by Harold R. Warner, Jr.l Rheology Research Center and Chemical Engineering Department, University of Wisconsin, Madison, Wis. 65706
A kinetic theory is developed for a dilute suspension of dumbbells which are joined by finitely extendible nonlinear elastic (FENE) connectors obeying the force-extension law F(c) = Ho/R( 1 - ( R / R O ) ~for ) R < Ro. Hydrodynamic. interaction is neglected or considered in “equilibrium-averaged” form. The equilibrium distribution function is calculated and steady-state shearing flow and small-amplitude oscillatory shearing flow are analyzed using this theory. The results for the various material functions computed vary between the Hookean spring dumbbell limit and the bead-string limit.
B e c a u s e of the success of the rigid dumbbell kinetic theory (Bird, et al., 1971) in describing qualitatively a wide range of shear-rate dependent viscoelastic phenomena, the development of a dumbbell model with a finitely extendible nonlinear elastic connector (hereafter referred to as FENE dumbbells) is a useful method of understanding the transition in behavior between the rigid dumbbell and the Hookean spring dumbbell. The finitely extendible connector is closer to the behavior of a macromolecule since a long chain molecule is neither rigid as is a rigid dumbbell nor infinitely extendible as is a Hookean spring dumbbell. Since little is known about the force-extension law for a macromolecule, the force-extension law chosen for the following calculations is one which expedites the mathematical manipulations required to obtain a perturbation solution. Although admittedly crude, the dumbbell model is considered since a dumbbell model seems to contain the essential physical ideas which are contained in models with a greater number of beads and connectors. However, the dumbbell model avoids the coupling of the equations which occurs in the multi-bead systems and which generally presents overwhelming difficulties in the mathematics. This paper is divided into four parts: development of the basic equations, calculation of the equilibrium distribution function, application of this theory to steady-state shearing flow, and analysis of small-amplitude oscillatory shearing flow using the FENE dumbbell model. The two flows which are studied most frequently are considered here. Other flows l Present address, Research and Development Department, Alantic-Richfield Co., Box 2819, Dallas, Tex. 75221.
can be analyzed using this theory; however, the math manipulations required to obtain even a perturbation solution to the equations for these flows would be very complicated. 1. Kinetic Theory-Basic
Equations
A. Basic Equations. The basic equations for the kinetic theory of a dilute suspension of dumbbells with an unspecified connector are well known (Bird, et al., 1971; Giesekus, 1966). Below are given the two equations, in the form derived by Bird, et al. (1971), which have to be solved to obtain the material functions for he various flows to be considered. R]$
2kT b
- --
-$
1‘ 3 R
- 2- P@)$})
r
(1)
’5-’5
Equation 1 is called the “diffusion equation” and is a combination of the equation of continuity for the distribution function $(R) and the equations of motion for the two beads of thedumbbell. The distribution function $(R) is such that $(R) d R represents the probability of finding a dumbbell within the orientation range dR. The other symbols in this equation are R, the orientation vector for a dumbbell; K, the tensor specifying the homogeneous flow field; k, the Boltzmann Ind. Eng. Chem. Fundam., Vol. 1 1 , No. 3, 1972
379
u
‘ 0
0.2
0.4 0.6 (R/Ro)
Figure 1. Modulus-extension laws:
----,inverse Langevin connector
0.8
1.0
-, FENE connector;
11. Calculation of the Equilibrium Distribution Function +e
r,
constant; T, the absolute temperature; the friction coefficient of a bead (Hydrodynamic interaction, the perturbation of the solvent flow patterns by the beads of the dumbbell, is neglected in the derivation of eq 1. For any dumbbell model the switch from negligible hydrodynamic interaction to “equilibrium-averaged” hydrodynamic interaction is accomplished by replacing { by r/(l - q) where p is a constant which can be evaluated once the equilibrium distribution function fie is known. For the model considered here the p factor is evaluated following eq 6.) ; and F(O),the tension in the connector of the dumbbell. Equation 2 is the equation for the stress tensor .e. The other symbols in this equation are .es, the solvent contribution to the stress tensor (the solvent x i ) where q s is the contribution .es is equal to - q s ( x solvent viscosity and (K K t ) is the rate-of-formation tensor); n,,, the number density of dumbbells; J (= S$dV), the normalization constant; 6, the unit tensor; ( ), the expectation value defined by (9) = $Q$dV where dV is a differential volume element; and b/bt, the convected derivative defined by the equation
+
+
~ - (-R R )- b- (RR) - { K . (RR)] bt bt
- {(RR) . K+]
(3)
Equation 2b is obtained by combining eq 2a with a manipulated form of eq 1. The assumptions made in the derivation of eq 1and 2 are discussed by Bird, et al. (1971). B. Force-Extension Law for the FENE Connector. The force-extension law chosen for the FENE connector is F(c) =
1
-HoR (R/R,J2
of Hookean spring dumbbells has been studied extensively, a limiting condition as ROtends to infinity is readily available. In Figure 1 is plotted the behavior of 1/(1 - (R/RJ2)over the range of R/Ro from 0 to 1. The dashed line is the behavior of the inverse Langevin function; this function represents the modulus-extension law for a chain of randomly oriented links (Treloar, 1950). The inverse Langevin spring connector has been used by Stevenson (1970) in his calculation of the elongational viscosity of a dumbbell suspension. Because of the mathematically complex form of this function, the inverse Langevin connector model has not been used in other calculations. (The elongational viscosity for dumbbells obeying the force-extension law given in eq 4 has been calculated (Warner, 1971). This result will not be presented in this paper since the behavior of the elongational viscosity as a function of the elongational rate is similar to that predicted for the inverse Langevin spring connector (Stevenson, 1970).)
(R < Ro)
The boundary conditions on this equation are (1) $e is finite at R = 0, and (2) the normalization of $e is taken into account in the stress tensor equation by the J in the denominator. Since in the equilibrium situation $e should be independent of the angle variables 0 and +, a solution of the form $e = ICle(R)will be postulated. Substitution of ilie(R)into eq 5 and application of the specified boundary conditions to the general solution to the differential equation results in the following form for $e
(The q referred to following eq 2 is obtained by averaging the Oseen tensor 0 (Bird and Warner, 1971; Oseen, 1910) over all orientations, assuming that all orientations are equally probable
(4)
where Hois a spring constant and Ro is the ultimate length. The choice of this relationship is quite arbitrary since innumerable equations can be proposed which predict Hookean spring behavior for small extensions (R/Ro < 0.2) and yet have some length beyond which the connector cannot be stretched. The method used to choose this form was to select three possible forms for the force-extension law and then to calculate the equilibrium distribution function +e. The mathematical form of $e is a very good indicator of whether or not nonequilibrium calculations such as steady shearing flow can be easily done. Another useful feature of eq 4 is that 8s ROtends to infinity, this form reduces to Hooke’s law: F(O) = HoR. Since the kinetic theory of dilute suspensions 380 Ind. Eng. Chcm. Fundam., Vol. 1 1 , No. 3, 1972
To calculate the equilibrium distribution function eq 4 is substituted into eq 1 with K = 0 and b$/bt = 0 (here S = sin 0) :
As b tends to infinity, eq 6a, after notation differences are accounted for, reduces to the form given in the footnote to eq 2 of Bird and Warner (1971) for the Hookean spring dumbbell. Also this limiting form is equal to the one given by Zimm (1956) in his eq 8 except for a factor of l/(ij - kl)’’’ which enters his expression because he is considering a multibead model.) In Figure 2 are plotted $e curves for various values of b. It should be noted that as b becomes large, the shapes of the equilibrium distribution function curves resemble the shape of the Hookean spring equilibrium distribution curve obtained from $e(H) = exp(HoR2/2kT). As b tends to zero, the equilibrium curve has a value of unity until R/Ro is nearly 1; this limit is the bead-string limit (no force in the connector until the “string is taut”).
where S = sin 0 , C = cos 0, s = sin +,and c = cos 4, The time the friction factor of constant X = rRo2/12kT contains one bead of the dumbbell, and Ro, the ultimate length of the dumbbell. The parameter b = HoRo2/kT is a dimensionless grouping of model parameters containing Ho, the spring constant in the Hookean limit. (A second time constant T = r/4Ho can be defined for this model. This time constant is the one usually associated with the Hookean spring dumbbell model. For the FENE dumbbells X and T are related through b as follows: A = ( b / 3 ) ~ . ) As Ro tends to infinity, eq 10 reduces to the differential equation for the distribution function of the Hookean spring dumbbell. The solution to that differential equation was first reported by Kramers (1944)
r,
Figure 2. Equilibrium distribution function FENE dumbbells from eq 6
curves for the
111. Steady-State Shearing Flow
A. Flow Description. We here consider the simplest shearing flow, namely the steady-state flow between parallel flat surfaces, one of which is moving a t a constant rate in its own plane relative to the other. Experimentally, flows of this type are achieved using either a concentric-cylinder device or a cone-and-plate apparatus such as the Weissenberg rheogoniometer. The velocity profile to be considered is: v, = KY and v u = v z = 0 where K is a constant called the shear rate. Three material functions are defined from the stresses produced by the shearing flow Ty,
(7)
= -qK
72,
-
Tug
=
-eK2
(8)
Tyy
-
721
=
-PK2
(9)
Here r Y z is the shearing stress, T,, - r y Yis the primary normal stress difference, and r y y- r z 2 is the secondary normal stress difference. q is the viscosity function and e and j3 are the normal stress functions. The material functions q, e, and 0 are experimentally measurable and the techniques employed have been discussed thoroughly elsewhere (Lodge, 1964). The methods of measuring q, the non-Newtonian viscosity function, and e, the primary normal stress function, are fairly well established and for most macromolecular fluids 7 and e have been found to be positive and are monotone decreasing even functions of K . The procedures for measuring 0,the secondary normal stress function, are less well understood, and a t present even the sign of 0 is in doubt. Generally, it is thought that the magnitude of 0 is, say, less than 20% of that of e and the best measurements to date seem to indicate that 0 is negative. For the dilute-solution concentration range very few reliable data are available, even for the viscosity function, the simplest of the three to measure. B. Shearing Flow Equations-for the FENE Dumbbell Theory. 1. Diffusion Equation. The differential equation for the distribution function $ is (Bird, et al., 1971; Warner, 1971) :
The form of GH obtained by expanding the second exponential term will be used as a limiting condition on the distribution function $ of the FENE dumbbell when the perturbation solution of eq 10 is obtained for the steady shearing flow. From eq 11 and from the form obtained for the distribution function in several other flow problems (Bird and Warner, 1971; Bird, et al., 1970; Kramers, 1944; and Stevenson, 1970), it is reasonable to postulate that $ can be separated into the product of two terms: one, the equilibrium distribution function $e, and the other a function which will be called $' and will depend on K .
iL
=
$'(R,~,+,K)
$e
(12)
At equilibrium $' will be unity. Substituting eq 12 into eq 10 and using the expression for $e given in eq 6 results in
a*'
Ro2
+~~2-
% ( R 2 %) _bRo2R _ _ a$' _ Roz - R2 bR
a
R2SdO =
sc +o R. a2$' --
( {
R2S2
X I )
(3XK) S2Sz[ R
=]
blC' - Ro2 - R2 dR
f
or A$' = (3XK)Q$"
(13)
Two types of solutions for $' will be obtained: first, a perturbation solution which is good only for very small values of X K , and second, a numerical solution which is valid for a wider range of X K values. It would be useful to have an asymptotic solution for X K + m , but to date no solution of this form has been found. 2. Expressions for the Stress Tensor Components. The stresses in the suspension are given by (Bird, et al., 1971; Warner, 1971) : Tux
- Ty2,s
=
Ind. Eng. Chem. Fundam., Vol. 1 1 , No. 3, 1972
381
=o
(16b)
-- -
Here
(&) =
S2T’SR’ 0
0
&(R,e,~)$(R,e,~,K)R2SdRded~
0
Two limiting cases, b eq 20 and 21: case 1;b
7
and
J
=
~2TIT~”’ $R2SdRd0d+
The quantity T ~is the ~ solvent , ~ contribution to the stress and is equal to - q s ~where q s is the solvent viscosity; no is the number density of dumbbells. The spherical harmonics Prim, where the argument of Pnmis understood to be C = COS e, are those defined by Hirschfelder, et al. (1964). The a and b parts of eq 14-16 are obtained from the a and b parts of eq 2, respectively. The result of eq 16b that p = 0 (known as the “Weissenberg hypothesis”) is a consequence of the expression for the stress tensor given in eq 2b of Warner (1971) and is applicable to all dumbbell suspension theories which are derived using the assumptions given by Bird et al. (1971). (Bird and Warner (1971) have shownth a t by using the complete Oseen tensor instead of the “equilibrium-averaged’’ form when including hydrodynamic interaction effects, the rigid dumbbell kinetic theory gives a nonzero value for p.) C. Perturbation Series Calculations. The form of eq 13 suggests that for small shear rates a power series solution for $’ should be possible. Therefore, we postulate $‘ = 1 f (3XK)$i’
+ (3Xt~)~$z’+ . . .
(17)
Since $ = $e a t equilibrium, eq 12 requires $’ to be unity when K = 0. Therefore, the first term of the power series is set equal to 1 to satisfy this constraint. By substituting eq 17 into eq 13 we find that $ k t can be found by solving the system of differential equations =
h$k’
where $o‘ is 1. The first few
(18)
O$k-l’
=
- q8
=
~ k T r
0 = 2 ~ k T r ~
(22)
This is the Hookean spring dumbbell limit and all terms dependent on the dimensionless shear rate T K vanish. For case 2; b 0 (the bead-string limit)
-
In this case terms dependent on K do not vanish. Hence, the viscosity and the primary normal stress function decrease as the shear rate increases. D. Numerical Calculations. I n this section numerical calculations using Galerkin’s method are described. This work is an extension to a three-variable problem of the Stewart and Splrensen (1972) computation of the steady shearing flow material functions of the rigid dumbbell. The equation for the distribution function $’ is A$’
- (3XK)Q$’ = 0
(13)
Since this equation is homogeneous, a normalization condition is necessary. The condition imposed here is that AO000, defined below, be unity. After a change of variable from R to z (= (R/Ro)Z) is made in eq 14-16 and in the A and Q operators defined in eq 13 and the h and Q operators are multiplied by z(1 - z), eq 13 is solved by expanding $‘ in spherical harmonics (Pnmsm and Pnmc,) and Jacobi polynomials (Jj(a,d,z), defined by Morse and Feshbach (1953), which are orthogonal over the range 0 to 1)
1 (R/Ro)2PzZ~z 6 1
+
+
( R / R o ) ~ P ~ ~cp22c2 ~ 3(2b 7)
[
(R/Ro)~-
21 (R/Ro)~]
(19)
After the distribution function is substituted into eq 14b and 15b and the integrations are performed, the following expressions are obtained for q and e 1
b+5
2 + 7 (2
2b
-
‘)1+ . .} .
b+9
18ndcTX2
e = @ + 5 ) ( b + 7 ) (l
382
0, can be obtained from
-
1
-2
and b
are found to be
$k‘
*1‘
2
03
03
3
+
g(xK)2[2(b
(20)
+ 9) (b + 11) -I
7
Ind. Eng. Chem. Fundam., Vol. 1 1 , No. 3, 1972
7
Here the parameters a and d in the Jacobi polynomial J j (a,d,z) are determined from the weight function ~”‘(1 z)‘-’. When the stresses are calculated using eq 14b and 15b, the weight function ~ ’ ’ ~ ( 1 z)’j2 from the integral in the numerator has been used to determine the a and d parameters. Once the $ ’ ( N ) is known in terms of Jj(a,d,z),it is a simple matter to convert the solution into a form which contains Jacobi polynomials which are orthogonal with respect to the weight function of the integral in the denominator. In eq 25 Prim terms with odd values of n or m are not included because the form of $’ obtained as a perturbation series earlier in this chapter included no such terms and because of symmetry arguments (Stewart, 1972). Finally, the COefficients Anmilof the finite series in eq 24 are fitted using Galerkin’s method (Mikhlin, 1964). into the left-hand side of eq 13 Inserting the series gives A $ ’ ( N ) - (3XK)Q$’(N)= S $ ’ ( N ) (26)
1.0 0
E
.95
6
Y
I .I
I.
AK
IO.
I
I
0.1
0.2
I
0.5
I
I
1.0
2.0
I
100.
Figure 3. Numerical results for the viscosity function (7 - v8) for the FENE dumbbells. The dashed curve is the rigid dumbbell theory prediction for (7 - 7,) (Stewart and Serensen 1972)
Figure 4. Numerical results for the primary normal stress difference function (0) for the FENE dumbells. The dashed curve is the rigid dumbbell theory prediction for (0) (Stewart and SGrensen 1972)
The residual, is to be made small; it would vanish if +'cN) were the true solution. Galerkin's principle is applied by setting S+'(N) orthogonal to each term finmi' of eq 25
Figure 5. Data of Noda, et a/. (19681, for the intrinsic viscosity of a series of monodisperse poly-a-methylstyrene samples in toluene (a, M = 694,000; a, M = 1,240,000; 8, M = 1,460,000; 8, M = 1,820,000). The solid curves are: A, rigid dumbbell theory (Bird, et a/., 1971); B, FENE dumbbell theory, b = 0.1; C, FENE dumbbell theory, b = 1.0; D, FENE dumbbell theory, b = 10.; E, FENE dumbbell theory, b = 31.6; F, Hookean spring dumbbell theory (Bird, et a/., 1971)
digit convergence is achieved. Equation 16a was used as an indicator of the convergence since the result of this integration should be zero. Over the range of X K and b values for which the curves are drawn, the value of the integral of eq 16a was 0.01 or smaller. ii comparison of the FENE dumbbell curves on the two plots with the rigid dumbbell results, the dashed curves, allows one to conclude that a t high values of X K the FENE dumbbell tends to the same power-law slope as the rigid dumbbell. The same power-law slope seems to be approached for a wide range of b values. In Figure 5 are displayed the intrinsic shear viscosity data of Koda, et al. (1968), for a series of monodisperse poly-CYmethylstyrene samples in toluene. The polymer samples ranged in molecular weight from 694,000 to 1,820,000. The dimensionless shear rate P is defined as
P
=
{ [qIotlsJf/RT} K
(28)
(i = 0, l ; j = 0,N - n - l ; m = 0 , n ; n = 0, N )
For plotting the FENE dumbbell curves the time constant X has been determined from the equation
This operation results in a homogeneous system of ZNi2 linear algebraic equations in Z N i 2 unknowns. Since the determinant of the coefficient matrix for this system of equations is found to be equal to zero, a nontrivial solution exists. Replacing the equation obtained from the integration of $ooOO in eq 27 by the condition that AooO'J= 1 removes the linear dependence from the system of equations and results in an equal number of nontrivial equations and unknowns. The solutions to these equations are the coefficients and these coefficientsare functions of X K and b. By increasing the value of i2', the function +'(x) can be determined as accurately as desired. The value of N needed to obtain convergence increases with increasing XK and decreasing b, so that automatic computation is a necessity. Solutions have been obtained on a Cnivac 1108 computer by solving the system of linear algebraic equations a t various values of XK and b. Once the +'(K) is known, the material functions 7 - qs, e, and P can be calculated using eq 14b and 15b or 14a-16a. In Figures 3 and 4 are plotted the numerical results for q - q s and 0 for various values of b over the range of XK values. The curves are plotted only over the regions where a t least two
This expression is obtained by manipulation of eq 20 in the zero-shear-rate limit. The F E S E dumbbell curves in this figure are much closer together than those in Figure 3 as a result of dividing the XK values by the factor (b 5)/3. From Figure 5 we can conclude that the F E S E dumbbells are superior to the Hookean spring dumbbells at predicting curves qualitatively similar to these experimental data. Also the FEKE dumbbell curves are qualitatively better than the rigid dumbbell curve since a range of curves is possible from the FENE dumbbell theory; however, none of the curves quantitatively fits the data. Because of the slow convergence of the numerical computations which limited the calculation of the viscosity function to a 50y0 decrease from its zero-shear-rate value, only the following results xere obtained concerning the role which the increase in the axial extension ratio plays in determining the viscosity function. While 7 - q s for the FEXE dumbbells is decreasing to half its zero-shear-rate value, the root-meansquare length (R2/Ro2)1/2 is found to increase to 1.10 to 1.15 times its equilibrium value. This tends to indicate that a t
+
Ind. Eng. Chem. Fundam., Vol. 1 1 , No. 3, 1972
383
Table 1. Percentage of FENE Dumbbells Which Are near Their Ultimate Extension length Ro during a Steady Shearing Flow’
b
0.1
Percentage of dumbbells in range 0.8 < ( R / R ~< 0.9 (R/R~ 0.9 1.o
