T H E U N I V E R S I T Y OF M I C H I G AN COLLEGE OF ENGINEERING Department of Chemical and Metallurgical Engineering Technical Report A STUDY OF THE TAYLOR-COUETTE STABILITY OF VISCOELASTIC FLUIDS Chester Miller J. D. Goddard ORA Project o06673 submitted too NATIONAL AERONAUTICS AND SPACE ADMINISTRATION GRANT NO. NsC-659 WASHINGTON, D.C. administered through~ OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR November 1967

This report was also a dissertation submitted by the first author in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The University of Michigan, 1967. ii

ACKNOWLEDGMENT This is a report on work on the rheology of non-Newtonian fluids which was supported by the National Aeronautics and Space Administration Grant NsG-659 to The University of Michigan, College of Engineering. iii

TABLE OF CONTENTS Pa.ge LIST OF TABLES vi LIST OF FIGURES vii NOMENCLAUTRE ix ABSTRACT xiii CHAPIER 1. INTRODUCT ION 1 1.1 Preliminary Comments 1 1.2 Summary of Viscoelastic Fluid Properties 3 1.3 Kinematics 6.1.4 Constitutive Equations 11 1.4.1 Definition of the Stress Tensor 11 1.4.2 Fluid Models 12 1.4.3 Definition of a. Noll Simple Fluid 18 1.4.4 Principle of Fading Memory 21 1.4.5 Viscometric Flows 23 1.5 Discussion of Previous Work 26 1.5.1 Theoretical 26 1.5.2 Experimental 38 1.6 Objectives 43 2. EXPERIMENTAL WORK 44 2.1 Samples Studied 44 2.2 Sample Preparation 45 2.3 Stability Apparatus 46 2.4 Stability Experiments 48 2.5 Viscosity Mea.surements 50 2.6 Results and Discussion 56 3. THEORY OF SMALL PERTURBATIONS ON VISCOMETRIC FLOWS 70 3.1 Equations of Motion 7.1 3.2 Derivation of the Stress Perturbation S(1) 73 5.5 Derivation of H(1) (ti ) in Closed Form 77 3.4 Viscometric Perturbations 83 3.5 Rheology of Viscometric Perturbations 86 iv

TABLE OF CONTENTS (Concluded) CHAPTER Page 4. TAYLOR-COUETTE STABILITY 90 4.1 The Primary Flow 90 4.2 The Stress Perturbation j(1) 92 4.3 The Disturbance Equations 98 4.4 Low-Shear-Rate Approximation 104 4.5 Previous Analyses 106 4.6 Special Cases 108 4.6.1 Negligible Inertia 108 4.6.2 Plane Couette Flow (Simple Shear Between Parallel Planes) 109 4.7 Solving the Disturbance Equations 110 4.7.1 First Analytic Method 110 4.7.2 Second Analytic Method 112 4.7.3 Numerical Method 114 4.8 Qualitative Discussion 118 4.9 Presentation and Discussion of Numerical Results 123 4.9.1 The Newtonian Case 124 4.9.2 Case of a Fluid Defined by S+pI = 2r (y)E 125 4.9.3 Approximate Analysis 130 4.9.4 Effect of Varying P2, P3, N4, and 85 133 5. SUMMARY AND CONCLUSIONS 140 APPENDIX A. Ht (t') IN CLOSED FORM 143 B. EXPERIMENTAL DATA 147 C. CRITICAL EIGENFUNCTIONS 154 D. COMPUTER ANALYSIS 157 BIBLIOGRAPHY 169 v

LIST OF TABLES Table Page 2.1 Systems Investigated 45 2.2 Comparison of Experimental Studies 68 4.1 Material Functions in Previous Analyses 107 4.2 Effect on Tc of a Small Change in ~i(Tc) 122 4.3 Effect on cc of a Small Change in ~i(Tc) 123 4.4 Values of Tc and Ec for Newtonian Case 125 4.5 Numerical Results for a Fluid Defined by (4.110); a = -1 127 4.6 Fractional Reductions in Tc and cc 130 4.7 Material Constants used to Approximate Huppler's24 Data 131 4.8 Comparison of Approximate Theory with Experiment 132 4.9 Effect of Varying P2 to P5 for P-75-XH CMC 133 4.10 Effect of Varying P2 to P5 for Natrosol 250-H 134 vi

LIST OF FIGURES Figure Page 1.1. Schematic diagram of simple shear flow. 4 1.2. Schematic of a viscosity curve for a viscoelastic fluid. 5 1.3. Geometrical specification of contact forces. 11 2.1. The stability apparatus. 47 2.2. Transition to a vortex pattern. 49 2.3. Schematic of plate-and-cone device. 51 2.4. Schematic of rheogoniometer. 52 2.5. Schematic of the method of temperature control. 55 2.6. Viscosity vs. shear rate for solutions of 7MT CMC. 57 2.7. Viscosity vs. shear rate for solutions of P-75-XH CMC. 58 2.8. Viscosity vs. shear rate for solutions of Natrosol 250-H. 59 2.9. Critical Taylor number as a function of concentration for solutions of 7MT CMC. 61 2.10. Critical Taylor number as a function of concentration for solutions of P-75-XH CMC. 62 2.11. Critical Taylor number as a function of concentration for solutions of Natrosol 250-H. 63 2.12. Critical wave number as a function of concentration for solutions of 7MT CMC. 64 2.13. Critical wave number as a function of concentration for solutions of P-75-XH CMC. 65 2.14. Critical wave number as a function of concentration for solutions of Natrosol 250-H. 66 vii

LIST OF FIGURES (Concluded) Figure Page 4.1. Plot of Tc against 1 for a fluid defined by (4.110); a = -1. 128 4.2. Plot of Ec against 1 for a fluid defined by (4.110); a = -1. 129 4.3. Mapping of the 52-P3 plane into the Tc-ec plane for a 0.51 solution of P-75-XH CMC; a = -1, a = 12, 4 = 2, =, = 0. 135 4.4. Mapping of the P2-P3 plane into the Tc-Ec plane for a 0.9% solution of Natrosol 250-H; a = -1, a = 12, A4 = 2, =5 = 1, l s = O. 136 4.5. Mapping of the P4'-5 plane into the Tc-EC plane for a 0.5% solution of P-75-XH CMC; a = -1, a = 12, 2 = 3 = t9 = O. 137 4.6. Mapping of the N4-P5 plane into the TC-Ec plane for a 0.9% solution of Natrosol 250-H; a = -1, a = 12, =2 = P3 = 9 = 0. 158 C.1. Critical eigenfunction *c vs. X for Newtonian case; a = -1. 155 C.2. Critical eigenfunction Vc vs. X for Newtonian case; a = -1. 156 viii

NOMENCIATURE The principal symbols are presented below. Those quantities having only temporary or limited significance are not mentioned here, but are clearly defined in the text. The dyadic notation is discussed in Section 1.3. Symbol Definition a a = rl/d CMC C arboxymethylce llulo se D Differential operator = d/dX D/Dt Material time derivative, (1.16) ~ /, t Jaumann derivative, (1.17) d d = r2 - rl E Rate of strain tensor, (1.14) El Elasticity number = Tro/pd2 HEC Hydroxyethylcellulose Ht(t') Deformation history tensor, (1.32) I Identity tensor MW Molecular weight O Zero tensor p Pressure R Reynolds number = pd2Ql/r Rm Modified Reynolds number = R/5o Rt(t') Rotation tensor, (1.33) and (1.34) ix

NOMENCLATURE (Continued) Symbol Definition r Radial coordinate r1 and r2 Radii of inner and outer cylinders respectively S Stress tensor, (1.18) T Taylor number, (1.1) for large-gap case and (1.60) for small-gap case Tc Critical Taylor number Tc,N Critical Taylor number in Newtonian case Tm Modified Taylor number = T/52 t Present time t' Time previous to present time (tt < t) u Perturbation on radial velocity component u* u* = u/2]d V Eigenfunction, (4.40) v Perturbation on G-velocity component v* v* = v/21d Iv Velocity vector w Perturbation on axial velocity component w* w* = w/ald X X = (r-rl)/d Z Z = z/d z Axial coordinate x

NOMENCLATURE (Continued) Greek Symbol Definition a a = (0i2-01)/O Po Po = To/r PBi Rheological parameters, (4.47) - (4.55) r Transpose of velocity gradient tensor Shear rate A Rheogoniometer tangential stress armature displacement CE dimensionless wave number = -d/k Cc Critical wave number Bl Viscosity TIo Zero-shear viscosity 1,Tjz2,rl3, and B14 New material functions defined in (4.30) G Angular coordinate Gc Rheogoniometer cone angle X Height of vortex cells Perturbation variable p Density al and a2 First and second normal stress differences respectively a3 and o4 New material functions defined in (4.30) T Characteristic time for fluid Eigenfunction, (4.40) xi

NOMENCLATURE (Concluded) Greek Symbol Definition 1 and SQ2 Angular velocities of inner and outer cylinders respectively Vorticity tensor Xu Angular velocity of rheogoniometer platen The superscript (o) refers always to quantities associated with the primary motion, while the superscript (1) denotes secondary flow parameters. xii

ABSTRACT Theoretical and experimental consideration is given to the stability of viscoelastic flow between concentric rotating cylinders. In the mathematical analysis, the general Coleman-and-Noll simple fluid with fading memory is employed as the rheologica.l model. Since all previous treatments of this stability problem have dealt with special cases of the simple fluid model, the earlier works are shown to represent special cases of the current investigation, General techniques are developed for treating small perturbations on "viscometric flows" of simple fluids, and these are subsequently applied in analyzing the pertinent stability behavior, under the usual assumption of neutral, axisymmetric disturbances. It is found that there are exactly eight material func.tions of the rate of shear which are necessary to define this problem. Three of these functions can be determined from existing laboratory tests, but the other five are entirely new, a.nd their form can be deduced theoretically only in the limit of vanishing shear ra.tes~ Another important result of the theoretical investigation concerns the definition of one of the stability criteria,, the Taylor number. It is shown that this quantity should be calculated using the appa.rent fluid viscosity, as evaluated at the average shea.ring conditions in the a.pparatus. Stability tests and viscosity measurements were conducted on aqueous solutions of cellulose-derivative polymers whose three determinate material functions have been reported in the literature at specific concentrations. In all instances, it is found experimentally that the flow is less stable than in the Newtonian case and that the wavelength of the disturbance is greater. By comparing these findings with the results of numerical stability calculations, one concludes that certain well-known approximate models do not suffice to describe the observed stability behavior. X1iii1

CHAPTER 1 INTRODUCTION 1.1 PRELIMINARY COMMENTS The tangential flow of any fluid between two long coaxially rotating cylinders (Couette flow) represents a curved motion which is inherently unstable due to centrifugal effects. In 1923, G. I. Taylor39 first considered this stability phenomenon for the case of incompressible Newtonian liquids, and showed both theoretically and experimentally, that at a certain critical speed of the inner cylinder, a steady, laminar disturbance appears in the flow. This disturbance takes the form of a set of cellular toroidal vortices spaced regularly along the cylindrical axis. In recent years, two dimensionless parameters, the Taylor number T and the dimensionless wave number ~, have come to be associated with the above effect. The Taylor number, which is a measure of the ratio of centrifugal forces to viscous forces, is given by T = 4( a1 + (a+l R2 where a = rl/d (1.1) d = r2 - r1 R (a Reynolds number) = pd201/ and where Q1 and Q2 are respectively the angular velocities of the inner 1

2 and outer cylinders, rl and r2 are the corresponding radii, p is the fluid density, and j is the viscosity. The dimensionless wave number is defined by C= Td/X (1.2) where X is the height of the vortex cells. Couette flow remains stable as long as the Taylor number is maintained below a critical value Tc. At T = Tc the characteristic vortex pattern develops, and the critical wave number Ec then describes the vortex spacing. In the case of viscoelastic fluids, non-Newtonian effects such as shear thinning, normal stresses, elasticity, and perhaps other phenomena which have not as yet been discovered may influence the values of Tc and cc determined for the Newtonian case. These possibilities are considered herein. The theoretical portion of the present study has been devoted to treating the Taylor-Couette stability of a Coleman-and-Noll9 simple fluid with fading memory. This is the most general mechanical constitutive equation currently available and has had remarkable success in predicting all known aspects of viscoelastic behavior. Experimental consideration has been given here to the stability of aqueous solutions of three different polymers; for the systems investigated, the critical Taylor number and dimensionless wave number are presented as functions of polymer concentration. The remainder of this first chapter will be devoted to familiarizing the reader with the properties of viscoelastic fluids, presenting the

3 preliminary theory, and discussing previous work. In Chapter 2, the experimental work of the current investigation is presented. A theory of infinitesimal secondary flows of "simple fluids" is developed in Chapter 3. This is then applied in Chapter 4 to the Taylor-Couette stability problem; here the disturbance equations are derived and solved, and the theoretical and experimental results are compared. Finally, the conclusions of this study are given in Chapter 5. 1.2 SUMMARY OF VISCOELASTIC FLUID PROPERTIES The following is a brief summary of some of the unusual effects exhibited by viscoelastic fluids. It is presented in order to acquaint the reader with the type of material being considered and also to establish some nomenclature. Stress Relaxation and Recoil.-Stress relaxation and recoil are the principal phenomena distinguishing viscoelastic fluids from viscoinelastic fluids. Stress relaxation is characterized by the gradual, as opposed to the instantaneous, decay of stress in a fluid which is held at fixed strain after experiencing a rapid deformation. Recoil, on the other hand, is the partial recovery of strain occurring in a deforming viscoelastic fluid when the driving force causing the motion is suddenly removed. Shear Thinning (Viscosity) and Normal Stress Effects.-Consider the motion of a viscoelastic fluid between two horizontal parallel plates which are everywhere separated by a distance d. Let the lower plate be held fixed while the upper one is moved at constant speed, V = yd (Fig. 1.1).

Syy 1t V =yd — Sxx -Ilrrot —— x Sxx Sxx W4 t Sx Sxyx Sxy Syy Fig. 1.1. Schematic diagram of simple shear flow. At steady state, the velocity profile will be linear and given by vy = 0 (1.3) Vx = 7y where the quantity 7 is referred to as the shear rate. Viscoelastic fluids are found to yield the following stress pattern for simple shear between parallel plates7: Sxy = yr(Y) Sxx - Szz = 1(y) (1.4) Syy - Szz = 2(7) Syz = Sxz = 0 where r is the fluid viscosity, and ao and a2 are the so-called first and second normal stress differences, respectively. In a Newtonian fluid, there are no normal stresses and, in addition, the viscosity is independent of the shear rate. (The reader should note that, with regard to the

5 subscripts on the a's, a notation opposite to that of Coleman, Markovitz, and Noll7 has been employed here; it was felt that the quantity commonly de'fined as the second normal stress difference (Syy-Szz) should be named 02 rather than al.) A typical viscosity curve for a viscoelastic fluid is shown in Fig. 1.2 (which, as a rule, involves several decades of magnitude of i and y). Power Law Log Y Fig. 1.2. Schematic of a viscosity curve for a viscoelastic fluid. As can be seen from this figure, three distinct regions are present. At very low rates of shear the viscosity approaches ~o the "zero-shear viscosity." For extremely high shear rates, the "upper limiting viscosity," I,, is approached. The intermediate region is characterized by a viscosity which decreases with shear rate, a phenomenon referred to as shear thinning. The familiar "power law" viscosity model is frequently found to apply in this region. The two normal stress differences al and 02 have been studied extensively for a great many molten polymers and polymer solutions. It has been found that ol is always positive in sign and considerably larger in

6 magnitude than 02, which usually assumes a value quite close to zero. Weissenberg,45 in 1948, set forth the so-called "Weissenberg hypothesis" which claims that a2 0 for polymer solutions in general. Today, this hypothesis is no longer accepted as fact by most rheologists, since recent experimental evidence points to the conclusion that 02 > 0. We therefore may write 1 >> 2 >. (1.5) Weissenberg Climbing Effect. —When a stirring rod, placed vertically in a beaker containing a viscoelastic fluid, is rotated about its axis, the fluid will begin to climb the rod. This is very different from what is observed for a Newtonian fluid. Swelling of a Jet. —When a viscoelastic fluid emerges from a tube, the resulting jet of liquid increases in diameter indicating, to some extent, the action of the normal stresses present. Phase Lag in Oscillatory Shear. -When a Newtonian fluid is placed between two infinite parallel plates and the upper plate is oscillated, it is found (in the case where inertial effects are negligible) that the shear stress is in phase with the rate of oscillation. However, for a viscoelastic fluid, the stress will lag behind the rate of oscillation, 1.3 KINEMATICS Differences in the response of various materials under identically imposed boundary conditions are the result of differences in the stress

7 strain laws, or constitutive equations, for the materials. In general, these laws can involve all the kinematic quantities which define the states of the same particle of material during its previous history of motion through space. Therefore, in our discussion of the preliminary theory, we shall first consider, briefly, the kinematics of fluid motion. The Gibbs dyadic notation is employed throughout the present work to represent vectors and second-order tensors. (For a good summary of this notation, the reader is referred to the appendices of either Transport Phenomena by Bird, Stewart, and Lightfoot2 or Principles and Applications of Rheology by Fredrickson.17 A more comprehensive treatment may be found in Vector Analysis with an Introduction to Tensor Analysis by Wills.47) The following conventions will hold here: (1) Vectors and second-order tensors will be denoted by lower and upper case letters, respectively, with a tilde beneath each. (2) The transpose or conjugate of a second-order tensor A will be given by At while its inverse, if it exists, will be A-1 (3) I will denote the unit tensor or idemfactor (identity tensor): A I = A (1.6) (4) 0 will denote the zero tensor: A O = 0 (1.7) (5) The matrix of the physical components of a tensor will be indicated by employing the symbol > Thus in cartesian coordinates, for

example, Axx Axy Axz A < = Ayx Ayy yz (1.8) Azx Azy Azz while, in cylindrical coordinates, Arr ArG Arz Ar AGr AGz (1.9) Azr AzG Azz Further, (6) A = At implies that A is a symmetric tensor. (7) A = -At implies that A is an antisymmetric tensor. (8) A-1 - At implies that A is orthogonal. (9) Z denotes the gradient (vector) operator. (10) tr A shall denote the trace of A a scalar: tr A = (A I). (11) A ~ A is given by A2 Now that we have established this notation, let us proceed with the discussion of kinematics. The kinematic state of a fluid at some timet (taken to be the present time) is determined if we know the velocity, acceleration, etc., of every particle of the fluid at that time. We shall assume that the fluid is s c~nrnhuum ana tratna the xifnematib variables associated with a point fixed in space may be regarded as continuous functions of the spatial coordinates.

Of primary concern here is the motion which occurs relative to a particular fluid particle since the Theological behavior of a material is presumably independent of the motion as a whole (the "principle of local actiont42 ). Let P and P' be two material points of a fluid, which at time t are located at points x and x+dx with respect to the same frame of reference. The velocity at point P relative to point P' may be written as dv = r dx (1.10) where r is the transpose of the velocity gradient tensor: X= (Jt. (1.11) In cartesian coordinates, avx vx v 6V 5v 5v r(Yx y z (1.12) avz 6vz 6vz Let us now express r as the sum of a symmetric and an antisymmetric part: r = 1 (r+rt) + 1(rr) (1.13) e\.- 2 VV 2 /-' For the symmetric part, we write _2 (r)rt) (1.14)

10 and for the antisymmetric part, - = - (rr- ). (1.15) =%.- 2 2 The tensor is a measure of the rate at which strain occurs; it is called, indeed, the rate of strain tensor. Q represents the average rate of material rotation near a fluid particle and for this reason, it is referred to as the vorticity tensor. There are two time derivatives which are also frequently employed in considering the Theological behavior of a fluid. These are defined as follows: (1) The material derivative, DA/Dt of a tensor A is the derivative of A as seen by an observer who is "translating with a given fluid particle." From the laws of partial differentiation we have DA 6A - - -+ v V A (1.16) Dt 6t t # - (2) The Jaumann derivative, rA/at of a tensor A is the derivative of A as seen by an observer who is not only translating with the fluid, but also rotating with its angular velocity. This derivative is given by iA c= tDi A A + A, s (1.17) k m-a Dt f~ u f This completes our discussion of some of the basic aspects of the kinematics of fluid flow.

11 1.4 CONSTITUTIVE EQUATIONS The mechanical constitutive equation or rheological equation of state for a material is an expression of the relation between the state of stress present within the material and the kinematics of the motion causing this stress. We shall begin our discussion of constitutive equations by first defining the stress tensor,S. which enters into the determination of the dynamic state of a body. Next, the various types of "fluid models" currently being considered as constitutive equations for viscoelastic fluids will be presented. A definition of the general simple fluid of Noll30 will then be given. This definition will be augmented by a discussion of the principle of fading memory. Finally, the special class of fluid motions called "viscometric flows" will be considered. 1.4.1 Definition of the Stress Tensor Let a differential element of surface dA be located in the interior of a deforming body of material and let the orientation of dA be specified by giving the unit normal vector n (Fig. 1.3). Furthermore, let the resultant contact force, acting on the side of dA towards which n is directed, be denoted by sdA where s is the stress vector. n sdA Fig. 1.3. Geometrical specification of contact forces.

12 In general, s will depend on n. One can show by a force balance that this dependence may be expressed in the form s = S ~ n (1.18) where S is defined as the stress tensor. Equation (1.18) expresses the fact that the stress tensor S transforms any unit vector into the stress vector on a surface normal to that unit vector. By a balance of moments, it is furthermore possible to show that the stress tensor S is symmetric. (More precisely, this condition, which is often attributed to Cauchy, defines a "non-polar" material.42) 1.4.2 Fluid Models As understood here, "fluid models" are particular, as opposed to general mathematical statements of the constitutive equations for certain idealized materials. Some models, such as that of the Newtonian fluid, are of special importance in that they can describe the rheological behavior of real materials over a wide range of flow conditions. However, most models, due to their limited nature and sometimes arbitrary method of formulation, represent nothing more than mathematical curiosities. All previous Taylor-Couette stability analyses on viscoelastic fluids have involved consideration of one type or another of special fluid model. Hence, we discuss fluid models here in order to put these analyses into the proper perspective with regard to the present more general treatment. The rheological equation of state for an incompressible Newtonian fluid is given by

13, + pI = 2iroE (1.19) where p is a pressure and qo is the viscosity, a material constant.7 It is a well-established fact that this equation gives an excellent description of the mechanics of a great many fluids; unfortunately there are also numerous examples of real fluids for which it is quite inadequate. In 1945, Reiner35 attempted to account for more complicated material effects in fluids by assuming that the stress tensor S could be expressed as a polynomial in the rate of strain tensor E. By making use of the characteristic equation for E, he was able to show that the most general equation obtainable for an incompressible fluid would then be S + pI = 2j1E + 4o2E2 (1.20) where a, and c2 are scalar functions of the invariants of E. Because of the later work of Rivlin,36 Eq. (1.20) has come to be known as the ReinerRivlin model. The Reiner-Rivlin fluid gives a definite improvement in some respects over the Newtonian model. A fluid obeying (1.20) would, for example, not only exhibit a shear-dependent viscosity but also normal-stress effects. However, the two normal-stress functions ol and a2 resulting from the Reiner-Rivlin equation are found to be equal to one another, a situation contrary to what is observed experimentally on most polymer solutions and melts. Furthermore, Eq. (1.20) is incapable of representing fluids which are elastic in nature since such materials exhibit stress relaxation; in

14 contrast, this model predicts an instantaneous decay of stress (to a hydrostatic pressure) when motion ceases (E = 0). As Fredrickson17 points out, there are no known substances, other than Newtonian fluids, which are described by the Reiner-Rivlin theory. During the past twenty years, considerable progress has been made in the development of fluid models which presumably can better represent viscoelastic fluid behavior. There currently appear to be four types which are enjoying prominence: differential models, integral models, Rivlin-Ericksen models, and "anisotropic fluid" models. The differential models are often founded on the assumption that there is some analogy between the stress-strain behavior of dilute solutions of randomly coiling macromolecules (polymer solutions) and the behavior of dilute suspensions of elastic spheres in Newtonian fluids (which have been studied theoretically). Characteristically, the differential constitutive laws are written in the form of ordinary differential equations in time which relate the stress tensor and its time derivatives to the rate of strain tensor and its time derivatives. Unfortunately, terms are all too often added in a rather arbitrary manner to these equations, in order to include certain higher order viscoelastic effects; for this reason, they cannot be expected to give more than a qualitative description of any existing fluids. One example of a differential type model is the equation

15 S + Tj1 Iw- 1i(S'E+E-S) + vl(E:S)I + C.o(S:I)Ij 2= o [E + T2 -242E + V2(E:)I] (1.21) where oy, o i, 1y e12, V1, V2, Ti, and T2 are assumed to be constants. This is the "Oldroyd31 8-constant model." The so-called Boltzmann principle of superposition (see e.g., Fredricksonl7) provides the basis for a certain class of integral models of viscoelastic fluids. In deriving these models, it is assumed that, for a particular fluid particle, a small strain at one instant of time t' results in a small stress at some later time t and that the two are related through an influence function ~(t-t'). The relation between stress and strain is then generalized to a "convolution" integral of the form t stress (t) = (t-t')d[strain (t')] (1.22) t =-? by further assuming that stresses are linearly superposible. When expressed in the properly invariant form (relating the stress tensor S to integrals of kinematic quantities), the integral models can exhibit normal stresses, shear dependent viscosity, and stress relaxation. However, there is no guarantee that any existing fluids will conform to these models. Fredrickson17 cites two examples of integral models which, in the present dyadic notation, take the forms: S + pI = 2 i (t-t')Ft(t')' E(t ) F t(t')dt (1.23) -00

16 and S + pI = 2 (t -t')Ft(t') E(t') F-t(tdt' (1.24) — 00 In these equations, Ft(t' ) represents the relative deformation gradient tensor,7 a quantity which has been shown20 to satisfy the differential equation DFt(t')..= (t')' Ft(t') (1.25) Dt' subject to the condition Ft(t) = (1.26) Equations (1.23) and (1.24) have been discussed by Walters,44 who refers to them as liquids A' and B', respectively. In 1955, Rivlin and Ericksen37 developed a new theory of nonlinear viscoelasticity. By assuming that the stress within a fluid at a given time depends on the first-order spatial gradients of the displacement, velocity, acceleration, etc., at that time and by making use of certain invariance conditions, these authors were able to derive a constitutive equation which, for incompressible fluids, is expressible in the form S + pI = f[A1,A2,...,An] (1.27) where Al = 2E (1.28)

17 An+iE An + An E (1.29) and where the tensors An are commonly referred to as the Rivlin-Ericksen tensors. As with the Reiner-Rivlin theory, the Rivlin-Ericksen approach is insufficient to account for complex "memory" effects, such as stress relaxation, in viscoelastic fluids. The Ericksenl4 anisotropic fluids represent the last type of fluid models to be mentioned here. Briefly, these have been founded on the following three assumptions: (1) There is a preferred direction, defined by a vector n, associated with each particle of a given fluid. (2) Variations in the magnitude and direction of n are governed by the fluid motion. (3) The stress at any point in the fluid is a function of the rate of strain tensor E and the preferred direction n. A typical constitutive equation for an Ericksen anisotropic fluid is given by S + pI = 24E + [141+2(E:nn)]nn + 243(E.nn+nn-E) where (1. 30) Dn ~ ~ = Q ~ n + o[E - (E:nn)I] ~ n Dt.. and where 4 o0, 41, 42, and 43 are material constants. The anisotropic fluid models might be expected to give good qlualitative descriptions of certain substances (such as suspensions of rod-like or ellipsoidal par

18 tidles in a Newtonian fluid) which retain a flow-induced orientation long after motion has ceased. However, for polymer solutions and melts, they show little advantage over any of the other types of models discussed. 1.4.3 Definition of a Noll Simple Fluid No1130 has recently developed a very general theory of fluid behavior which appears capable of describing the stress patterns in viscoelastic fluids for all motions that have been studied experimentally. The essence of this theory is embodied in the definition of an incompressible simple fluid,* a substance which satisfies the following two postulates: (1) The present stress at a material point in a simple fluid is determined, to within a hydrostatic pressure, by the history of the motion in the immediate neighborhood of that point. (2) No preferred configurations exist in a simple fluid. All the types of fluid models (except the Ericksenl4 anisotropic fluids**) discussed in the previous section are seen to satisfy these two postu~~e- Act E-xa eepeni jecial cases of simple fluids. There are numerous ways of expressing statements (1) and (2) in symbolic form; for the present purposes we choose the constitutive relation which has been suggested by Goddardl9: t S+ P [Ht(t' )l (1.31) *Henceforth, the term "simple fluid" shall imply "incompressible simple fluid." **One can show that, for the present analysis, anisotropic fluids may also be regarded as simple fluids.

19 where Ht(t') gives a measure of the history of the rate of strain. It is defined, for -oo < t' < t) by Ht(t') t(t') E(t') Rt(t') (1.32) where Rt(t'), an orthogonal tensor, satisfies the differential equation DRt(t') (t) ~ Rt(t') (1.33) Dt' subject to Rt(t) = I (1.34) Physical interpretation may be given to both Rt(t') and Ht(t'). Rt(tl) represents the average rotation suffered by material in the neighborhood of a given fluid particle during the time interval (t',t); Ht(t') is the rate of strain tensor, at time t', as seen by an observer in a coordinate frame that is both translating and rotating with a particular fluid particle. The symbol 4 in Eq. (1.31) stands for a functional, that is, an operator which maps a tensor-valued function into a tensor. It merely relates the present stress S to all values assumed by the deformation history Ht(t') as t' varies between -co and t (the present time). The form of the functional 2 determines the rheological characteristics of each particular simple fluid; thus, __ will vary from one fluid to another. Consider, as one example, the case of a Newtonian fluid for which

20 t S + pI= [Ht(t')] -= 2roHt(t')]t,t = 2o 0E'(t) (1.35) t'=-oO There is an important restriction which must be placed on Eq. (1.31) due to the so-called "principle of material objectivity."7 Briefly, this principle requires that the state of stress within a deforming material be independent of the motion of the observer. For the functional of (1.31), Goddard19 gives as the required condition, t t [Q _ Ht(t) ] = [Ht(t')] Q1 (1.36) t'=- ~ t'=-~ =00 for all constant orthogonal tensors Q and for each function Ht(t'). Before concluding the present section, it is perhaps appropriate to introduce a special property of the deformation history Ht(t') which we shall have occasion to use in the forthcoming analysis. Let us begin by taking the material derivative of Ht(t'): DH (t') Dt' Dt' Dt' t (t') ~ E(t') ~ (t')] Dt' ~ E(t') Rt(t') + Rt(t') ~ Dt) rDRt(t') +Rt(t') I E(t') ~ [ D * (1.37) Substituting (1.33) into this equation and combining terms, we then have DH t(t (t' ) D(t'?) Rtt(t') L, - Q(t, )E(t') + E(t')Q(t' ~ Rt(t') Dt' 1 Dt' (1.38) But from the definition of the Jaumann derivative (Eq. (1.17)), this expression becomes

21 DHt~t') _t Rt E(t') Dt' = Rt(t') t' ~ * Rt(t'139) Thus, the material derivative of Ht(t') is related to the Jaumann derivative of E. By taking successively higher order material derivatives of Ht(t'), it is furthermore possible to show that, in general, DnHt(t') Rt(t,) t' ((t40) Rt t(i 40) Dt 1.4.4 Principle of Fading Memory The principle of fading memory has been developed by Coleman and Noll9 in order to allow for stress relaxation and other memory effects in simple fluids. Essentially, this principle implies that, for a simple fluid with fading memory, deformations which occurred in the distant past have less effect on the present stress than deformations which occurred in the recent past. In order to put these intuitive notions of fluid behavior on firmer mathematical grounds, Coleman and Noll have made use of the "norm" (or magnitude) of a deformation history which may be defined as follows (for the present analysis): norm[Ht(t')] = IHt(t' )l = 0| h2(t-t')trHt (t')dt (1.41) Here, h(t-t') is a "weighting" or "influence function" which varies with the fluid under consideration; it is positive, continuous, and realvalued and goes to zero rapidly as t-t' becomes large. The norm is interpreted as the "distance" of a given deformation history Ht(t') from the "zero-history," Ht(t') - 0Q It possesses the im

22 portant property that the weighting function h(t-t') places greater emphasis on values of Ht(t') for small t-t' (recent past) than for large t-t' (distant past). Thus, a deformation history with a "small norm" is expected to give rise to only a small contribution to the present stress. Coleman and Noll9 have made use of the idea of a norm of a history in conjunction with the operation of Frechet differentiation of constitutive functionals: If X represents the constitutive functional for a simple fluid, it is Frechet differentiable at the function Ht(t') if there exists a functional'-b (which depends on 4 ) such that the following equation holds for all lHt(t') with finite norm: t t t [Ht (t')+AIHt(t'l )] = t, I AH()]t + - ) (t''=-0 - - [tHt t )] + [/Ht(t' )!] (1.42) where - is both linear and continuous in AHt(t') and where l|[!lHt(t' )|I] satisfies the equation AHt (t' )-O JIAHt(t')ll' That is, I has a Frechet derivative at Wt(t') if the difference between the values of. at Et(t') and E (t')+Agt(t') is given to a good approximation by a continuous, linear functional of AHt(t') whenever AHt(t') is small in norm. Let us now record, for later use, one additional condition which must be met by Ax as a consequence of Eq. (1.36). This is

t t [Q Ht(t' )Qt,,Q.AHt(t') ] = Q ~ [Ht(t'),Ht(t )] Q t =-00 t =_00 (1.44) for all orthogonal tensors i and for all functions Ht(tT) and AHt(t') which satisfy (1.42). 1o4.5 Viscometric Flows There exists a special class of fluid motions, called viscometric flows, which possess extremely simple deformation histories. Their common characteristic is that they are all kinematically equivalent to simple shear between parallel plates, except for a time dependent rigid rotation of each material particle. Viscometric flows are of interest in the present study mainly because the primary motion in the Taylor-Couette stability analysis is viscometric. Fredrickson17 for one, has cited some well-known examples of viscometric flows: (1) Simple shear between parallel plates. (2) Flow in a pipe of circular cross section. (3) Tangential flow between concentric rotating cylinders (Couette flow). (4) Flow between a plate and cone. (5) Combined tangential and axial flow (helical flow) in an annulus. It has been shown by Goddard and Miller20 that the kinematics of all viscometric flows may be described by the following set of equations: The characteristic equation for E is r3

24 E3 - 2 E (. 45) -4,~ where = 12(E:E) (1.46) is the shear rate. Furthermore, 2 +,y2E = O (1.47) W = O (1,48) ~)t and -D 0 (1.49) Dt These relations can be put in a more convenient form (involving the deformation history Z.t(t')) by employing Eq. (1.40) and recalling that the tensor Rt(t ) is orthogonal. The final result is g (t 2) = Ht(t').5) y = 2(Ht(t):Ht(t )) (1.-51) D2lt(t+ ) +2Ht(t ) = (1.52) Dt,2 D[pHt(tj )] 0 (1.53) Dt' Let us now solve Eq. (1.52), an ordinary, linear, second-order differential equation, for Ht(t'). We obtain Ht(t;) = E(t)cos y(t-t') 1 sin (t-t) (1.54) 7 ~)t

25 which holds for general viscometric flows. Coleman, Markovitz, and Noll7 have considered theoretically the mechanical behavior of simple fluids undergoing viscometric flows and have shown that, for all such motions, the state of stress at a material particle is determined by three material functions of the shear rate y. These "viscometric functions" are the viscosity r(y) and two normal stresses ol(y) and o2(y), all of which are even functions of y. With particular regard to steady laminar shear between parallel plates, this means that a simple fluid will give a stress pattern identical with that represented by Eq. (1.4) for real viscoelastic fluids. Now, for flows obeying (1.50)-(1.553), it is also possible to show that the constitutive functional ( for a simple fluid may be reduced to any one of the following equivalent forms: t s PI [t(t')] = 2l(Y)Ht(t') + +2(l+2)(t' ) t J t =t 72 t t'=t + (a2-1 ) D~t(t ) (1.55) aY2 Dt tt = 20(y)E + 72 E + (arm ) )t (1-56) 2 y2 t = ('(y)A1 + A + (2I ) A2 (1.57) 7 p, 272 where we recall that An are the Rivlin-Ericksen tensors (see Eqs. (1,28) and (>129)). As was mentioned above, Coleman, Markovitz, and Noll7 have determined the minimum number of material functions necessary for describing

26 viscometric flows of simple fluids. One of the objectives of the present study is to do the same for the Taylor-Couette stability problem. 1.5 DISCUSSION OF PREVIOUS WORK 1.5.1 Theoretical The stability of fluid flow between two long concentric rotating cylinders was first considered by Lord Rayleigh34 in 1916 for the inviscid case (S+pI = 0). He derived a simple condition for instability with respect to rotationally symmetric disturbances based on energy considerations. Rayleights criterion is as follows: if the magnitude of the circulation increases outwards, the flow is stable, whereas, if it decreases outwards, the flow is unstable. Assuming that, in real fluids, viscosity serves to maintain the steady flow but does not affect the occurrence of instability, Rayleigh's criterion leads one to conclude that, for cylinders rotating in the same direction, the flow is stable if Q2ri > Qlr2 where Q1 and 02 are respectively the angular velocities of the inner and outer cylinders, and r1 and r2 are the corresponding radii. For cylinders rotating in opposite directions the circulation decreases outwards in at least part of the field of flow and the flow is always unstable. In 1923, G. I. Taylor39 published his now-famous work on the Couette stability of incompressible Newtonian fluids. This theoretical treatment has since served as a model for many later stability analyses. Taylor began by assuming that, superimposed upon the primary Couette motion

27 v(o), there is a small secondary velocity perturbation v(l) which is a function of the radial and axial coordinates but not of annular position (axisymmetric disturbance). He employed this assumption in the NavierStokes and continuity equations, dropped those terms involving products of secondary quantities, and obtained a set of "disturbance equations" which are linear and homogeneous in the components of v(l). When combined with the homogeneous boundary conditions on v(l), these disturbance equations defined a characteristic-value problem for Ql, the angular velocity of the inner cylinder. The minimum value of 1l, over all allowable eigenvectors v(l) then gives the critical conditions for the onset of instability. In solving the disturbance equations, Taylor confined attention to the neutral mode of stability where (1) neither grows, decays, nor oscillates with time, and further assumed that v(l) is periodic in the axial direction (which amounts to resolving the disturbance into its normal spatial modes). From Taylor's calculations, it is found that the criterion for the onset of instability may be expressed in the following form: 2jt4 Tc = 2 (1.38) a[oo 0571 P+O.00056/P] where P 2+ + 0.652 (1o59) and where a and a are defined in Eqs. (1.1). Equations (1.58) and (1.59) give an excellent approximation for Tc, to terms of order 1/a and for a > -1. Taylor also determined that, for the case of cylinders rotating

28 in the same direction (a > -1), the critical wave number cc will have a value approximately equal to i. From Eq. (1.2), one therefore notes that the vortex cells will have roughly the same height as width, and will thus be nearly square in cross section. Chandrasekhar5 has recently (1954) reconsidered the Newtonian stability problem for the case where the radii of the cylinders are extremely large compared to the distance between them (1/a << 1). This represents the so-called small- approximation and gives rise to considerable simplification in the mathematical treatment of the problem. In particular, for a small gap, it is found that: (1) The primary Couette motion approximates simple shear between parallel plates with Y -- (Q2-Ql)rl/d, a constant. (2) The definition of the Taylor number reduces in form to T = -2caR2 (1.60) (3) The form of the disturbance equations and their solution are greatly simplified. Chandrasekhar's results agree quite well with those of Taylor for comparable situations. For the case of a stationary outer cylinder (a=-l), he finds that Tc = 3390

29 and = 3.12 _L.1 Interest in the Couette stability of non-Newtonian fluids apparently began in 1961 with the work of Graebel.21 This author investigated the behavior of Reiner-Rivlin fluids under circumstances where the functions al and O2 in Eq. (1.20) are constants. Here, the three viscometric functions are given by r(y) = ~l = constant J1(y) = a27 (1.6l) 2 G2(y) = a2y By employing the small-gap approximation to simplify the analysis,* Graebel found that the coefficient of cross viscosity a2 can strongly influence the results obtained for the critical Taylor number and critical wave number. He showed in fact that for a2 > O, Tc < Tc,N while Cc > Cc,N' where Tc,N and c, N represent the Newtonian values of Tc and cc. Although GraebelTs analysis has been performed for the somewhat unrealistic ReinerRivlin model, it does give some general indication of what one might expect from a real non-Newtonian fluid, that is, that the critical parameters may differ from their Newtonian values. Graebel22 has also presented a treatment of the Taylor-Couette stability problem for Bingham plastics. According to the standard model, *Most non-Newtonian stability studies, including the present, have employed the small-gap approximation, and this fact will be tacitly understood throughout the further discussion.

30 such materials exhibit no normal stresses in simple shear between parallel plates but do have a viscosity which depends on shear rate: =7 L +] (1.62) where u is a constant, with units of viscosity, and Y is the yield stress. Graebel defined the Taylor number as T = -2aa P dQj (1.63) and concluded from his solution to the disturbance equations that the flow is more stable than in the Newtonian case. One might suspect, however, that a Taylor number based on the apparent fluid viscosity Iy(7) evaluated at the gap shear rate y might give a truer measure of the ratio of centrifugal to viscous forces. Had Graebel employed this definition, his above conclusion would have been reversed. This illustrates an important problem which we shall attempt to resolve in the present analysis, namely, that of determining which viscosity, in general, is most appropriate for defining the Taylor number. Thomas and Walters40'41 have considered the stability in flow between coaxial rotating cylinders of a Walters B' liquid (Eq. (1.24)), which has for viscometric functions (Y) = No = constant ol(Y) = 2Koy2 (1o64) c2(Y) = O

31 where 00 no = ~(t)dt 0 and Ko = t~(t)dt 0 They found for the physically meaningful case of Ko > O, that Tc < Tc N while wc > ~c,N Chan Man Fong has treated this same stability problem for the related Walters A' model, in order to compare results with those of Thomas and Walters. The A' liquid may be shown to yield = =o = constant 01(7) = 0 (1.65) 2 (7) = -2Koy2 as viscometric functions. For values of Ko > O, Chan Man Fong determined that Tc > Tc N and ec < EcYN7 exactly the opposite of what Thomas and Walters obtained. Chan Man Fong summarized: "Comparing the conclusions of the present study with those of Thomas and Walters, it is seen that the stability criterion is very dependent on the particular elastico-viscous model considered." It is not surprising that the results of these two studies should differ so greatly, expecially when one considers that the A' and B' models exhibit extremely different normal-stress patterns. Indeed, the question natura:lly arises as to whether such diversity in behavior might not al

32 most entirely be attributed to differences in the viscometric functions r, al, and a2. It shall become apparent from the results of the present study that this is in fact the case with regard to liquids A' and B'. However, for simple fluids in general, there will be additional material functions which can influence flow stability. In 1964, Dattall published a theoretical analysis on the Taylor-Couette stability of the particular Rivlin-Ericksen fluid represented by S + pI = a + A + 3 + A (1.66) where al, a2, and a3 are material constants. Coleman and Noll8 have shown that this model, referred to by Truesdell42 as a "fluid of second grade," describes the limiting behavior of simple fluids with fading memory at extremely low rates of deformation (or, equivalently, of fluids with very short memories). If we choose for the moment to consider the special case of extremely slow viscometric flows, we conclude from (1.66) and (1.57) that 1 = lim r(y) = al 7+0 2 = lim 2l(Y)/Y (1.67) y+O a3 = lim (a2-al)/2y2 y+O Thus, the three viscometric functions for this model are given by =l2y) a2 (168) a2(7) = (a2+2a3)y2.

Several very important results have come from Datta's treatment. First of all, Datta noted that he could reproduce the disturbance equations derived by Graebel21 for the Reiner-Rivlin fluid and by Thomas and Walters40 for their "liquid B'" (for the case of short fluid memory) merely by suitable choice of i,,1, and a2 corresponding to the viscometric functions in these investigations; now, it is furthermore possible to show that the same could have been done for Chan Man Fong's6 treatment of the liquid A' (again for the case of short fluid memory). In fact, one can demonstrate that, even for fluids with long memories, the disturbance equations for the liquid A' and liquid B' analyses could be reproduced almost exactly from Datta's equations. These findings are highly suggestive of the possibility that the three viscometric functions Ar, cl, and a2 may play a large role in determining the stability criteria. Datta has also shown that, for a fluid obeying (1.66), the sign and magnitude of the second normal stress coefficient (c2+2a3) can drastically affect the values obtained for the critical parameters Tc and cc. In particular, a small increase in (a2+2a3) can result in a relatively large decrease in Tc along with a relatively large increase in cc. In 1966, Giesekus18 reconsidered the disturbance equations developed by Datta for the purpose of treating two special cases; these were, simple shear between parallel plates where a+oo, and the case of negligible inertia where T-+O. He showed that, in both situations, the parameter (a2+2a3) is again of considerable importance. Furthermore, as part of this same work,

Giesekus put forth an interesting suggestion which is worthy of note here. He postulated that it might be possible to approximate the stability behavior of certain viscoelastic fluids outside the range of applicability of Eq. (1.66) (i.e., at finite rates of shear) by replacing the material constants Gl, 2~, and c3 in Datta's disturbance equations by the functions which they approach at y = 0. That is 0 = rl(y) 2 = C(y)72 (1.69) a3 = [G2(7)-Cl(7)]/272 This same author stressed the fact that such a substitution would not yield a strictly valid stability criterion within the framework of the general simple-fluid theory but indicated that it was the best that one could do at that time. Based on the results of the present study, we now know that this was not the case. If Giesekus had instead substituted (1.69) into Datta's original fluid model and then rederived the disturbance equations, his results would have come considerably closer to what is actually obtained for simple fluids with fading memory. Leslie25 has considered the stability in Couette flow of the idealized Ericksen anisotropic fluid represented by Eqs. (1.30) for the case of 1l = 0 and Ipol > 1. The viscometric functions for such a material are given by Leslie as

35 r (y) = " + 13 + 12nln2 al(7) = nln2 (2n2+23s)y7 (1.70) 2(7 ) = nln2( 2nj+243 )y where 2n2 = (Po-1)/4o and 2n2 = (o+1l)/4o and where nl and n2 are respectively the x and y components of the preferred-direction vector n in simple shear between parallel plates (Fig. 1.1). The author defined the Taylor number by T = -2aa Ldiu(l+ce/2) (,7 )e and then solved the disturbance equations for the case of 4 = 12 = 43, 40 = 2, O = - 2, and a = 20. His calculations reveal that when nl > 0 and n2 > O one has Tc:L. 9000 and ec c 3.0 while for nl > O and n2 < O one has Tc L 5000 and cLc _ 5.0. It is interesting to note that these results would have been quite different had Leslie defined the Taylor number as in (1.60) and, moreover, calculated its value using the apparent fluid viscosity q(y) rather than p. Under such circumstances, he would have obtained Tc ax 2500 for nl > O and n2 > 0, and Tc Ad 1400 for nl > O and n2 < 0. These findings may be compared with the "exact" value of 2275 41 - 1 given by Thomas and Walters41 for the Newtonian case (for a = - 2) Recently, Daviesl2 has investigated the Couette stability of the idealized viscoelastic fluid defined by

S + pI - Lo[(S+PI):I]E = 2 J (t-t')t(t') E(t') t(t')t' _00 (1.72) which reduces to the Walters B' liquid when o0 0, and which has as viscometric functions (~Y) = 1o [l+k1oKoY2/To] a1(7) = 2Koy2 (1.73) 92() -= 0 His purpose was to elucidate the separate effects of (a) elasticity and normal stresses, and (b) shear-dependent viscosity, on flow stability. In order to provide a basis for comparison, he also considered the behavior of the model S + pI = 20o[1+2LoKo(E: E)/o]E (1.74) which exhibits the same viscosity variation as (1.72), but yields neither normal stresses nor elastic effects. Davies defined the Taylor number in terms of the zero-shear viscosity 0o; however, he recognized the fact that, in stability problems involving a shear-dependent viscosity, two definitions of T are actually possible, one based on a characteristic parameter with units of viscosity (e.g., qo) and the other based on the true apparent viscosity I(y). Unfortunately, the author made no attempt to resolve the question as to which definition, in general, is more meaningful. The results of Davies' analysis may be summarized as follows:

37 For Ko > O and 40 < O, (1) A fluid obeying (1.74) yields TC < Tc N and Cc < Cc N. (2) A fluid obeying (1.72) gives a value of Tc which is less than that for (1.74) and a value of cc which is greater. These findings are independent of the definition of the Taylor number employed. All of the stability analyses we have discussed up to this point have involved consideration of various and particular simple-fluid models. As such, they shall henceforth become special cases of the present more comprehensive treatment, in which the disturbance equations for a general simple fluid with fading memory are derived. It will be shown here that there are eight material functions of the shear rate which can influence flow stability, and the special forms of these functions will be listed for each of the models mentioned so far. In formulating the disturbance equations for the present study, it will be assumed, following Taylor39 and Chandrasekhar,5 that the onset of instability may be characterized by an infinitesimal, steady velocity disturbance superimposed upon the primary viscometric Couette flow. We shall thus be dealing with a total flow which is "almost viscometric" in nature, that is, viscometric except for a small perturbation. The general techniques necessary for treating such motions of simple fluids will be developed in Chapter 3. Metzner, White, and Denn9 have recently suggested that flows which are sufficiently "close" to viscometric motions might adequately be de

scribed by an equation of the form (1.57), which holds in the strictest sense only for perfectly viscometric flows. If this viscometric hypothesis were correct, then knowledge of the three viscometric functions would suffice to determine the stability behavior of a given fluid a priori; furthermore, no new information concerning the fluid could then be obtained from a stability test. It is therefore considered imperative that we determine here the validity of the "viscometric hypothesis." Beard, Davies, and Waltersl have recently presented the results of the only Taylor-Couette stability treatment that is not a special case of the present analysis. Consideration was given to the so-called overstable mode of behavior for the Walters B' liquid, where it was assumed that the disturbance, rather than being steady in time (neutral mode), fluctuates with a complex frequency. The authors found that, for certain values of the fluid parameters, the overstable mode can yield a lower critical Taylor number than the neutral mode. However, as Giesekusl8 has pointed out, this type of behavior is not necessarily characteristic of all simple fluids. The experimental results of the present study indicate that, for the polymer solutions tested, the neutral mode of behavior is preferred. 1.5.2 Experimental As part of his well-known work on the Couette stability of incompressible Newtonian flows, G.I. Taylor39 performed several experiments designed to test the validity of his theoretical results. The apparatus

39 employed consisted of an opaque inner cylinder contained within an outer cylinder of glass; the cylinders could be rotated in the same or in opposite directions. A colored dye was injected into the region near the inner cylinder in order to allow for observation of the onset of instability which was characterized by the formation of a compartmentalized cellular-vortex pattern. Taylor confined attention to the case of water, and found his experimental results to agree very closely with theory. In 1928, Lewis26 extended Taylor's findings to cover a wider range of fluids and conditions. The fluid motion in this study was followed by means of tiny suspended aluminum particles rather than by dye injection. Thus, the problem of diffusion of dye into the fluid was avoided, and each experiment could be repeated several times for the same sample. Lewis' results have, for the most part, confirmed the mathematical relations derived by Taylor. The first investigators to consider experimentally the stability of viscoelastic flow between coaxially rotating cylinders were Merrill, Mickley, and Ram28 in 1962. These authors determined the shear stressshear rate relation for several polymer solutions using a Couette viscometer with a stationary outer cylinder; high rates of shear were employed and the gap between the cylinders was quite small compared to t~.e rwd2i (.a _180o' The onset of instability was observed as an abrupt increase in the slope of the shear-stress vs. shear-rate curve. rUnfortunately, the findings of this investigation must be viewed with some skepticism since, from the authors' description of the experimental pro

40 cedure, helical flow was present within the apparatus to an unknown extent. Merrill et al., have found that (1) For the Newtonian fluids tested, the critical Taylor number was about 37% higher than predicted by theory. (2) The polymer solutions studied showed no significant deviation from Newtonian behavior, provided that the Taylor number was calculated from the "thinned," apparent viscosity of the fluid. Rubin and Elata38 have recently considered the stability in Couette flow of various dilute polymer solutions and report that, for all the systems examined, the critical Taylor number increases with concentration while the cellular spacing remains relatively constant at its Newtonian value. Here, the viscosity at zero shear rate has been used in calculating Tc. As part of this same study, the authors have compared the above experimental results with the predictions of several rheological models; their attention was focused upon the idealized fluids studied by Thomas and Walters,40,41 Datta,11 Chan Man Fong,6 and Leslie.25 In addition Rubin and Elata have performed several new calculations of their own based on Leslie's analysis. They conclude: "Ericksen's anisotropic fluid for nm > 0 and n2 > 0 might be a suitable model for the investigated fluids. Not only can this model predict an increase of Tc, without appropriate changes in cc, but both Pll-P22 (al-a2) and P22-P33 (a2) are positive, which is in accordance with most of the experimental results." These conclusions are erroneous for a number of reasons. First of all, for n1 > 0 and n2 > 0, it is seen from Eqs. (1.70) that, if y > 0, al and c2 are both positive while if y < 0, they are both negative. The

41 shear rate 7, throughout Leslie's treatment, may be shown to have been taken less than zero. Hence a2 was not positive for n1 > 0 and n2 > 0 as Rubin and Elata claim. This further reflects an error in the work of Leslie. In particular, Ericksen15 has derived the viscometric functions for the fluid model defined by (1.30), and his findings are consistent with (1.70) only if absolute value brackets are placed around the <'s in these equations. The three viscometric functions will then be even functions of the rate of shear, in accordance with the theoretical results of Coleman, Markovitz, and Noll.7 As a further source of confusion in their theoretical treatment, Rubin and Elata have, following Leslie, defined the Taylor number in terms of the characteristic parameter 4 rather than the actual viscosity (-y), It will become apparent from the results of the present study that, unless T is based on q(7y) the Newtonian and non-Newtonian results may not be properly compared. The two cases considered by Rubin and Elata are listed below along with the values for Tc, the critical Taylor number based on A, and Tc N, the Newtonian critical value. 40 = 1.5, 1 2 = p2, p3 = 0, o = -1, Ec: 3o1, Tc 4600 Tc, T 35500, Tc N = 3390 o = 1.2, 42 = 4 = -1, EC 1,3.0, Tc__ 16000, Tco _n 3700, Tc N = 3390. One notes that a somewhat less pronounced effect on Tc is obtained by correctly defining the Taylor number.

42 Giesekusl8 has conducted some "makeshift" experiments on viscoelastic fluids in Couette flow as part of his previously-mentioned theoretical work. For a 4% solution of polyisobutylene in decalin, he found that the onset of instability occurred at a Taylor number which is approximately 3 x 105 times smaller than in the Newtonian case. Furthermore, the critical wave number was about 25% larger than predicted by Taylor's analysis. Giesekus also studied a 1% solution of aluminum napthenate in decalin and obtained a critical Taylor number approximately 900 times lower than for a Newtonian fluid of the same apparent viscosity. Huppler24 has recently reported the results of measurements on the three viscometric functions ir, al, and c2 for solutions of several different polymers at various concentrations. In the present work, it was originally anticipated that these functions might play an important part in determining viscoelastic flow stability. Therefore, two of the polymers studied by Huppler have been chosen for experimental consideration herein. In particular, attention has been given to the Couette stability of highviscosity grade CMC (carboxymethylcellulose) solutions, ranging in concentration from 0% to 0.55%, and to solutions of Natrosol 250-H HEC (hydroxyethylcellulose), of concentrations between 0% and 1.0%. Furthermore, in order to compare theory with experiment, the stability of a 0.5% CMC solution and a 0.9* HEC solution have, based on Huppler's viscometric functions, been considered in detail analytically.

43 1.6 OBJECTIVES The objectives of the present study may now be summarized as follows: (1) To develop general methods for treating small perturbations on viscometric flows of simple fluids. (2) To establish, as e.g., Coleman, Markovitz, and Noll7 have done for viscometric flows, the minimum number of material functions necessary to analytically describe the Taylor-Couette stability of simple fluids with fading memory. (3) To determine the relationship, if any, between the above material functions and the three viscometric functions r, ao and a2. (4) To list the above material functions for all of the particular fluid models whose stability has been previously analyzed. (5) To establish which viscosity, in general, is most appropriate for calculating the Taylor number. (6) To perform stability experiments on polymer solutions whose viscometric functions have been measured experimentally; and (7) To study analytically the stability of viscoelastic fluids whose viscometric functions are known. (8) To determine, by a combination of theory and experiment, the validity of the'Viscometric hypothesis."

CHAPTER 2 EXPERIMENTAL WORK 2.1 SAMPLES STUDIED The experimental work of this investigation was undertaken for the purpose of determining the stability behavior of some real viscoelastic fluids in flow between coaxial rotating cylinders. Aqueous solutions of three cellulose-derivative polymers, whose characteristics are summarized in Table 2.1, were chosen for study; the viscometric functions A, oa, and 02 for the latter two of these have been reported in the literature by Huppler24 for various particular concentrations.* Since the present theoretical treatment has indicated that r, l,, and c2 should greatly influence flow stability, detailed numerical calculations have been performed on a O.5% P-75-XH CMC solution and a O.9% HEC solution based on Huppler's results. For these concentrations, the current theory and experiment may be compared and, as we shall see later in Chapter 4, one may then evaluate the validity of the "viscometric hypothesis" which postulates that, for sufficiently close approximations to viscometric motions, the flow criteria are determined entirely by I, sa, and c2. *Gratitude is owed to Dr. B. Duane Marsh and Professor R. B. Bird of the University of Wisconsin for furnishing samples of these two polymers. 44

TABLE 2.1 SYSTEMS INVESTIGATED Average Viscosity Conc. Polymer Designation Producer Grade Range CMC* 7MT Hercules not known Medium 0-2.5% CMC P-75-XH Union Carbide 3.5 x 105 High 0-0.55% HEC** Natrosol 250-H Union Carbide not known High 0-1.0% *CMC = Carboxymethylcellulose **HEC = Hydroxyethylcellulose 2.2 SAMPLE PREPARATION The volume of sample required for the stability apparatus was approximately two liters. The most concentrated polymer solutions were prepared by first weighing solid polymer on a trip balance, with a least count of 0.1 g. The polymer was then transferred to a gallon container and 2500 cc of distilled water were added; complete solution was accomplished by vigorous and frequent shaking. Following Lewis,26 aluminum particles were suspended in each fluid tested in order to facilitate observation of the flow patterns. The particles could not be mixed directly with the polymer solutions, however, since under these circumstances they tended to form "clumps" due to the high viscosity and the accompanying low wetability of the solutions. To avert such difficulties, approximately 0.5 cc bulk of particles was mixed with about one liter of distilled water, which served to wet them quite readily; 500 cc of this suspension was then decanted and combined with the 2500 cc of solution already prepared to yield an intimately dispersed

system totaling three liters. Less concentrated solutions were prepared by diluting those of higher concentration. Of some interest is the effect of the aluminum particles on fluid behavior. Now, Einstein13 has derived the equation n = N0(l+2.5i) (2.1) for the viscosity of a dilute suspension of small spheres in a Newtonian fluid. Here qo refers to the viscosity of the suspending medium while / is the volume fraction of spheres. For volume fractions on the order of magnitude used, Eq. (2.1) predicts a negligible effect on fluid viscosity. Hence, it is fair to assume that the material properties of the solutions studied were not appreciably influenced by the presence of aluminum particles. Huppler24 has observed no mechanical degradation in solutions of any of the cellulose derivatives tested. In the present work, degradation due to bacterial growth may also be assumed negligible since the stability behavior of the samples was examined immediately after their preparation. 2,3 STABILITY APPARATUS The stability apparatus used in this investigation was made available by Professor W. P. Graebel (a member of the Doctoral Committee), It consisted of an inner cylinder of aluminum, having an outside diameter of 3.52 in., contained within a transparent plexiglass outer cylinder, of inside diameter 5.98 in. (Actually, these diameters varied by approximately

+ 0.005 in. along the length of the apparatus, and the values quoted'here are those measured at the center.) The ratio of the inner cylinder radius to the width of the gap for the system was then a = rl/d = 12. In order to minimize end effects, the apparatus was originally designed so that its length, approximately 56 in., was relatively large compared to the diameter; the working height of fluid within the apparatus totalled approximately 30 in. (Fig. 2.1). Fig. 2.1. The stability apparatus. The fluid shown is a 2.5% aqueous solution of 7MT CMC. Vortex cells are visible, indicating that the inner cylinder is rotating above the critical speed.

48 nected to the cylinder shafts through rubber belts; the maximum attainable speed of rotation was approximately 400 rpm. In the present work, tests were performed with only the inner cylinder rotating; there were three reasons for this. First of all, the outer cylinder tended to "wobble" when driven, and it was expected that this might influence the steady flow. Secondly, it is frequently desirable to compare stability results at a fixed ratio of the rotational speeds which was automatically assured under these circumstances. Finally, the theoretical findings suggest that similar effects would have been experienced at other speed ratios (at least for cylinders rotating in the same direction). 2.4 STABILITY EXPERIMENTS The rotational speed of the inner cylinder was increased slowly until vortex cells, which appeared as alternate light and dark bands wrapped around the inner cylinder (Fig. 2.2), began forming in certain parts of the system (primarily at the ends). It was next observed that, at a slightly higher rate of rotation, the disturbance completely filled the apparatus. The average of these two limiting speeds was taken as the onset of instability and in no case did this differ by more than 3% from either of the limits. The rate of rotation of the inner cylinder was determined by counting the number of revolutions and dividing by the corresponding time interval. Revolutions were counted by means of a small protuberance, extending from the top edge of the cylinder, which was allowed to strike a flapper, thus emitting a clearly audible sound at each rotation. The time interval for several counts was measured using a stopwatch.

(a) (b) Fig. 2.2. Transition to a vortex pattern. Figure 2.2(a) shows the stability apparatus when the inner cylinder is rotating below the critical speed. Figure 2.2(b) shows the stability apparatus when the critical Taylor number is exceeded. The fluid is a 2.5% aqueous solution of 7MT CMC.

50 The wavelength of the disturbance could be determined quite accurately in each case by measuring the combined height of several vortex cells at the surface of the transparent outer cylinder. The results of these measurements were reproducible to within 4% of the average value. After each stability experiment, the polymer solution was drained from the system and its temperature was measured using a standard mercury thermometer. It is estimated that the temperature was known to +~1F. Fredrickson17 has presented the equation max- T - r(Q2-Q2) (2.2) 4k for estimating the temperature rise due to viscous heating in a Couette apparatus. In deriving this relation, it was assumed that the outer cylinder is maintained at the constant temperature To and that no heat transfer takes place at the inner cylinder. The quantity k in (2.2) stands for the thermal conductivity of the fluid, while Tmax refers to the maximum temperature attained in the system. Calculations reveal that, in the present work, the value of Tmax-To was never greater than 1.8~C. As Fredrickson17 has indicated, a temperature rise of 1 or 2~C is tolerable for most fluids in this apparatus. Therefore, it is reasonable to assume that viscous heating effects were unimportant herein. 2.5 VISCOSITY MEASUREMENTS In this study, fluid viscosities have been measured using a commercial version of the Weissenberg Rheogoniometer, which is essentially a plate

51 and-cone viscometer. A detailed description of this instrument may be found elsewhere,16 and hence we present only the basic features of its design and operation. The fluid to be tested is sheared between a fixed flat plate and a rotating cone as shown schematically in Fig. 2.3. This flow may be regarded as having a constant rate of shear y provided that the angle between plate and cone Gc is less than 4~ (see e.g., Fredricksonl7); in the present work, Gc was 1~37'. PLATE CONE~ SAMPLE Fig. 2.5. Schematic of plate-and-cone device. The torque M on the flat plate is related to the apparent viscosity through the equation M = D3y1(7) (2.3) 12 where D is the platen diameter. The shear rate y, for a small cone angle, is given by 7 = - /OG (2.4) where y has units of sece1, X is the rotational speed of the lower platen

52 in radians/sec, and Gc is in radians. Thus, from measurements of M and u, one can determine the shear stress-shear rate relation for a given fluid. A schematic of the entire apparatus is presented in Fig. 2.4. It is seen that the lower conical platen is driven by means of a 1-hp electric motor operating at 1800 rpm; variation of the platen speed is facilitated by a gearbox connected to the output of this motor. Gearbox settings of from 0.0 to 5.9 (in steps of 0.1) correspond to speed reduction ratios of from 1:10-0'0 (1:1) to 1:10-5.9 (1:792,000), respectively. An internal speed reduction of 1:4 is present within the rheogoniometer so that the maximum rotational velocity of the lower platen is 450 rpm; by Eq. (2.4), this corresponds to a maximum shear rate of 1670 sec-1 for the present study..-Z C LAMP TORSION BAR ARMATURE TRANSDUCER IR BEARING ROTOR THERMOCOUPLE CONE AND PLATE 1800 GEAR- I4 1 SPEED MOTO R Fig. 2.4. Schematic of rheogoniometer.

53 The torque on the upper flat plate is determined by the deflection of a torsion bar joined to that platen through a rotor. This rotor passes through an air bearing which serves to minimize frictional effects during testing. The torsion bar armature deflection is measured by means of a torsion-head transducer whose signal is displayed on a meter. For the 7.5 cm diam platens and the 1/8 in. diam torsion bar (torsion bar constant = 3.43 x 10-4 dyne-cm/0.001 in.) used in this work, the shear stress yq(y) was related to the above deflection through the equation YQ(7) = 311 A (2.5) where A is the armature displacement in thousandths of an inch and y)q(y) is in dynes/cm2. In practice, the apex of the lower conical platen is truncated slightly so that the two platens do not actually touch one another; this truncation may be shown to give a negligible contribution towards the total-torque measurement. The upper flat plate is positioned vertically by means of a gap-setting transducer (not shown in Fig. 2.4) which places the hypothetical tip of the cone at the surface of the flat plate. In the present work, viscosity measurements have been made immediately after each stability experiment by using sample taken directly from the stability apparatus. A wide range of shear rates, which included that occurring at instability, were employed and care was taken to assure that, near the critical shear rate, the temperature differed by less than 0.20C from its recorded value in the stability test. The method of temperature

54 control, illustrated in Fig. 2.5, was as follows. If the initial platen temperatures (AB), determined by means of a thermocouple embedded in the upper platen, were below that desired, the platens were doused with hot tap water until their temperatures rose a few degrees above this (BC). They were then wiped dry and sample was loaded into the region between them (CD). The system was then allowed to cool slowly (for several minutes) until the temperature dropped to about 0.5~C higher than at instability (DE); at this point, testing was begun. During the entire course of measurements, the temperature was observed to fall an additional 10C (EF). Calculations have shown that the several minutes necessary for the platens to cool was sufficient for the system to reach thermal equilibrium. Thus, it is reasonable to assume that the thermocouple gave an accurate representation of the sample temperature. It was initially observed that air currents generated by the air bearing caused the platens to cool too rapidly during testing; furthermore, these currents gave rise to significant evaporation of the sample from the gap between plate and cone. To protect the system from such effects, the platens were enclosed in a wooden test cell with a plexiglass window. (This apparatus was designed and constructed by Dr. L. F. Carter4 as part of his doctoral work, and is described in detail in his Ph.D. thesis.) Fredrickson17 has presented an equation for calculating the effects of viscous heating in a plate-and-cone device. This is expressible in

C w~~~~ 0.5~C L, 0.50C F LLJ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~J w aw I — - - - Temperature in Stability Experiment AB TIME Fig. -25. Schematic of the method of temperature control.

56 the form Tmax - To - (2.6) 2k where k is the thermal conductivity of the fluid and Tmax is the maximum temperature attained in the apparatus. In deriving (2.6), it has been ass2med that both platens are kept at the constant temperature To and that no heat transfer takes place at their outer edge. For the present work, the maximum temperature rise due to viscous heating may be shown to have been less than 1.5~C in all cases. Therefore, such effects were not considered significant herein. 2.6 RESULTS A-ND DISCUSSION Flgures 2. 6-2. 6 prce.senut he-.rets-z2 as of the viscos y m2Teasurem2ents obtained in this study; the viscosity r has been plotted vs. the shear rate 7 at fixed values of temperature and composition. It is seen rrom these graphs that the aqueous solutions of all three polymers tested exhibited shear thinning. The dashed lines in Figs. 2.7 and 2.8 represent equations which were used in the present theoretical work to approximate the experimental find~g or 7ff2<zc24 a- zts l o& 2 taitDms and temperatures shown. The lines are seen to fall very roughly into appropriate positions with respect to the current experimental results; that is, those solutions which flanked Huppler's in concentration also flanked his in viscosity.

l. 250/o,25.6 oc 0 1.67%,2 6.7oC ~cn 1.o%,2 44. 5 c W2L3zOrj r0-0b- 0I 0 _ 0, 0.5 %125.6 C 4a. 4 ICTE-2tHFAR RATE (SEC') Fig. 2.6. Viscosity vs. shear rate for so1lutions of 714 CMC.

i0o 103 - 0 i 0,n,0 0.55%,24 50C Cf) I 0o.500/o,250 0 -',o,.,>_~_'~ 0.367%,24.5~C ~ t: O. 183%,24.5~C a. - I I -1010 102 103 104 SHEAR RATE (SEC-1) Fig. 2.7. Viscosity vs. shear rate for solutions of P-75-XH CMC.

59 10I w 100 N~ 1.0%,27.80c -- N 0.9/%,25 0C CJ) 0.8%,26.70C 0 10 0.48%,24.50C z I01 102 1o3 104 SHEAR RATE (SEC-1) Fig. 2.8. Viscosity vs. shear rate for solutions of Natrosol 250-H.

60 For each of the investigated fluids, the critical Taylor number as a function of concentration is plotted in Figs. 2.9-2.11; here, Tc is defined by Eq. (1.1), and the apparent viscosity at the critical shear rate r(y) has been used in its calculation. The effect of polymer concentration on the critical wave number cc is shown in Figs. 2.12-2.14. One can employ Eqs. (1.58) and (1.59) to calculate the theoretical Newtonian value for Tc for the present stability apparatus. It is found that, for a = -1 and a =12, TcN = 3620. We also recall that Taylor39 has shown that CcN N- - = 3.14. These values are indicated as broken lines in Figs. 2.9-2.14. The only Newtonian fluid considered in the current investigation was pure water which corresponds to the point at zero concentration in the above-mentioned plots. For this case, the critical Taylor number was found to be 41i40 which is about 14% higher than predicted by theory. However, this value compares quite favorably with the value, 37% higher than 28 theory, obtained by Merrill et al., for the comparable situation, and represents fair agreement between theory and experiment.* It is worthwhile noting that the point at zero concentration appears to lie on a smooth curve with the other experimental points and should thus provide a valid basis of comparison for the non-Newtonian results. The critical *One might argue that the present discrepancy between theory and experiment could be attributed to the fact that, in practice, only finite secondary motions of fluids can be observed whereas, according to the theory, the secondary flow is assumed to be infinitesimal in magnitude. Although such may indeed be the case it should be pointed out that Taylore9 and Lewis26 both achieved good confirmation of Taylor's mathematical findings using experimental procedures similar to those employed here.

L3", 4000~~- NOW low4000 w 4~~~~~~~~ H 2000 1000..No 0- - TheoreticG' NewtOn'0\ VOIue 0 1.0 ~~~~~~~2.0 0 PERCENT CMC Fig. 2.9. Critical TaYlOr inumber asa function Of coiceitration f'or solutions Of

5000 Co 4000 z a: o 3000 IJ 2000 F - -.- - Theoretical Newtonian Value o 1000 0 0.1 0.2 0.3 0.4 0.5 0.6 PERCENT CMC Fig. 2.10. Critical Taylor number as a function of concentration for solutions of P-75-XH CMC.

63 5000 4000 - 2000 -~-.- -. Theoretical Newtonian Value -I 1000 2000 0...4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 1.2 PERCENT NATROSOL 250-H Fig. 2.11. Critical Taylor number as a function of concentration for solutions of Natrosol 250-H.

64 4.0 U 3a, w w2.0: -- -- -- Theoretical Newtonian Value cr) o I I 1.0 0 1.0 2.0 3.0 PERCENT CMC Fig. 2.12. Critical wave number as a function of concentration for solutions of 7MT CMC.

65 4.0 3.0 U w w 2.0 -- --- -- Theoretical Newtonian Value - 1.0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 PERCENT CMC Fig, 2,13, Critical wave number as a function of concentration for solutions of P-75-XH CMC.

66 4.0 3.0 co W z w 2.0 -| -- - Theoretical Newtonian Value a:: 1.0 I! l I 0 0.2 0.4 0.6 0.8 1.0 1.2 PERCENT NATROSOL 250-H Fig. 2.14. Critical wave number as a function of concentration for solutions of Natrosol 250-H.

67 wave number obtained herein for the case of pure water was cc = 3.14. As can be seen from Figs. 2.9-2.14, the critical Taylor number and critical wave number both decreased with concentration for all of the systems investigated. Thus, the flow was less stable than in the Newtonian case, and the length of the vortex cells was larger. A comparison of these findings with the results of several other stability studies on viscoelastic polymer solutions is presented in Table 2.2. One notes immediately that the only study to report a decrease in the critical wave number is the current one; Rubin and Elata38 found no change in this quantity for the systems they studied, while Giesekus18 obtained increases in cc. The critical Taylor number was observed to decrease in our study, giving qualitative agreement with the findings of Giesekus; however, the magnitude of this decrease was much smaller in the present instance. By way of contrast. Rubin and Elata report increases in Tc, while Merrill et al.,28 have found, in most cases, no change in this quantity. Solutions of several different polymers have been used in all of the above investigations; thus, it appears that the stability criteria are highly dependent on the systems studied. In the theoretical portion of the present work, the small-gap approximation has been employed in order to facilitate the mathematical treatment of the Taylor-Couette problem. Unfortunately, the current experimental stability apparatus was not a small-gap system in the strictest sense of the term. This is illustrated by the fact that the theoretical Newtonian value for Tc (corresponding to the present system) was 3620 as

68 TABLE 2.2 COMPARISON OF EXPERIMENTAL STUDIES Investigator TcTc9N Cc-cN Rubin and Elata38 Greater than zero Equal to zero Giesekusl8 Much less than zero (almost equal to -TcN) 28 Zero in most cases Merrill, et al. Not determined -- * Less than zero in some cases This study Less than zero Less than zero compared with 3390 for the small-gap case.5 A comparison of theory and experiment is still possible however, if one assumes that the fractional reduction in Tc is approximately the same in both the small-gap and largegap cases. Under such circumstances, one would have Tc/Tc N a 0.86 for a 0.5% solution of P-75-XH CMC and Tc/Tc N,O.69 for a 0.9% solution of Natrosol 250-H; the corresponding values for the critical Taylor number for the small-gap case would then'be (roughly) 2900 and 2350, respectively. Before concluding the present discussion of experimental results, one last factor should receive some consideration. In the current theoretical analysis, it has been assumed that the neutral mode of behavior is exhibited by the fluid of interest at the onset of instability; i.e., the disturbance takes the form of a set of cellular vortices, steady in time. Now, as we have noted earlier, the experimentally observed vortex pattern (obtained using aluminum particles) was characterized by the appearance of a steady array of alternate light and dark bands wrapped around the

69 inner cylinder. Thus, the assumption of neutral stability appears to *be valid for the systems studied here.

CHAPTER 3 THEORY OF SMALL PERTURBATIONS ON VISCOMETRIC FLOWS As was discussed in Section 1.4.5, Coleman, Markovitz, and Noll,7 have made use of the general constitutive theory for simple fluids to analyze the mechanical behavior of these materials in viscometric flows. It has been shown that, even at very high rates of shear, the state of stress at a material point is determined by three viscometric functions l, cl, and c2. Unfortunately, most other motions of simple fluids are too complex to be dealt with by general theory, save for those cases where the deformation histories are extremely slow. As a result, rapid flows have, up to the present time, gone largely untreated. In this chapter, we shall consider a class of motions which are not only amenable to general analysis'but which also possess finite rates of strain. These are defined by the presence of a small velocity perturbation superimposed upon some primary viscometric flow. The techniques to be developed here will have obvious application to certain linear stability analyses, such as the current Taylor-Couette problem. However, it is fair to say that there exist many other instances of practical importance for which the methods discussed may find use. Let us briefly cite, as some examples, the following: (1) Combined steady and oscillatory shear between parallel plates where the amplitude of oscillation is small. 70

71 (2) Problems in which the primary viscometric stress distribution is not strictly compatable with the equations of motion so that secondary flows arise. Among these are: (a) Axial motion in a channel of noncircular cross section. (b) Flow between a plate and cone. (c) Torsional flow between rotating flat disks. It should be noted that Pipkin and Owen33 have recently derived constitutive relations appropriate for the description of small perturbations on viscometric flows. The present independent development is more general, perhaps simpler to follow, and may be applied more readily than theirs. 3.1 EQUATIONS OF MOTION The equations of momentum and continuity, governing the flow of any incompressible material, are given respectively by Dt + vVv = pb + VFS (3.1) and V-v = O (3.2) where b is the body force/unit mass (assumed here to be due only to gravity) and p is the density of the material. Let us suppose that, within a particular simple fluid, there exists a primary velocity field v(o) and an associated stress field S(o) which satisfy (3.1) and (3.2) exactly* and let these fields be disturbed slightly *The discussion of this section may be easily modified to deal with problems of the type (2) listed above.

72 by a small perturbation on the motion. Our objective in this section shall be to derive the linearized equations which must be satisfied in order for the new total flow to also be consistent with (3,1) and (3.2). For this purpose, we write v = v(~) + Sv(1) (3-3) where i has a value between 0 and 1 and where v(l) is some arbitrary, spatially smooth velocity field having the same order of magnitude as v(o) The new total-stress field S may then be assumed expressible in the form = S(e)+ tS(l) + ((3.4) where ~( ) is a tensor-valued function with the property that lim =(n) 0. (3.5) Substituting (3.4) and (3.5) into (3.1) and (3.2), we obtain oD() v( I) + V [ D( V(o + v2 D(1) v(l) P Dt Dt Dt Dt pb + V. s() + S(1) + (i (3.6) and V.v(~) + tV.v(l) = O (3.7) with D(~) + v(O).V (3.8) Dt =.

73 and D(1) (1) V (3*9) Dt where, in accordance with Eq. (3.8) and the discussion of the material derivative presented in Section 1.3, it is to be noted that the operator D(o)/Dt defines time-differentiation along a pathline of the primary motion v(~) Recalling now that the primary flow satisfies (3.1) and (3.2) identically, one has D() () + D(1) (o) + D(1) V(l V. [S()+ )(])/~] Dt Dt Dt (51. (3.10) and v.(l) = o. (3.11) Equation (3.11) represents the required perturbation on the equation of continuity. Assuming that D(l)v(l)/Dt is bounded at all points in the field of flow, we obtain our perturbation on the equation of momentum by taking the limit of (3.10) as t+0. We obtain thus r(o) (1) D(l) (o)] V S(1) P Dt..v + v = (3.12) L Dt Dt 3.2 DERIVATION OF THE STRESS PERTURBATION S(l) An expression for the stress perturbation S(l)is obtained by consideration of the constitutive equation for a simple fluid:

74 t S + pI = [Ht(t')]. (1.31) t I — 00 We first assume that the deformation history H (t') is analytic in its time arguments and then expand this quantity in a Taylor series about the present time t: 00 Ht(t') = (t'-t)n ~DnHt(t' 1) n=O n' L'Dt _n n-=O n' Dtn=> tT=t (3.13) Combining (3.13) with (1.40), we then have 00oo Ht(t') = (t'-t)n nE (3.14) n. otn n=O For the velocity field represented by (3.3), Eq. (3.14) gives Ht(t' ) = H(o)(t') + H(l)(t') + ~2H(2)(t') + (3.15) with H(~)>(t')= ) (tt) _ E(O) (3.16) It n.' 13 tn n=O _ E(1) + (t'-t)n (1) n (o)m Ht( 1)(t')- = (t'-t n (tn m~ ~~~~n=O ~n=l m=O (x) (1) (o)n-m-l ) (.17) Ot t n-m-l and so forth, where the H(i) are independent of 5 and where, for any tensor field A, (o) D() A - (o). A + A ~ (o) (3.18) 8 t A A - N r

75 and A = D A - (). A + A +..() (3.19) pt " Dt e' The quantities E(~) and Q(o) refer, of course, respectively,to the rateof-strain and vorticity tensors associated with the primary flow v(~) while E(1) and (1l) are obtained from v(l). It is to be noted that H() (t') is the deformation history due to v(~) alone. Let us now make the assumption of fading memory, so that the functional t may be Frechet differentiated at the primary deformation history H(o),H-t (t'): t t + pi= [H o)(t')] + / [Ht (t'),Ht(t')-H (t')] 3+' p = r:(t')] t'=-oo t'=0oo + 8(l t(t'1)- t,) i(t ) (3.20) where, as before, 65 is a linear and continuous functional in [Ht(t ) (o) -Hit (t')] and where iim ( (t) )] Ytit(t') -t ( ~)(t) it(ti )(ot') ||t t'-~o (' 1(()H"t ( t)- H(~ ) (t'\ ) t The norm of Ht(t )-H(~) (t) is given by t /2 whic')h,~' = Lttr[Ht(t o)-H (t)] h tt)dtJ (3.22)

76 |H (t')-H(o) t (1) W (t) (2) (t (t)) + t(H)( t') )(t,) +H (t )t (t )) +...]h2(tt')dt' /2 LH 1(t 12+2 tr[ tl)(t' ).H 2)(t')] t t t — 00 (x) h2(t-t')dt i'+... (3.23) If it is now assumed that the weighting function h(t-t1) and the functions'=(n)(t') are such that all the integrals appearing in (3.23) are bounded, we have lim t(t')-)(t') IH(l)(t') (3.24) 4 +0 z Let us next substitute Eq. (3.15) into (3.20) and recall that ~ is linear in Ht(t' )-t (t'). One may then write pi = t t =-_ t =-W (3.25) + X n 5s [Ht~)(t'),Hn)(t')]+ (I]Ht (t')-H()(t') ~ t -— 00 n=2 But, from (3.24), we have that (for n > 2) 3n [Ht (t')Htn)(t')] lim t =ItHt (t') ~(ll+0 |Ht (t')-t~)(t' || [n-l (o) rn-1 ~ [t ~)(t'),Hn)(t')] - im t- t,oo (3.26) ~-~0 llHT(1-)II~H (t, l)l

77 Therefore, the summation term on the right-hand side of Eq. (3.25) may be grouped together with the term ((lHt(t')-H(o)(t')lj to yield t t S + pI = [H~)(tl)] + e o. ["~) (t')Ht)(t'))-H~)(t')!) t'=-z1 tv=-o (3.27) where Ei: also satisfies an equation of the form (3.21). Comparing (3.27) with (3.4), we see finally that S(1) is expressible as t ~S (1) + pH I = F ))H( 1 )(t ) ] (3.28) t _00 where p\l) is a pressure. Owing to Eq. (1.44), the relation 50 [Q.. ~ ) (t).QQ. 1 ) (t'). Q ] = Q. ( [H (t')JTH t t'Q=_00 (3.29) must also hold for all orthogonal tensors Q and for all functions Htl)(t') (o) and Hit (t') which are consistent with (3.28). 3.3 DERIVATION OF Ht(1 (t') IN CLOSED FORM We now demonstrate that, for all viscometric primary flows, Eq.(3.17), our defining relation for Ht (t'), may be written in a convenient closed form. Let us begin by defining a "primary rotation tensor" Rt)(tl) which satisfies the equation DR.()= (~) * (t') ~3.30) Dt Y

78 subject to R(O)(t) = I. (.1) As a consequence of (3.30) and (3.31), it is easily verified that, for any tensor field A(t), the following relation holds: (3.32) We now consider the first (single) summation in (3.17). By employing (3531) and (3532), we have 000 co (t? t)n ")() (1) (t,_t)n rD(~)n t n! Xtn -A n. L Dtln n=O n=O Dn (t')' E(l)(t') ~)(t')] (3.33) But this expression represents the Taylor series expansion for E(l)(t' R(o )(t') about t'=t. Therefore, for any primary flow, one may write H(l)(t) = R~) (t )-E(l)(t ). R O)(t ) 02 n-' (t,_t)n (o) M ) ()n-m-L n t n.n Ctm et Dtn-m-l n=0 m=O where the quantity Rto) (t'))E( )(t')-Rt~)t') must be evaluated at positions on a pathline of the primary motion since, as has been pointed out earlier in Section 3.1, the operator D(O)/Dt' (appearing in Eqs. (3.3O)

79 and (3.33)) defines differentiation along a primary line of flow only. Now, it has been shown in Section 1.4.5 that the two relations O&(o)2 E(o) +7(0)2 E(o) = (335) and D(o) Dt 7() = (3.36) apply to all viscometric primary motions; as will soon be seen, Eq. (3.34) can be simplified further by making use of these. Defining the tensor quantities K -(l) E(o)(t) (3.37) L (1 oE( (t) (3.38) Dt:)Nm (~( = (3.40) 9 t L Dt and making use of Eqs. (3.35) and (3.36), one finds that the sums involving m in Eq. (3.34) may be written as n=l: K n=2: ) K + L -t n= F (~) K + L (o) -7 j~ ~~~ ~t

8o _ (0)2 K _ Y L M & N __ _ (o )2 M- ~ ) n=6: - 7(o)t -o 2 (0 )2 K +t L N) +=' -P7(o L -t6 Ot4 t'Vt + l _ Y(O (0)4 4LL + 0)3 2 L ()7 _2(0 )2 M _ _ _ - (o) _ and so forth. We turn attention here to only those terms in (3.1345) which involve the tensor K. The sum of all such terms Kt(tT), say, can readily be shown to be represented by ~(tct), = Ktt -y 70) t2 2 ti + 7 ()2 )4 + _ t(t4 )dt4dt3dt2dtl+. (3.4 ) where Kjt(to) iS defined as Kt(t' )= (t't)n+l,(o)n o (K.4) n2O (n+l)4 ttn

Let us now make the substitution s = t'-t into (3.41) and (3.42). We obtain for them, respectively, (Kt(t+s) = Kt(t+s) - 7 t(t+s2)dS2ds( O O + 7(o) f ( t+s4)ds4dS3dS2dsl+... o o o o (3.43) and 00 ~t(t+s) S (O)n K. (3.44)'(t+s) = X (n+l)' gtn n=O Next, taking the Laplace transform of 1t(t+s) with respect to the new variable s, we have (o)2 + (O)4 ] (t s [-+(t+s)] Ll q2 q4 q2 gi [~Kt(t+s)] (3.45) q+ ()2 where q stands for Laplace's operator. The derivative property of the Laplace transform then yields y:d (t~s~ q t(t+s)1 = qt (t+s)]. (3.46) But since, from (3.44), Kt(t+s)s= = O., Eq. (3.45) becomes [1t(t+s)] = 42 ds L (t+si (3.47) This expression may now be inverted to yield a convolution integral:

82 Kt(t+s) d= n 5L Kt(t+sl) cos y(~)(s-sl)ds1 (3.48) or equivalently 4Kt(t,) =- [ (~)(t'C-tl)dt, (3.49) where n-KtTt =0y (t'-t)" 8' K (3.50) dt' K t n' )tn n=0 Now, once again employing the "primary rotation tensor" Rt~)(t') (as has been done in Eq. (3.33)), we obtain from (3.50), Ldt'ttJ t(1) jt t o )( t ) Finally, then, R( t' = E (t R~) R (tl)cos Y(O)(t'-tl)dtl. (3.52) The same type of analysis as presented above may also be applied to the series of the terms in (3.34) involving the tensors L M and N. One ultimately arrives at an expression for Ht )(t') which holds in all cases where the primary motion v(o) is viscometric. The details of these derivations are left for Appendix A, but the final result is

H( (t R )(t') 1)(t' ) ~'(t'1) t (t') - E- (Itz'RR(~) (t')cos y(~)(t'-tl)dt1 t E(f l 0) +, ( ~Rt (tl)L- E(o)(ttj R(~)(tl)sin 7(~)'Y t ~o Pti>t(o) (tl-t')dt, +t H(~)(tl) D(lj() [(tl-t')sin 7(O)(tl-t')]dt t' DtD tD () H(o)(t 2 D )(~ ] [sin y(O)(tl-t,) y(O)(t-t')cos y(o)(tl-t')]dtl (3.53) which, owing to Eq. (1.54), may also be written in the more compact form l)t (tR) R(~)(t )' ER, (t' -' 1) t R(0)(t''R 0) (t)dt..4 (3-54) 3.4 VISCOMETRIC PERTURBATIONS In the remaining two sections of this chapter, we shall consider a class of motions which will be called "viscometric perturbations." These are flows for which the primary velocity field v(o) is viscometric and, in addition, the new total motion v(o)+ v(1) may also be regarded as such. The present section will be devoted to investigation of the kinematics of these motions, while in Section 3.5 we shall examine their rheology.

84 The following set of equations has been listed in Section 1.4.5 as describing the kinematics of all viscometric flows: Ht(t') = t Ht(t') (1. 50) 4 0 (2) D(t')2 t (t ) = (1.52) Dt'2't D [H t(t')]2 = 0 (1.53) Dt' and Dt. (1.49) Substituting (3.3) and (3.15) into these, we have [H~ ) (t')]-3 + (o) (to') (l+ ) (t 1)+(t')2H () (t') 4 ~t( (t') + H) (t' (0)2 (~)( ) r~y'l)y(~)S~)(t* +7(~) j0 ) ( t' ) (3.55) D(~) _H()(t') + 7(0)2 )(t')+i D(1) D(~) D(~) D(1l () DtV2 +t -t+ Dt' Dt' DU + D( 2 Hl)(t' ) + 24y(l)y(o?)(tj) + =(O)2 0)(t 7) f... = 0A (3.56) - [H t)(t)] + [o)(t) t D(1) (tH)(t') 12 D o) (. ) (t') )(t')(3.5)

85 and D(o) (o) + D(l) (o) D(~) 11)]+. (3.58) D~ Dt YDtJ where (l) 2 tr[E(O).E(l)] 757 -i _(0 0 -I and where, of course, the symbol... indicates terms of order t2 and higher. Recalling now that the primary flow v(o) is viscometric, we may drop those terms involving only primary quantities (those with no coefficient involving f) from (3.55)-(3.58). By dividing the resulting expressions by t and taking the limit as i approaches zero, one then obtains [H(o)(t')]H(1)(t')+H( )(t' )H(l)(t').H )(t')+H(l)(t') [H )(t')] [2y(l)y(O)HO)(t )+y(O)2Hl)(t)l) (t')] (3.60).........(1) D(O) (1) ( 0).+ H (t') + 27 (t') ~+of the form v D(~ D()v(~)ft r) which kDt' Dt Dt' Dt' oH (ti)] [H(o)(t,)]H _ D )(t')'t(t)(t)H (t)] 0 Dt' Dt' I (3.62) D(1 )Z(o) D(o) (1) + = 0 (3.63) Dt Dt

86 (1) H(~)(t') satisfies (1.49)-(1.53), (2) +0, and (3) Eqs. (3.60)-(3.63) are satisfied by H(l)(t'). 3.5 RHEOLOGY OF VISCOMETRIC PERTURBATIONS We have proven in Section 3.2 that the stress perturbation i, due to a small perturbation on any primary flow, will be given by 5(1) +( (o) (1) s) + p(I [Ht(t'),H (t')] (3.28) t'= -O where 6g is a linear functional in Hl)(t'). Let us now investigate the particularly simple form which this functional assumes whenever the velocity field v(~)+0v(l) corresponds to a viscometric perturbation. For all viscometric flows, it was shown in Section 1.4.5 that the stress tensor S may be written as t S + pI [(t')] = 2nHt(t')]t',t +2(a,+a2) Ht(t 2 t t'=t t'= - 72 + 2 LDt Htt(t')' (1.55) I2 t - t t'=t After substituting Eq. (3.15) into this expression, subtracting those portions of the resulting equation due to v(O) alone, dividing by l, and taking the limit as t+0, one finally obtains

t 87 y(D) 7 [H W)(t')7H W)(t')] 2qHt ( t')t sL~+~l~l~, /t /t, / t/t=t + 2 ( 2-) LDH(l)(o) (t)' )+H (1) t1IHo(0)(t?)] (o )2 Dt' t (t )J + _) D(~) (1) L_)]t =t + 2y(1) [y(O)(o1+cy) - 2((j+a2)[H) (tT)]t2=t y(O)e 21 L + (1) [Y(~)(cr'-cr)-2(c2-al)] [ t ~(t) (3.64) 7(o )3 ( D2 i Jt'=t where primes denote differentiation with respect to 7 and where r=, l a;r', oj, and 2 are all evaluated at Y = y(). Hence, we see that, whenever v(O)+4v(l) is representative of a viscometric perturbation, the stress perturbation S(1) may be obtained simply by a direct substitution into Eq. (3.64). Let us now consider viscometric perturbations for which the derivatives D(l)/Dt of all primary quantities are zero. In such instances, Eq. (3.60) remains unchanged, but Eqs. (3.61)-(3.63) become (o )2 -I!) - D)2 +.(~)2 (1))(t') + 27y()y (3 65) LDt'o (0) (t' ) H( t (t ) + )]= (.66) Dtt and (o) (1) = o (3.6y) Dt

88 Since terms involving D(l)/Dt no longer appear in (3.65)-(3.67), the corresponding stress perturbation is found from (3.64) by dropping the term (a2-C ) D() (o ) We then have ~t ( + () =(1)(t')] = 2H4()(t,)]t S(1) + p(1)I F= + [o)(t )), t I =-0O 2 a,+ t t(0) +) (t)o + 2,[~la2) [ ~)(t')Ht(t' + t(t')]t =t:,(o)I, (o)-2-2 D() 1tt + 2-y (0)(Cry+ 7(o)2 LDt Jt =t t7(o) - 2 ( cl+ a2) ] [7H(O)(t1 ) ]2 /A- t'=t + (1) (o)( a)_-2 (a2 )] Dt )(t' ~(o)3 [ Dt'=t + 2y (1) (o) (t' ]t= t. (3.68) Thus, (3.68) will yield the appropriate stress perturbation whenever v(~) is viscometric and, in addition, Eqs. (3.60) and (3.65)-(3.67) hold for the perturbed flow. We next suppose that there exists a perturbation on a viscometric deformation history which may be written in the following form: Ht (t') = H(1)(t')] + [l)(t)nv (3.69) where [H(l)(t')]v satisfies (3.60)-(3.63) while [H(l)(t')]my does not; accordingly, the subscripts v and nv will refer here to "viscometric" and

89 "nonviscometric" respectively. Since the functional b of Eq. (3.28) is linear in Ht(1)(t'), the quantities [H(l)(t')]v and [el)(t')]nv may be treated separately, in so far as their respective contributions to the total stress perturbation are concerned. It therefore follows immediately that the stress perturbation due to [Ht (t')]v is obtained by substitution of [H(1)(t')]l for H 1)(t') in (3.64); similarly, if [Ht (t')]v satisfies Eqs. (3.65)-(3.67) and (3.60) (but not (3.61)-(3.63)), then the stress perturbation due to [H(l)(tl)]v is obtained from (3.68) on replacing Ht (t') by [Ht (t')]v. -— t tI

CHAPTER 4 TAYLOR-COUETTE STABILITY 4.1 THE PRIMARY FLOW We now apply the general theory developed in Chapter 3 to treat the current Taylor-Couette stability problem. It is supposed that a simple fluid is contained between two long coaxial cylinders of radii rl and r2 which are rotating at angular velocities Q1 and ~2, respectively. Cylindrical polar coordinates (r,O,z) are employed, with the axis of the cylinders lying along the z-axis. The undisturbed Couette motion consists of a velocity field v(o) with physical components given by v() =' v() )(r), (O) = v (4.1) Now, as we have indicated earlier, tangential motion between concentric rotating cylinders represents a viscometric flow; thus, from Eq. (1.54), we have (o) (o) H( )(t') = E()cos y(o)(t_-t) _ 1 sin y()(t-t') (4.2) A-t c ~~~y(o) t where, in the present instance, one may readily show that 7() = r d(v(~)/r) (4.) dr o 1 o E(O) (o) 1o 0]o (4.4) _v 2Loo 90

91 and ___ * ~~'/ (4-5) At o o o The corresponding stress distribution is given by S(o ) = (o ) I(o ) e(~)-s(~ ) = Gr -Szz a2(7 ) (4.6) Srr zZ S(o) S(o) =0 Gz =rz One can theoretically obtain an equation for the velocity of the primary flow v(O)(r) by substituting the above stress components into the equations of motion and applying the appropriate boundary conditions. In particular, we have SrG) = e(o)> = c/r2 (4.7) v(O)/r = 21 at r = r1 (4.8) v(O)/r = D2 at r = r2 where C is a constant. It is seen from Eq. (4.7) that, in general, the shear rate y(o) is a function of radial position within the gap between the cylinders. Therefore, since the viscosity function r(y) may vary with shear rate, there exists no general solution to (4.7) and (4.8). However, the solution to these equations can be approximated if the gap

92 width d = r2-rl is regarded as small compared to either r1 or r2. Under such circumstances, C/r2 will remain relatively constant as r varies between r1 and r2, and one finally obtains v(o) n1 + (r-rl) (n2-n1) (4.9) r - d y(o) A- (n2-1n)rl/d (4.10 ) and dv(o) Y() (4.11) dr 4.2 THE STRESS PERTURBATION S (1) Let us now assume, following Taylor,39 that the primary Couette motion v(o) is disturbed slightly by a steady, rotationally symmetric velocity perturbation ~v(l), for which 0+O. We write accordingly (1) u v(l) v v(1) w (4,12) r 9 where u, v, and w are functions of r and z only. In this section our objective will be to determine the stress perturbation S(1) due to the altered motion. One can first calculate the deformation history perturbation H(l)(ty) by substituting (4.1) and (4.12) into Eq. (3.54). The result is expressible in the form (H1)~(t') = [Ht( )(t')]v +4 [H(1)(t')]nv1+ [( )(t )]nv,2+ [H )(t )]nv,3 + [H(1) (t')]na a (4.13)

where, in cylindrical coordinates, 2 )r ar r A(vr) -r a(u r) av aW av 0 )r r _z -1 a(vr) -r (u/r) v sin y(~o)(t ) r1 —(v/r) r 0 (O)(t-t?)cos w(o)(t-t') 2 r o r r(u/r) -r 6(v/r) >w | 1 0 -r (3r 6r 6z _v _w lo 0 0 0 0 r -(vr) 0 1 (~)(t-t' )sin (~)(t -t ) 1 0 0 - 6(ur) _t r n(/r, 1 0-: 2()tt)iry~(-'41) O~0 1 0 O O O (4.15)

[Htl)(t' )]nv 2 <'- E ~ w A1>t nv.,2 0K O (4.1a6) 0 0 [HH 1J~~ ~(o [E 0(t 0)nv, > 0 a0 2cos(~)(t-t' )-cosy(~)(t-t' ~t n I 20 2 O | O' 2sinz (t-t' )-siny(Z)(t-t') 2 (o) () (o) (t-t (o0 2cos ()( t-t') -cosy (~)(t-t )2sin (t-t')-sin (~t-t') 2 2 (4.17) [H)(t)]nv (t-t) ) (~)(t-t') sin(~)(t-t') dr 0 0 sinZ( )(t-t' ) -cos (~)(t-t' ) O0 2 dr cos7(~ ) (t) (t-t' 0 o0 0 0 (4.18)

95 The subscripts v and nv are used here, as in Section 3.5, to denote "viscometric" and "nonviscometric" respectively; that is, one can easily (1) verify that [Hlt (t')]v satisfies Eqs. (3.65)-(3.67) and (3.60) with (1) = r, (v/r) (4.19) 6r Therefore, based on the discussion of Section 3.5, the stress perturbation S(1), due to [H(1)(t')]v, follows immediately from (3.68): Srr -SZz ]v = r + r vr r 2 ((1)1)v _ (u/r) (v/r) [S-$$SzzV = -r ( + r) 6r 6r [S(1)] = r _(u/r) (a-a a (v/r) (rn+T1y(o)) (4.20) [Ssr 2(0)v [s(l)]v + v a2 rz v =: (1) I _]6 [SZG v -= r7() ( l We are next faced with the problem of determining the remainder of the stress perturbation arising from the "nonviscometric" terms in H(1) (t,) As will presently be shown, it will be necessary to define several new material functions of the primary shear rate 7(0) (other than the three viscometric functions q, ao, and a2) in order to account for these stresses. Let us begin by considering the contribution towards S(l) due only to [H)t (t')]nv 1. Making use of Eq. (3.28), one notes that the appropriate stress perturbation Sn) will be given by ~nv, 1 wl egnb

96 t [S(1)+P(1)Inv 1 = [H(o)(t'),Ht (tt ) (4.21) t' =-00 where it is recalled that A is linear in [It (t')]nv 1 and satisfies the relation t t Q ~ 5 [H(~)(t' )~[H(1(t' )]nv l]' Q = t [QQ't~)(t')'Qt t' =-O t -00 -i[Ht l(t )n ].it] (4. 22) for any orthogonal tensor Q. We now choose Q in (4.22) to be the orthogonal tensor whose matrix is 1 0 0 Q, 1 0 (4.23) O 0 - Substituting (4.23) into (4.22) and employing Eqs. (4.2), (4.4), (4.5), and (4.15 ), one finds that t t Q ~}[H ~J(t,)~[H~( t'), E)t [H(1)(t')],[ H=)(t')[ 1]. t =00 t'=-' (4.24) Hence, the stress perturbation S v) due to the history [H(1) ]nv st %nv~l due [Ht (t?)]nv,l must have the property that Q S(1) s() (4.25) i 1 -nv 1,' - = nv(4 The matrix of Qi Sv(l). Q is readily calculated to be

97 S(1) S S "S — 3nvl ~ rr Sr rz S(1) (1) ((426) S$ - S.z (4.26) S (1) _ s) s(1) Srz Gz zz nvl It follows then from (4.25) and (4.26) that s(1) = [S (1) = o (4.27) rS z nv l z nv, l That is, S(1) must have the form.~nv 1 v~ZZ nv, nvrG nv,l r Or s(1) (1) [g SrO l=r ( r (4.29) o 0 s(1) rr nvl r Nowur purpose ince the new material functional s linear in (t it is seen from (415) that the components of S( wll be it to say proat this pointional to'nv1 that, by (ur)/rr;again moreover, the "constants" of proportionality will be material functions of the rate of shear y(o) Hence one can write (1)] l6(ur) c3 rG nv,1 r 6r o (1) (1)] _ (ur) (4.29) [e zz nvl r 6r 1 (1) _ (ur) rr zz nv,1 r -

98 A treatment similar to that presented above may also be employed to determine the stress perturbations S(1 S(1) and (l)arising from nv,2 Snv,3' -nv,4 arising from [H()(t)] nv [H)(t)]nv, and [HXl)(t')]nv 4, respectively. As a final result, we obtain the total stress perturbation S(1): S()+ p(1) = r (ur) + r (v/r) + 1 6(ur) dr(2 + Srr 6r r 2 r 6 r dr 7(O) _____ 2y___ _____ dy(~ ) _ _ (1) (1r1+ r ra + ( + u__ ~r r r dr o(o) r () r 6(u/r) ( r-2) + r 6(v/r) (,+ 1 ()(0) (l) r 27(0 a) 1 6(ur) dy(o) r 6r 7(o) dr +(l) = 2n v 02 U rz ar az y(o) az S 1 w u +) avL + v u C z 6r Y(o) az az (O) where the new material functions r11, 12, 13, T 43, 03q,43, \4 2, and X3 are all even functions of y evaluated at y = y(o) 4.3 THE DISTURBANCE EQUATIONS Appropriate expressions for the perturbations on the equations of momentum and continuity for this problem are found by substituting (4.1) and (4.12) into Eqs. (3.11) and (3.12), respectively, and employing cylindrical polar coordinates. We obtain (continuity) 1 (ur) + (4._1) r~~ -r (431

99 (momentum) (s) as(Q) S()() -2pv[(O)/rl rr + Srz + rr e ar az r Pu v() v() 1 6(r2Sr9 as z (4.2) ar L + rj = +.... (452) 0 (rSrz) _S_ r ar az Now, in order to facilitate the succeeding analysis, we make the usual simplifying assumption that the annular gap between the cylinders is smaall compared to their radii. The primary shear rate 7(o), which may then be presumed constant, is given by Eq. (4.10), and, in addition, the components of the stress perturbation reduce to (1) + p(L) a u +2 + V Orr + = [N+2 ] + 2 See + P W ar + x S(1) p() 23 (.33) W(l) _ U u 3+( r-2)/2] + r [1v+, (O Srz + o + Br U() aw aV C au ~z =e ~ +~ (0) x0+ — (1) _ -W by u C4 Sz ar - aZ aZ where terms involving d7(~)/dr have been dropped from (4.30). Furthermore, for the small-gap case, the equations of continuity and momentum become

100 au + aw = 0 (4.34) -2Pv + (r-ri)] S + pu(02- d)r - a + (435) d 6r 6z as(1) S(1) z zz o r az Substituting (4.33) into (4.35) and making use of (4.34), we then have -2pv 01+(02-01) (,r-1 ),= a_(1) a22u a2u t a2v 7(o) az2 r1 ar r1 ar - Crs- (&112)/ _ +U i(o0)] 62V(4.36) - 0j) + 62w + (23-14) a2 + a'2 2v )z 6r2 6I~Z2 7z(o) 5zr To reduce these equations to dimensionless form, we define the following quantities: X = r-r u = u/rld r2-ri d Z = z/d v* = v/1d = (02-s21)/j1 w* = w/Q1d (4.37) a = r1/d R = pd2Q/n(y (~o)) p. = p(1)/;=B(y(o)) T = -2caR2

101 The Reynolds number R and the Taylor number T are specified here in terms of the apparent viscosity rq(y(~)) rather than by some other "characteristic" point on the viscosity curve for the fluid. Any other definition of these quantities would serve only to introduce an additional and superfluous parameter into the analysis. Thus, we have established that the most appropriate and meaningful definitions of R and T are based on the apparent viscosity evaluated at the gap shear rate. From (4.34), (4.36), and (4.37) one now obtains au* + Hi_ 0 (4.38) +x + z.2v*(l+ctX)R = - a 2ux + _ + + 2 (~1-nl) au* L_ av* ia aX ra ax _ j u* = + y 2 + 2 r3( 2)7 2u* 4 &2u' 2R L B JC2 az2 L a 0) J X2 y(o) az2 O = _ 0 + 2w* (2n3-n4) a2w*+ a2v*X (4.39) az ax2 a Z2 7y(o) axaz Next, resolving the disturbance into its normal modes, we make the usual assumption that the disturbance is spatially periodic in the z-direction and write u* = ~e(X)sin eZ v* = -cV(X)sin cZ (4.40o) w*x = d cos EZ dX p* = cP(X)sin cZ (4L41)

102 where c is the dimensionless wave number of the disturbance. One can readily verify that (4.40) satisfies Eq. (4.38), the equation of continuity, identically. Substitution of (4.40) and (4.41) into (4.39) then yields 2V(l+caX)R + -a D2V- 2 2 = -DP + ILa D2= _ c2 (- ) D* n IZ (o) r l rTa a _' __. D2V - - T + [aY3-(Cl+C2)/21 D2 (- 4 C 2R~.Y. 0 BDo) (4.43) and 0 = -c2P + D3V - D(2*34) (c2Df - 2 2DV (4.44) where D - d/dX We can now eliminate P between Eqs. (4.42) and (4.44); the final disturbance equations are then given by [B1D2-2]V = 2D (4.45) and LD4'4Z2D2 + C3~ 2D+-5Se =]2 2LR(l+CX)+0D2+02 + aj V (4.46) where the fis are rheological parameters defined as* *For convenience, we henceforth drop the superscript on the primary shear rate y(o) with the understanding that y and y(o) are now one and the same.

103 Pi d In[ry] (4.47) d In y =2 F3 (ar+a2)2 (4.48) =3 = 4/ Y (4.49) =4: (2T3r"2 -14)/T (4.5o) AL =c 4/1 (4.51) ~7 = d2/n (4.51) depending on the shear rate y4 which are necessary to define the TaylorCouette stability of a simple fluid in the narrow-gap limit. Three of these are the viscometric functions l, CL, and 92, but the other five 29, 53, 549 85, and 59 are entirely new. One notes, from (4.48) and (4,49), that 52 and 53 are odd functions of y while, from (4.50), (4.51), and (4.55), we see that eqa 5t, and 5h are even functions of y. Equations (4.45) and (4.46) represent two ordinary, homogeneous, differential equations in two unknowns, a, and V. Together with the homogeneous bouna5ry conditions,

lo4 = D = V = O at X = 0,! (4.56) they determine a characteristic-value or eigenvalue problem for the Taylor number T, as function T(e) of the dimensionless wave number c (for fixed a and a); the smallest value, Tc, of T(c) determines the critical conditions Tc, cc at which instability sets in. 4.4 LOW-SHEAR-RATE APPROXIMATION There are presently no known methods, either theoretical or experimental, for determining values for the five new material functions 52 to P5 and P9 at finite rates of shear. However, as we shall see below, the variation of these quantities in the limit as y-*O can be predicted, to terms of first order in y, by making use of simple-fluid theory. We recall from Section 1.5.1 that the fluid model considered by Datta,ll S + = czA1 + 2 = A + + 3A (1.66) has been shown to describe the limiting behavior of simple fluids with fading memory at extremely low rates of deformation. Here = lim Ir(y) 7+0 a2 = lim 1(y)/72 (1.67) y+0 a3 = lim (C2-C1)/272 Y+0 One finds from Datta's disturbance equations that

o105 P2 = P3 = (a2+a3)y/a1 (4.57) 4= 2, = 12, 5 = 0 Therefore, in the limit as y+O, we have =2 = 3 = (a1+o2)/2n7 (4.58) B4= 2 = 1, = 0 An obvious suggestion which one might put forth at this point would be to approximate the material functions B2 to B5 and B9 by Eqs. (4.58), even at finite rates of shear. If such a procedure were valid, then the solution to the disturbance equations would be governed entirely by the viscometric functions u, l,, and a2; hence, knowledge of these three functions would suffice to determine the flow stability of a particular fluid a priori. One may now conclude the following: If at the point of instability, the approximation (4.58) were reasonable, absolutely nothing new concerning the fluid of interest (beyond i, ol, and CT2) would be learned from a stability experiment. Conversely, if the approximation were not valid, then new information would be obtained. In fact, any difference between the results of a stability experiment and the solution to the disturbance equations with the approximation (4.58) would presumably be attributable to deviations from (4.58). The validity of the "viscometric hypothesis" is likewise contingent upon whether Eqs. (4.58) represent a reasonable approximation. One will recall that this hypothesis requires that Eq. (1.57) describe motions,

o106 such as the present one, which are close to viscometric flows; it is readily verified that Eqs. (4.58) are precisely what would have resulted if Eq. (1.57) had'been chosen as the original Theological model for the stability analysis. In certain cases, a stability experiment can be used to advantage in determini:ng which of two fluid models is more appropriate for representing the constitutive'behavior of a particular viscoelastic fluid. The primary requirement for these models, of course, should be that they both adequately describe the experimental data on ~, Al, and 02; in this way, the material functions BL and 56 to 58 would be the same for both. If the models then differed sufficiently in their predictions of f2 to 5s and Se the corresponding stability analyses would yield different results. Based on a comparison of the theoretical findings with experiment, a choice of the'best model could then be made. 4.5 PREVIOUS ANALYSES As has been pointed out earlier, all previous theoretical treatments of the small-gap Taylor-Couette stability problem have involved consideration of the'behavior of various idealized simple-fluid models; these treatments will henceforth represent special cases of the present, general analysis. In order to illustrate more clearly the relationship'between past studies and the current work Table 4.1 has been prepared. This table presents a listing of the eight material functions ~, 0l, c2, B2, 3, F4. 55 and 9 for each of the previous investigations.

TABLE 4.1 MATERIAL FUNCTIONS IN PREVIOUS ANALYSES* Investigator Model q al a2 52 s 84 85 89 Chandrasekhar5 Newtonian o 2 1 Eq. (1.19) ReinerGraebel21 Rivlin al a2 22 a22 (aC+a 2)/2my (C3al+2)2ry 2 1 0 Eq. (1.20) Bingham 2 Graebel22 Plastic ~ 0 0 0 0 2 1 0 Thomas & Walters B' Walters40 Eq. (1.24) To 21oy2 0 (cl+Cr)/2f Y (al+2)72y 2 1 0 Walters A' Chan Man Fong6 WaltEq. (1.2) o 0 -2Koy2 (1a+52)/2rIy (z+la2)/2rq7 2+3SoY2/ o 1+3S0o2/o ~ Fluid of Dattall Second Grade al a2y2 (a2+2a3)72 (cl+a2)/2,y (Cl+a2)/2]7 2 1 0 Eq. (1.66) Modified (al+a2 /2 y Walters4l Eq. (1.24) 2KOy2 0 (CraIc2)/2r7 (1i+92)/2q7 2 1 Ericksen 2 3So~_/to** Leslie25 icksen +3 nln2(t2n2 nln2(42n 0 n(- (44+n21 +n42 0 Eq.Leslie An(1sotr0) + 2 +24s)y +24)y 4/r)/n2 2nBooL2)/2a +22nIt)/ 0 Eq. (1.72) *It is interesting to note that, due to the simplicity of many of the models considered, the approximations shown in Eqs,(4.58) are found to occur quite frequently in this table. **S-o | f 2t(t)dt 0

108 One study not included in Table 4.1 is that of Giesekusl8 who attempted to generalize the theoretical analysis of Dattall for finite rates of shear. He substituted 1 = a2 = a/y2 (4.59) a2+2%3 = J2/72 into Datta's disturbance equations and obtained Pi = 1 6 = - 2/ y P2 = 83 = (G1+92)/217 P7 = 92/17 (4.60),4 = 2 B8 = (2aj+crp2)l l P5 2 = o By comparing these equations with Eqs. (4.47)-(4.55), one notes that most of the expressions in (4.60) are clearly incorrect, except of course at y = 0. 4.6 SPECIAL CASES In this section, we briefly discuss two special cases. The first is that of negligible inertia where T-+O and the other is plane Couette flow where a+oo. 4.6.1 Negligible Inertia There are numerous ways in which the disturbance equations for the current study could have been written. However, in the form we have pre

109 sented them (Eqs. (4.45) and (4.46)), the terms due to inertial effects (those containing the Taylor number T and the Reynolds number R) have been completely separated from the "rheological terms" (those involving the 5's). 18 Giesekus8 was the first investigator to consider the possibility that, for a particular viscoelastic fluid, instability might occur at an extremely small Taylor number, where inertial effects are negligible. Under such circumstances, Eqs. (4.45) and (4.46) become [P1D2-_2]V = [ 2D2 -S3C2] [D4 P4(!2D2 + f C2D+G]5E = G2D 6+ ]a+(4.6 1) Apparently, a modified Taylor stability criterion could no longer be regarded as applicable to these new disturbance equations since they do not even contain the Taylor number. Instead, an entirely new criterion results, in which the rheological coefficients determine the point of instability. 4.6.2 Plane Couette Flow (Simple Shear Between Parallel Planes) The disturbance equations which describe this situation are obtained by allowing a-+oo while y is held fixed in Eqs. (4.45) and (4.46). We then have [ D2-c2]v = [-R,+ 2D2-f3c2Ls (4.62) [D4 _4~2+8564 ] ~ = ~2 [6DI2+T7E2 ] V

110 where R = pd2yl. (4.63) For the case of a Newtonian fluid, these equations yield only the trivial solution, ~ = V = 0; however, for certain idealized viscoelastic fluids, Giesekusl8 has shown that solutions other than the trivial one are possible. 4.7 SOLVING THE DISTURBANCE EQUATIONS In the present work, three methods have been used to solve the characteristic-value problem defined by Eqs. (4.45), (4.46), and (4.56). The first and second of these were analytic techniques, and were employed only to obtain approximate expressions for the Taylor number as a function of the dimensionless wave number; the third method, a numerical technique, was used to obtain more precise results. 4.7.1 First Analytic Method The first analytic method, which can be attributed to Chandrasekhar5 is illustrated by considering the solution of the simplified equations [P1D2-e2]V = T + 22[D2-_2] (4.64) [D2_E2]2t = c2[2R(1+X)+6D + +7E + ajv (4.65) which result when the approximations shown in Eqs. (4.58) are made.

111 We begin by writing V = ~ Am sin mTX (4.66) m=l thus satisfying automatically the boundary conditions on V. Substituting (4.66) into (4.65) and solving for p, one then obtains co, = -2RE2 r Am rm cosh cX+Dm sinh eX + EmX cosh eX m=1 [ (m )2+C2 ] + FmX sinh cX + l+cyX + 56(mit)2 -7Cj sin m+TX + Im cos m (4.) 2R where 4{a Im: m.r - 2 (4.68) and where Cm, Dmn Em, and Fm are determined from the boundary conditions on ~: D = Dp = O at X = 0,1. (4.56) Equations (4.66) and (4.67) are next substituted into (4.64) to yield 00 oo, Bm[pl(m )2+C2][(mIt )2+]2 BmlC 42R cFcosh [X m=l T2 m=lT + (Dm + R cE) sinh eX + EmX cosh ~X + FmX sinh eX + (1 22R) ( + ax + p6(mit)2p7 ) sin mitX + Jm cos mXj (4.69) T / 2R where Bm = Am/[(mit)2+c2]2 and (4.70) J. = 2? R [(mi)}2 + + 4ozm 2iR

112 Multiplying this equation by sin n.tX and integrating between X = O and X = 1, we then obtain an expression of the form 00 K BmKmn = 0 (n=l,2,... (4.71) m=l which represents an infinite system of linear homogeneous equations in the variables Bm. The condition under which this system has a nontrivial solution is that the determinant of the coefficients Km should vanish. Since in the present treatment, we are seeking only an approximate solution for T as a function of E, we set only the first element K11 of the K-matrix equal to zero. The first approximation, after considerable reduction, is then given by [g2+e2 ]2[ 1t2+ c2] T= E2[l+4a/2] 1 8- c22E] (e+coshe:' (2] (+E2 + 2__/ D7El - I+/ ] [e2+C2] (e+sinhe 2R(l+a/2 (4.72) 4.7.2 Second Analytic Method The second analytic technique employed was in many respects, quite similar to the first. It was essentially a Galerkin method (see e.g., Crandall10), and has been used by Walowit43 in treating the Newtonian stability problem in the case of a radial temperature gradient. Consider here the complete disturbance equations: [1D.2-e2v V = + 2D2 - B3E] (4.45)

ll3 LD4_ 4C2D2 + G9:2D+PEj v = E2 L2R(1+X)+P6D2+P7E2 + a V (4.46) Let us write 00 V = - AmVm and (4.73) 00 = ) Bmfm m=l where Vm = X(l-X)mmand (4.74) lm = [X(1-X)]2 Xm-l One notes that all boundary conditions on both J and V are thus satisfied. Substituting (4.73) and (4.74) into (4.45) and (4.46), we then have 0oo LAm(lD2-c2 )Vm + Bm [R+ P2D2 - 3ca 2j = 0 m=l (4.75) PBmD4 42D2+_42D 4- ~m+Ame2R(l+cX)+]6D+7c + VVj3 = o m=l These equations are next multiplied by 4n and Vn, respectively and integrated between X = O and X = 1 to yield expressions of the form

114 00 [CmnAm+DmnBm] = 0 (n=1,2,...) m=l (4.76) Equations 6) represent a doubly infinite set of linear homogeneous m=l Equations (4.76) represent a doubly infinite set of linear homogeneous equations in the variables Am and Bm. By setting the determinant of the system equal to zero, we obtain the required characteristic equation for T. As a first approximation, we have* 28 [504+12,4.2+p5c4] [1Op1+C2] 27 28 (7 2 [1+/2] 1- R ( p2+3C + T 0-67 IlLT t- 2R(l+a/2) J This same method may also be used to solve Eqs. (4.61) for the case of negligible inertia and Eqs. (4.62) for the case of plane Couette flow. One ultimately finds that the first approximations to the characteristic equations for these situations are given respectively by 28 [504+12p4c2+P5c4] [lO1+21 - 1 (4.78) 27' c2 2 28 2 IL 28) 3 3 and in 8 -524+12r, +DSC41[lOpl+C]_ 8 2+]e 2 (4.79) 27 C2FY72E 2_ 3 4.7.3 Numerical Method The numerical method used in solving the disturbance equations was *The similarities between Eqs. (4.77) and (4.72) should be noted.

115 apparently first successfully applied by Beard, Davies, and Walters.l It consists of replacing the current boundary-value problem with a related initial-value problem. In any eigenvalue analysis, there is an arbitrary amplitude associated with the eigenfunctions. This indeterminancy can be removed in the present instance simply by imposing the additional initial condition D2 = 1 at X = O. (4.80) We next assume values for D3* and DV at X = 0: D3 ]]X= = G1, DV]X=O = G2 (4.81) where G1 and G2 represent first approximations. An initial-value problem is now defined by Eqs. (4.45) and (4.46) with boundary conditions * = DV = V = o, D2 = l, D3 = G, DV = G2 (4.82) at X=O. Once the coefficients in (4.45) and (4.46) are specified, one can integrate the system of equations numerically by using a standard Runge-Kutta technique (see e.g., Carnahan, Luther, and Wilkes3) from X = 0 to X = 1. In the present work, the method of obtaining these coefficients was as follows: (1) Making use of Eq. (4.10), rewrite the definition of the Taylor number T in the form T = -2Ia d2 - _ 2 (4. 83)

where T is an appropriately chosen characteristic time for the fluid, qo is the zero-shear viscosity, and EA is the elasticity number, a dimensionless constant: E = r. (4.84) pd2 (2) Define the modified Taylor number Tm and the modified Reynolds number Rm by T = T 2 2 T -2aaR2 (4. 85 ) m = T~JV) = a LE Z-2c=a' (3) Express the fluid viscosity q in the form r = io/5o (7T) (4.86) where Po defines a new dimensionless function. (4) Express the material coefficients.i(y) (i=1,...,9) in the form ~i = i( 7). (4.87) (5) Assume a value for the modified Taylor number Tm: Tm = G3. (4.88) (6) Calculate the corresponding modified Reynolds number Rm: Rm = [] = G/-2aa (4.89) [ a] [Ef]

117 (7) Calculate the dimensionless shear rate yT 7T = [Ca][E]Rm. (4.90) (8) Evaluate the material functions ~i(7T) (i=0O,...,9) at the dimensionless shear rate. (9) Calculate the actual Taylor number and actual Reynolds number from the equations T = Tm,2 (4.91) R = Rmo The solutions for a, D*, and V obtained by integrating the disturbance equations from X = 0 to X = 1 will depend, presumably in a continuous fashion, on the values chosen for the three quantities G1, G2, and G3. Of course, for arbitrary choices, the boundary conditions (4.56) will not necessarily be satisfied at X = 1. Therefore, one must determine appropriate corrections to G1, G2, and G3 in order that these conditions may be met. Such corrections are obtained in the present work using a Newton-Raphson technique;3 in particular, we write, to terms of the first order, 6(D*) 6(Dt) 6(Dt) = Dr(92 a__ a__ aG X=1 _ _ V av, ac v aGv X=l

where 6Gi (i=1,2,3) represent the required changes in the Gi. Iteration on GI, G2, and G3 is repeated until satisfactory convergence of *, Dt, and V at X = 1 has been achieved. The nine partial derivatives in Eq. (4.92) have been evaluated herein by a finite-difference process. That is, each of the quantities G1, G2, and G3 was individually given a small perturbation, and the corresponding differences in V, D*, and V at X = 1 were noted. By making use of the methods outlined above, one can determine the Taylor number for any value of the dimensionless wave number e. The smallest T over all C represents the critical Taylor number Tc at which vortex cells first appear in the flow. In this study, Tc has been found numerically by use of the "golden-section" minimum searching technique (see e.g., Wilde46). Numerical calculations have been performed on the IBM 7090 digital computer located at The University of Michigan Computing Center. A short description of the programs and subroutines employed is presented in Appendix D. 4.8 QUALITATIVE DISCUSSION Let us now examine the qualitative effects of the various rheological coefficients B1 to P9 on the critical parameters Tc and cc. For this purpose, we consider Eq. (4.77), our approximate solution to the complete disturbance equations.

119 One notes immediately that the variables P8 and P9 are absent from (4.77). Therefore, to a first approximation, s8 and P9 have a relatively small influence on the stability criteria. This fact has been verified numerically; in one numerical case, P8 was varied from -2.4 to +2.4 while, in another case, 9s was varied from 0 to 2; in both instances, it was found that Tc and ~c changed by less than 0.05%. Hence, we see that, for moderate values of 58 and P9, the critical parameters Tc and cc are determined primarily by six material functions (Q, a2, and 52 to B5) of the shear rate rather than by eight. The corresponding disturbance equations for this situation are given by [PID2C2]IV = L + 2D2 - P3C>2 [D4 -4c2D2+5 E4] = C2[-2R(l+aX)+P6D +7c 2]V In studying the effects of the other material coefficients P1 to In on the stability criteria Tc and Cc, we begin by noting that Eq. (4.77) may be written in the general form T = f[Pi(T),TR(T),c] (i=l,...,7) (4.94) for fixed values of the parameters a and a. Taking the derivative of T with respect to c, we have dT If dDi dT f dT +f dR dT +f. -—. —+ — + -— + - (4.95) d- - i dT dc aT d aR dT d c or equivalently

120 dT = r[af/a l. (4.96) dl -_ rfdi + af +_ 3 dR| L LE pi dT aT R dTJ The critical Taylor number Tc and the critical wave number cc are found by setting dT/dc equal to zero. From (4.96), one obtains -_ [pi(Tc),TcRc(Tc),cc] = 0 (4.97) while from (4.94), Tc = f[1i(Tc),TcRc(Tc),c1] (4.98) where Tc = -2oaR. (4.99) Equations (4.97) and (4.98) represent two simultaneous equations in the two unknowns Tc and ec. Let us now allow Pi(Tc) to vary independently of Tc* and write d af [pi(Tc)Tc Rc Ec ] dpi(Tc) -, ra_2f 62f dTc 62f dR dTc 62f dc =_ _a _f + ++ - 0 6aEa:i E:cT dpi(Tc) c~eR dTc dpi(Tc) aC2 d~ i(TC)0 (4.100oo) and dTc _af af dTc af dRc dTc af dEc d~i(To) La:i aT d:i(To) aR dTc d~i(To) ae d~i(T~)JTo,~c (4.101) *This can be accomplished by varying the material parameters.

121 Substituting (4.97) and (4.99) into (4.101) and solving for dTc/d~i(Tc), we have ~~dT0: [3~~f/8~ 2; ~(4.102) di (Tc) R c _ 1_ f +R 2T Tc c Moreover, from (4.100), 6a2f +a2f 62f Rc \ dTc dc - LEaTi +ka'aT + cR 2Tc) dpi(Tc)J j (4.103) d3 i(TC) [C/a[2f/ JTCYEC The terms 6f 6f RC an' 2f a2f R 1 -+ -R 2T TcCc R TCTEC, adLCT +C6R 2Tc JTcec in these equations can be evaluated using Eq. (4.77). One obtains (F P6-P7 2) 2Rc 28 f af Rc 1]1 2RC(1+CK/2) Tc (- 2+2 )3 aT 6aR2HTcJTc,% 28 2 +72) 2R 2Rc(l1+/2) Tc 3- Cc (4.o104) and 2Rc 57_ C 62f +2f Rc r 3EC 2Rc(l++C/2) I~+ -- 1 = _ I R + I LWc:T 6E:R 2TcJTC,, K 2Rc 28 +2 ) 21 + 2 6-5 82) L Tc (5- 2+P3) E + 2R(1+ ~ (4.105) (4.105)

122 Now, in the present work, the expression given by (4.104) was, in all instances, small in comparison with unity. Thus, it appears from (4.102) that the sign of dTc/dpi(Tc) was determined entirely by the sign of [ f/6Si]TcCc. This quantity has been calculated for each of the cases i=1,...,7 and the final results are presented in Table 4.2. TABLE 4.2 EFFECT ON Tc OF A SMALL CHANGE IN Bi(Tc) Parameter dTc/dpi(Tc) 51 Greater than zero p2 Greater than zero F3 Greater than zero p4 Greater than zero P5 Greater than zero e6 Less than zero P7 Greater than zero In investigating the influence on cc of a small change in Bi(Tc), we consider here (for simplicity) only deviations relative to the Newtonian case, ~l = I, P2 = s3 = 56 = = O, 4 = 2,J 5 = 1. Under these circumstances, from (4.103) and (4.105), one has dec [62f/6cpi1l (4.106) dai(Tc) [a2f/c2] JTc, Cc Recalling now that at the minimum point of a function, the second derivative is positive, it is seen that ro2f 0. (4.17)

123 Therefore, for [a2f/ai~]lTcEc > O, we have that dec/d0i(Tc) < 0, and vice versa. Table 4.3 presents the results of a study on the sign of f/c ]Tcc for each of the cases i 1,.'..7. TABLE 4.3 EFFECT ON cc OF A SMALL CHANGE IN Pi(Tc)* Parameter dcc/di (Tc) Pi Greater than zero P2 Zero p3 Less than zero P4 Less than zero P5 Less than zero P6 Zero I7 Less than zero *Relative to Newtonian case. 4.9 PRESENTATION AND DISCUSSION OF NUMERICAL RESULTS In this section, the numerical results of the current investigation are presented. Consideration is given first to the Newtonian case, where the values obtained here for the critical parameters Tc and cc are compared with those found by other authors. Next, we examine the stability'behavior of a fluid defined by the general constitutive model S + pI = 2jl(7)E. Thirdly, we employ the approximations presented in Eqs. (4.58) to solve the disturbance equations for the case of two fluids whose viscometric functions r, a, and e2 have been determined experimentally by Huppler; 4 the current experimental and numerical results are then compared and certain definite conclusions are reached concerning

124 the validity of the "viscometric hypothesis." Finally, we consider numerically the effect on Tc and Ec of varying the parameters P2, P3, 54, and P5. 4.9.1 The Newtonian Case The disturbance equations describing this situation are given by [D2-E2]V = T _ (4. 108) [D2-_2]2x = -2c2R(l+aX)V which may be expressed in the form [D2-c2 Iv = (4.1o9) [D2-_2] 2 = _TE2(l+aX)V if the variable 4 is suitably redefined. It is observed from equations (4.109) that the Reynolds number R = 4T/(-2oaa) (or equivalently the radius-to-gap ratio a = rl/d) does not occur explicitly in the small-gap Newtonian stability problem. Thus, for this case, the stability criteria Tc and Ec are functions only of the parameter a. In order to illustrate the accuracy of the present numerical procedure, Table 4.4, which contains a comparison of the current findings with those of several other authors (for the Newtonian case), has been prepared. One notes the excellent agreement between the present and past results.

125 TABLE 4.4 VALUES OF Tc AND cc FOR NEWTONIAN CASE Critical Investigator Parameters a = -0.5 c = -1.0 a = -1.5 -2.0 Thomas and Tc 2285 3404 6431 Walters40 Cc 3.12 3.12 3.20 Chan Man Fong6 Tc 2290 3414 Cc 3.10 3.10 Beard, Davies, Tc 2275 and WaltersI c Chandrasekhar5 Tc 2275 3390 6417 18677 Ec 3.12 3.12 3.20 4.00 Harris and T - 3390 6414 18663 Reid23 ~c -- 3.13 3.20 4.00 This study Tc 2275 3390 6414 18669 e: 5.12 35.13 3.20 4.00 4.9.2 Case of a Fluid Defined by S+pI = 2(y)E A fluid obeying the general constitutive law S + p1 = 2(7)E (4.110) exhibits neither normal stresses nor elasticity; nonetheless, this model provides a worthwhile basis for comparison with real viscoelastic behavior since it possesses a viscosity which depends on shear rate, a property common to all viscoelastic fluids. The usual Bingham plastic model, the Ostwald-de Waele "power-law" fluid, the Reiner-Philippoff fluid, the Eyring fluid, and the Ellis fluid are all special cases of

126 (4.110). Hence, the need for any future stability treatments using these or other such simple models should be obviated by the present analysis. By redefining the variable V, it is possible to express the disturbance equations for (4.110) in the form [PjD2-E21V = (4.111) [D2_-212 = -To2[l+aX]V We see then that the stability criteria for this model will be governed by the two variables a and P1. Owing to the limited amount of computer time available in the present work, numerical calculations have been restricted to the case of a = -1. However, suitable approximations for Tc and Ec for other values of a can be determined using Eq. (4.77). In particular, it is found that, for -1 < a < O Tc(a ) Tc(al) Tc(-~Igf) Tc(-1,1) and (4.112) F(C(aOL) 1 EC(- -.l) where Tc(a,l) and Tc(-1,l) refer to the critical values associated with the Newtonian case. Table 4.5 presents the results of the nxrmerfc Z c7C~JZo ~rferQeb kerein for the model defined by (4.110), and Figs. 4.1 and 4.2

127 show plots of Tc and cc vs. P1 (for a = -1). It will be observed that both Tc and cc increase with increasing P1, in qualitative agreement with the predictions of Tables 4.2 and 4.3. TABLE 4.5 NUMERICAL RESULTS FOR A FLUID DEFINED BY (4.110); cZ = -1 B1 Tc Cc 2.00 5009 3.50 1.00 3390 3.13 0.99 3373 3.12 o.98 3356 3.12 0.97 3338 3.11 o.96 3321 3.10 0.95 3304 3.10 0.90 3217 3.07 0.80 3040 3.01 O.6o 2672 2.85 0.40 2276 2.74 0.20 1822 2.30 0.10 1542 2.00 From the present experimental findings, the values of P1 occurring at instability have been calculated for each of the polymer solutions investigated, and these values have then been used in conjunction with Figs. 4.1 and 4.2 in an attempt to predict the experimentally observed stability behavior. Table 4.6, which contains a comparison between theory and experiment, presents the corresponding fractional reductions in Tc and cc. One notes that the theoretical results based on Eq. (4.110) are in qualitative agreement with experiment, in that the critical parameters both decreased with polymer concentration. However, the predicted decrease in ec was, in all cases, somewhat smaller than observed.

128 6000 F" 5000 D 4000 z \ 3000 I_ 2000 ct 1000 0 1.0 2.0 Fig. 4.1. Plot of T~ against p1 for a fluid defined by (4.110); ac = -1.

129 4.0 3.0 w 2.0 -J oI... 0 1.0 2.0 Fig. 4.2. Plot of Ec against P1 for a fluid defined by (4.110); c = -1.

13o TABLE 4.6 FRACTIONAL REDUCTIONS IN Tc AND cc Experimental Theoretical* Substance Conc.(Wt.) TN TTcN 7MT CMC 0.5 1.0 0.97 1.0 1.0 1.0 0.92 0.93 o.g94 0.97 1.67 o.88 0.90 o.89 0.95 2.5 0.45 o.85 o.80 o.91 P-75-XH CMC o.183 1.01 o.88 o. 82 0.92 o.367 o.95 0.75 o.80 0.91 0.55 0.79 0.71 0.70 0.89 Natrosol 250-H HEC o.48 o.86 0.82 0.80 0.91 0.80 0.74 0.72 0.67 0.85 1.o o.63 o.6g o.64 o.82 *Based on Eq. (4.110) as a constitutive model. 4.9.3 Approximate Analysis Metzner, White, and Denn29 have suggested that the equation S + pI = (y)Al + i A2 + (&2&1) A2 (1.57) /S- ~ ~ 2 2 72' or equivalently S + pI = 2n(y,)E + 2( ) E+ 2 (1.56) "'''Y2 Y2 jt might be useful in describing motions, such as the current one, which are extremely close approximations to viscometric flows (viscometric hypothesis).* If this were the case, then, as has been indicated earlier in Section 4.4, the "unknown" material functions B2 to B5 and B9 would -be given by the low-shear-rate approximation (4.58), and, in addition, the Taylor-Couette stability'behavior of a fluid would be governed entirely *The reader will recall that (1.57) is applicable in the strictest sense only to purely viscometric flows.

131 by r,, ca, and c2. In order to explore these possibilities fully, the following course of action has been taken. First of all, the viscometric data of Huppler24 on solutions of 0.50% P-75-XH CMC and 0.9% Natrosol 250-H HEC have been represented (in the power-law region) by empirical expressions of the form n = No/17YTJ (l = CEo0 17 I//7T IP1 (4.113) 2 = c2o01Yl/IYT P2 where T is a characteristic time,* f0 is the zero-shear viscosity, and c1, C2, PO, P1, and P2 are dimensionless material constants; the particular values employed for these quantities are shown in Table 4.7. TABLE 4.7 MATERIAL CONSTANTS USED TO APPROXIMATE HUPPLER'S24 DATA Substance T(sect) ro(poise) c1 c2 Po Pi P2 P-75-XH CMC 0.76 7.42 0.275 0.0812 0.446 0.13 0.19 Natrosol 250-H 0.148 13.7 1.99 0.127 0.61 0.423 0.375 Next, based on Eqs. (4.58), (4.113), and Table 4.7, the disturbance equations for the two fluids of interest have been solved for the case *In the present work, the characteristic time T has been taken to be the reciprocal of the shear rate occurring at the point of intersection of the zero-shear viscosity line with the power law line in Fig. 1.2. It is found then that for?7TI < 1, rl _0 while for lyr i > 1, ~ is approximated by the power law.

132 where a = -1 and a = 12.* The numerical results obtained are presented in Table 4.8 along with the predicted small-gap values of To and Cc, determined from the stability tests. TABLE 4.8 COMPARISON OF APPROXIMATE THEORY WITH EXPERIMENT Theoretical Experimental Substance T Tc Cc Tc ec P-75-XH CMC 1771 3.73 2900 2.2 Natrosol 250-H 1001 5.25 2350 2.2 As can be seen, a considerably greater reduction in stability is found theoretically than was actually observed. Even more noteworthy is the fact that the critical wave number cc decreased from its Newtonian value of 3.13 in the experiments while, in the calculations, this quantity was greater then Cc N. One could argue that such discrepencies in results might be attributed to small deviations of Eqs. (4.113) and Table 4.7 from the exact viscometric functions. However, calculations have shown this not to be the case. For example, it was necessary to actually reverse the sign on 92 before the experimental and numerical findings could be brought into even rough agreement.** Thus, we have *It is recalled that values of a = 12 and a = -1 correspond to the current stability apparatus and experiments. **Dr. B. Duane Marsh27 has indicated in a private communication that Huppler's results for a2 have been checked by an independent technique, and that the data from the two methods agree to within 10* for HEC solutions, giving in each case the same sign.

133 demonstrated that Eqs. (4.58) do not represent a reasonable approximation for 52 to P5 and 59 at finite rates of shear and, furthermore, that the "viscometric hypothesis" is unacceptable in the present instance. 4.9.4 Effect of Varying 2, 53, 54, and 85 In order to illustrate the fact that agreement between the present theoretical development and experiment is possible for certain values of the parameters 52 to P5,* these quantities have been allowed to vary, while ri, 1,, c2 (and thus 51, 56, 57, and P8) have retained their functional dependence as given in the previous section. The disturbance equations for the two fluids of interest were then solved numerically for each of the cases listed in Tables 4.9 and 4.10. From the first enTABLE 4.9 EFFECT OF VARYING P2 to P5 FOR P-75-XH CMC P2 P3 P4 - 5 Tc Cc O 0 2 1 2067 3.47 -1 0 2 1 1863 3.50 1 0 2 1 2290 3.44 O -1 2 1 1798 3.95 0 1 2 1 2301 3.15 0 2 2 1 2511 2.92 -2 1 2 1 1882 3.22 2 -1 2 1 2223 3.86 0.5 0.5 2 1 2301 3.28 -0.5 -0.5 2 1 1838 3.70 -1 1 2 1 2082 3.19 -3 3 2 1 2032 2.85 O 0 2 2 2349 3.00 0 0 2 0.5 1851 4.03 O 0 3 1 2407 3.26 0 0 1 1 1706 3.71 *In Section 4.8, we have shown that the dependence of Tc and cc on 59 is very weak.

134 TABLE 4. lo EFFECT OF VARYING P2 TO P5 FOR NATROSOL 250-H P2 P3 P4 P5 Tc Cc 0 0 2 1 1925 3.36 1 0 2 1 2133 3.31 -1 0 2 1 1738 3.40 o -1 2 1 1683 3.88 0 1 2 1 2133 3.03 0 2 2 1 2314 2.79 1 -0.5 2 1 2007 3.54 2 -1 2 1 2080 3.8o -2 1 2 1 1748 3.11 0.5 0.5 2 1 2138 3.15 -0.5 -0.5 2 1 1718 3.60 -1 1 2 1 1932 3.o6 -3 3 2 1 1869 2.72 o 0 3 1 2228 3.12 0 0 1 1 1600oo 3.64 O 0 2 2 2158 2.89 0 0 2 0. 5 1743 3.93 tries in these tables, it is apparent that replacement of P2 to P5 by their respective Newtonian values, of 0,0,2, and 1, yields considerable improvement over the theoretical results shown in Table 4~8. Figures 4.3 and 4.4 indicate graphically the effect on Tc and cc caused by varying the parameters 52 and B3, with 4 and B5 being held fixed at 4 = 2, S = 1. Similarly, Figs. 4.5 and 4.6 show the dependence of the stability criteria on P4 and P5s with P2 = P3 = O. In these figures, the experimentally predicted small-gap values of Tc and cc are indicated using the symbol ~, and the numbers in parenthesis represent points in the original 82-s3 or 54-5s planes which have been mapped into

.0,1)(05,0.5) (1,O) 2000 l (-2,1) (i..0)(.Q5..Q5) (0,-I) H 10002 2.5 3.0 3.5 40 Ec Fig. 4.3. Mapping of the.2-3. plane into the Tc-Ec plane for a 0.5% solution of P-75-XH CMC; a = -1, a = 12, P4 = 2, 5 = 1) 9 = 0.

3000 ~o (0.2) (0o,o0) OI) (0.5,0.5)(I,0) (2,-) O,-0.5) Tc 2000 (-3,3) (I, I (-2,1) (-1,0) (-0C.5 0.5) 21000 2.5 3.0 3.5 4.0 cc Fig. 4.4. Mapping of the P2-P3 plane into the TC-EC plane for a 0.9% solution of Natrosol 250-H; a = -1, a = 12, P4 = 2, 15 = 1, B9 = 0.

(2,32) TC 2000 (2,0.5) H 1000O O2 3.0 4.0 5.0 EC Fig. 4.5. gapping of the P4-P5 plane into the T0-ec plane for a 0.5% solution of P-75-XH T C; a = -1, a = 12, y 2 = P3 = P5 = 0.

QP ~~~~(3,1) (212) Tc 20002 (2,0.5) 10001 2.0 3.0 4.0 5.0 CC Fig. 4.6. Mapping of the P4-P5 plane into the T0-Ec plane for a 0.9% solution of Natrosol 250-H; a -1, a = 12, P2 = P3 = Pq = 0.

139 the Tc-EC plane. Although no special attempt was made to achieve agreement between theory and experiment in any of the cases shown, it is obvious that certain appropriately chosen (but non-unique) combinations of the parameters 52 to 55 could readily accomplish this.

CHAPTER 5 SUMMARY AND CONCLUSIONS The stability of viscoelastic flow in the narrow channel between two long concentric rotating cylinders has been considered both theoretically and experimentally. In the mathematical analysis, the general Coleman-and-Noll9 simple fluid with fading memory is employed as the rheological model. Since all previous treatments of the small-gap stability problem have dealt with special cases of this model, the earlier works are shown to represent special cases of the current investigation. General techniques are presented for treating small perturbations on viscometric flows of simple fluids, and these techniques are subsequently applied in analyzing the pertinent stability behavior, under the usual assumption of neutral, axisymmetric disturbances. It is found that there are exactly eight material functions of the undisturbed rate of shear which are necessary to define the problem. Three of these are the familiar viscometric functions ir, l,, and 92, but the other five are entirely new, and their forms can be deduced, at present, only in the limit of vanishing shear rates. Furthermore, numerical calculations reveal that, in practice, only six of the functions should actually have a substantial influence on flow stability. Another important result of the theoretical investigation concerns the definition of one of the stability criteria, the Taylor number. It 140

is shown that this quantity should be calculated using the apparent fluid viscosity, as evaluated at the shearing conditions in the apparatus. Stability tests and viscosity measurements are reported for aqueous solutions of cellulose-derivative polymers whose viscometric functions have been presented in the literature, for certain concentrations. In all instances, it is found experimentally that the flow is less stable than in the Newtonian case, and that the wavelength of the disturbance is greater. By comparing these findings with the results of numerical stability calculations (based on empirical expressions for the viscometric known, approximate fluid models do not suffice to describe the observed stability behavior, and, in particular, that the "viscometric hypothesis" suggested by Metzner, White and Denn29 is unacceptable in the present instance. Additional calculations, however, reveal that the five newlydefined material functions of the present work can be appropriately (but not uniquely) chosen to give agreement between theory and experiment. At present, the measurement of the three viscometric functions ~, 04, and U2 constitutes an important part of the laboratory testing of viscoelastic fluids. This is mainly because a knowledge of these functions enables one to predict the stress distributions associated with an entire class of fluid motions, namely, viscometric flows. By way of contrast, the five new material functions of the current analysis may very possibly be peculiar only to the Taylor-Couette experiment. Therefore, since a stability test does not suffice to determine the varia

142 tion of the new material functions in a unique manner, it appears that there is but one important rheological application for this test. That is, it represents one additional means for experimentally testing simplified rheological models to determine whether they are adequate for describing more complicated flows of real viscoelastic fluids; the current theoretical development provides a basis for interpreting the results of such a test. Although the present work eliminates the need for further theoretical studies on the Taylor-Couette stability of simple fluids (at least for the small-gap, neutral-stability case), it would be interesting, as further theoretical work, to determine whether the new material functions found here are common to other stability problems.

APPENDIX A H(l)(t' ) IN CLOSED FORM In this appendix, we present the remaining details of the derivation of H(l)(t') in closed form, which was begun in Section 3.3. Let us start out by defining (t'), Mt(t'), and Nt(t') to be the summations of all terms in (3.34) involving the tensors L. M, and N respectively. It is readily shown that A2 t, tv Lt(t') = t(t') - yj) ((t2)dt2dt' +Y W) S - Lt (t2)dt2dt1 t t? t1 t2 t3 + )t(o ) /jf t /t (t4)dt4dt3dt2dtl+.. (A-l) t t t t.t(t' ) (t') 27(0 )2 1 a t t t1t t2 t3 t t.tT t2 t2 t3/ 7.+ 357(0o )4 ~J Jn J Nt(t4)dt4dt3dt2dtl.. (A-3) t t t t where A t 7 - tt )n+2 (~ ) LL(t') = (t-2t). _ Lt(t) (A-4) L+ ((n+02)/ tn 143

144 00o n Mt (t'-t) n+3 o( M(t) (A-5) (n+3) ~ )tn n=O Nt (t') N(t) (A-6) n=0 (n+4)4'ttn Substituting s = t'-t into (A-1)-(A-3), we then have Lt L (s) (t+s) - () (t+s2)ds2dsl o o + y s s /s s3 (t+s4)ds4ds3ds2dsl+... (A-7) O O O O o o N(t+s) = - t(t+S) )- 27(0O)2 (t+S2)dS2ds + 3y7(0)4 / (t+S4)ds4ds3ds2dsl+... (A-8) 0 0 0 0 o o o Next, taking the Laplace transforms of Eqs. (A-7)-(A-9) with respect to the new variable s, one finds that 1 2(0) X~ 2+(o) [4 (t+s) (A-) [(+ q 2+( A t()t+s)2 d ( A-ll) a Ixc~~~0 0 0 0L

145 If t(t+s)_1,. d_.t(t+s)] (A-12) [q2a+y(o)2]2 dS4 where q denotes Laplace's operator. These expressions may be inverted to yield convolution integrals: L(t+s) d= L t(t+sij sin y()(s-s) ds (A-13) s-t ( (t) st )dt (A-16) (t+s) = (t+sl) (t-ts)sin ((tl)d (A-17) - / t_ L tud 273 (0) N (t+s) -S $Id4 s -t(t ) = -4 2[siny~)()-3')-y(~-t' )cosy3(_)l(A- 5) or equivalently Now, by again making use of the primary rotation tensor )(t1) ~-t(t.(A-6 ), that,Mt tt) S I 3 (t'-t')sin Y (A-17) Mdt (0) Nt (t') d4 ~ _ [sinT() (t1-t' )-_T(~)(tz-t' )cosT(~)(t]-tl)ldt1 (A-18) NowI by again making use of the primary rotation tensor Rt(~)(t') (A-6) ~ that

146 dt'2 [(t')]O t )(t). l) E(o) (t' Rt, (t) dt 2 x t I' ~ -j x (A-19) d3 ]~)(t')'E(~)(t')R(~)(t')] D(1)7 (~)2 d) = t )( (t).R(~ (t )I L ) dt 3 Dt' Dt' H(o) (t ) D (l)2 (~)) +' tDt' d4 0)0(x) s)in y(o)(t) t)dt (OA-22 dt, 4 It D'Dt' (t) D(O)Ho) )(t ) D(1) (~ )) (t - 7(~ )(t-tj t (~)(-t tl1t')1t1 (A-24) Dt' Dt' t [ (t ) 1 R (O( )tl) k (). (o) (oj) yRo)(t ) o)(tl) D'(tl-t')sin y(~)(tl-t')dt (A-23) _ y(o)(tl-tI)cos y (~)(tl-t')]dt, (A-24)

APPENDIX B EXPERIMENTAL DATA The experimental data of the current investigation are reported in this section. First the results of the stability tests are presented, and then the viscosity data are given. Units are as follows: concentration in weight percent, time in seconds, angular velocity in revolutions per minute, height of vortex cells in cm., armature displacement in 0.001 in., shear rate in sec1, shear stress in dyne/cm2, and viscosity in poise. 147

148 Stability Data 1. Pure Water Temperature = 790F, revolutions = 3, time = 54.9~0.3, angular velocity = 3.28+0.02, height/20 vortex cells = 11.7~0.1, height/cell = 0.5850.00oo5, cc = 3.14+0.03, Tc = 4140 2. 7MT CMC Concentration 0.5 1.0 1.67 2.5 Temperature (OF) 78 76 80 78 Revolutions 40 100 100 100 Time 57.7~0.6 51.0+0.6 27.0+0.3 17.2+0.5 Angular velocity 41.6+0.4 117.5+1.1 222+3 348~10 Shear rate 52.3+0.5 148+2 279+3 438+13 Critical viscosity 0.111 0.33 o.64 1.41 Tc 4140 3790 3630 1880 Height/20 cells 12.0+0.2 12.6+0.2 13.0~0.1 13.7~0.1 Height/cell 0.600+0.010 0.630+0.010 o.650+0o~o.oo5 0.685+~o0.005 cc 3.04+0.05 2.91+0.05 2.82~0.03 2.68+0.02 3. P-75-XH CMC Concentration 0.183 0.367 0.55 Temperature (OF) 76 76 76 Revolutions 100 100 100 Time 48.5+1.0 27.8+0.5 20.5+0.3 Angular velocity 123.8+2.6 216~4 293+4 Shear rate 155.5+3.2 271+5 368+5 Critical viscosity 0.327 0.590 0.880 Tc 4200 3950 3270 Height/20 cells 13.3+0.5 15.5~0.5 16.5+0.5 Height/cell 0.665+0.025 0.775+0.025 0.825~0.025 Cc 2.76+0.10 2.36+0.08 2.22+0.07 4. Natrosol 250-H HEC Concentration o.48 o.80 1.00 Temperature (OF) 76 80 82 Revolutions 100 200 200 Time 36.8+~o.5 39.2~0.5 30.9~0.5 Angular velocity 163+2 306+4 388+6 Shear rate 205+3 384+5 488+8 Critical viscosity 0.470 0.950 1.32 Tc 3550 3080 2590 Height/20 cells 14.3+0.5 16.3+0.5 17.0+0.5 Height/cell 0.715+0.025 0.815+0.025 0.850+0.025 Cc 2.56+o0.09 2.25+0.07 2.16+0.06

149 Visco sity Data 1. Pure Water Viscosity of water at 79OF32 = 0.00872 2. 7MT CMC 0.5% Solution (25.6oC) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.5 0.018 52.8 5.59 0.104 1.4 0.024 66.4 7.45 0.112 1.3 0.029 83.6 9.00 0.108 1.2 0.038 105.3 11.8 0.112 1.1 o.o46 132.9 14.3 0.107 1.0 0.059 167.0 18.2 0.109 0.9 0.074 210.0 23.0 0.110 0.8 o.og090 264.0 27.9 0.106 0.7 0.112 334.o 34.8 0.104 o.6 0.138 420.0 42.8 0.102 0.5 0.172 528.0 53.4 0.101 0.4 0.211 664.o 65.5 o.986 0.3 0.256 836.o 79.5 0.950 0.2 0.314 1053.0 97.5 0.929 0.1 0.387 1329.0 120.0 0.903 0.0 o.480 1670.0 149.0 0.892 1.0% Solution (24.5OC) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.5 0.061 52.8 18.9 0.358 1.4 0.075 66.4 23.2 0.350 1.3 0.093 83.6 28.8 o.344 1.2 0.117 105.3 36.3 o.346..1 0.144 132.9 44.7 0.336 1.0 0.176 167.0 55.6 0.333 0.9 0.214 210.0 66.5 0.316 o.8 0.260 264.0 80.6 0.306 0.7 0.318 334.0 98.6 0.295 o.6 0.378 420.0 117.0 0.285 o.5 0.462 528.0 143.0 0.271 0,4 0.560 664.0 174.0 0.262 0.3 0.670 836.o0 208.0 0.249 0.2 0.800 1053.0 248.0 0.236 0.1 o.g30 1329.0 288.0 0.216 0.0 1,14 1670.0 354.0 0.212

150 1.67% Solution (26.7oC) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.5 0.143 52.8 44.4 0.840 1.4 o.176 66.4 54.6 0.823 1.3 0.212 83.6 65.8 0.786 1.2 0.258 105.3 80.0 0.762 1.1 0.312 132.9 96.8 0.728 1.0 0.377 167.0 117.0 0.700 0.9 0.455 210.0 141.0 0.672 0.8 0.540 264.0 168.0 o.636 0.7 0.655 334.0 203.0 o.608 0.6 0.775 420.0 240.0 0.571 o.5 0.910 528.0 282.0 0.535 0.4 1.08 664.o 335.0 o.505 0.3 1.28 836.o 397.0 0.475 0.2 1.52 1053.0 471.0 o.448 0.1 1.77 1329.0 549.0 0.413 0.0 2.05 1670.0 636.o 0.381 2.5% Solution (25.60C) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.5 0.435 52.8 135 2.56 1.4 0.535 66.4 166 2.52 1.3 o.630 83.6 195 2.34 1.2 0.735 105.3 228 2.17 1.1 o.885 132.9 275 2.07 1.0 1.04 167.0 322 1.93 o. 9 1.23 210.0 382 1.82 0.8 1.44 264.0 447 1.69 0.7 1.68 334.0 521 1.56 o.6 1.93 420.0 599 1.46 0. 5 2.22 528.0 689 1.31 0.4 2.52 664.o 781 1.18 0.3 2.90 836.o 90o 1.08 0.2 3.40 1053.0 1055 1.00 0. 1 3.95 1329.0 1225 0.920 0.0 4.45 1670.0 1370 0.820

151 3. P-75-XH, CMC 0.183% Solution (24.5~C) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.5 0.077 52.8 23.9 0.453 1.4 0.090 66.4 28.9 0.435 1.3 o.106 83.6 32.8 0.392 1.2 0,126 105.3 39.1 0.372 1.1 0.147 132.9 45.6 0.343 1.0 0.171 167.0 53.0 0.318 0.9 0.199 210.0 61.7 0.294 0.8 0.229 264.0 71.0 0.269 0.7 0.260 334.0 80.6 0.242 o.6 0.300 420.0 93.0 0.222 0.5 0.342 528.0 o106.o 0.201 0.4 0.390 664.0 121.0 o.183 0.3 o.448 836.o 139.0 0,166 0,.2 0.515 1053.0 160o. 0.153 0.1 0.585 1329.0 182.0 0.137 0.0 o.678 1670.0 210.0 0.126 o0.367 Solution (24.5~C) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.5 0.194 52.8 60.2 1.14 1.4 0.225 66.4 69.8 1.05 1.3 0.260 83.6 80.6 0.964 1.2 0.298 105.3 92.4 0.880 1.1 0.342 132.9 106.1 0.797 1.0 0.387 167.0 120.0 0.718 0.9 0.445 210.0 138.1 o.657 o.8 0.510 264.0 158.1 0.599 0.7 o.580 334.0 180.1 0.539 o.6 o.658 420.0 204.0 o.486, 0.5 0.740 528.0 230.0 0.436 0.4 0.825 664.0o 256.0 0.386 0.3 0.925 836.0o 287.0 0. 44 0.2 1.03 1053.0 320.0 0.305 0.1 1.16 1329.0 36o0.0 0.270 0.0 1.31 1670.0 406.0 0.243

152 0.55~ Solution (24.5~C) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.5 0.375 52.8 116.3 2.20 1.4 0.425 66.4 132.0 1.99 1.3 o.490 83.6 152.0 1.82 1.2 0. 558 105,3 173.0 1.65 1.1 o.630 132.9 195.5 1.47 1.o 0.710 167.0 220.0 1.32 0.9 0.790 210.0 245.0 1.17 0.8 0.885 264.0 274.0 1.04 0.7 O.982 334.0 305.0 0.913 o.6 1.090 420.0 338.0 o.80o O.5 1.22 528.0 378.0 0.716 0.4 1.36 664.o 422.0 0.637 0.3 1,50 836.o 46,.o o.556 0.2 1.67 1053.0 518.0 0.493 0.1 1.84 1329.0 571.0 0.429 o.o0 2.02 1670.0 626.0 o.375 4. Natrosol 250-H HEC 0.48* Solution (24.5OC) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.5 0.132 52.8 40.9 0.776 1.4 o.156 66.4 48.3 0.729 1.3 0.181 83.6 56.1 0.671 1.2 0.208 105.3 64.5 0.614 1.1 0.239 132.9 74.1 0.557 1.0 0.275 167.0 85.3 0.511 0.9 0.315 210.0 97.6 o.46s 0.8 0.358 264.0 111.0 0.420 0.7 0.402 334.0 125.0 0.374 0.6 o.460 420.0 143.0 o.340 o.5 0.520 528.0 161.0 0.305 0.4 o.590 664.o 183.0 0.276 0.3 0.670 836.0 208.0 0.249 0.2 0.755 1053.0 234.0 0.223 0.1 0.840 1329.0 260.0 o.196 o.o 0o.g955 1670.0 296.0 0.177

153 o.8% Solution (26.70C) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.5 0.465 52.8 144 2.73 1.4 0.538 66.4 167 2.52 1.3 0.600 83.6 186 2.22 1.2 0.675 105.3 210 2.00 1.1 0.745 132.9 231 1.74 1.0 0.822 167.0 255 1.53 o.9 o.91o 210.0 282 1.34 0.8 1.00 264.0 310 1.18 0.7 1.11 334.0 344 1.03 0.6 1.21 420.0 376 0.896 0.5 1.34 528.0 416 0.787 0.4 1.47 664.o 456 o.687 0.3 1.61 836.o 500 o.598 0.2 1.76 1053.0 546 0.520 0.1 1.92 1329.0 595 0.447 0.0 2.10 1670.0 651 0.390 1.0% Solution (27.8oC) Gearbox Armature Displacement Shear Rate Shear Stress Viscosity 1.0 1.44 167 447 2.68 o.9 1.56 210 484 2.30 0.8 1.69 264 525 1.99 0.7 1.83 334 568 1.70 0.6 1.96 420 608 1.45 o.5 2.12 528 658 1.25 0.4 2.31 664 717 1.08 0.3 2.51 836 779 0.931 0.2 2.70 1053 838 0.797 0.1 2.95 1329 915 o.688 0.0 3.18 1670 986 o.590

APPENDIX C CRITICAL EIGENFUNCTIONS The critical eigenfunctions *c and Vc are defined herein as the eigenfunctions associated with the critical parameters Tc and Ec. For any particular simple fluid, *c and Vc can be determined using the present numerical techniques, by integrating the disturbance equations from X = 0 to X = 1 for T = Tc and c = cc. As an example, Figs. C.1 and C.2 show plots of Vc and Vc vs. X corresponding to the Newtonian case with a = -1. 154

.03l.02 \J1.01 ~0 0.1 0.2 0.3 04 05 06 07 0.8 0.9 1.0 Fig. C.1. Critical eigenfunction *c vs. X for Newtonian case; a = -1.

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 X Fig. C.2. Critical eigenfunction Vc vs. X for Newtonian case; a = -1.

APPENDIX D COMPUTER ANALYSIS A main program and three subroutines (BETA., SOLVE., and INVERT.) were used in this work to (a) numerically solve the characteristic value problem arising as a consequence of the current Taylor-Couette stability analysis, and (b) determine the critical parameters Tc and Ec. The basic techniques employed have been presented in Section 4.7.3. All programs were written in the "MAD" computer language, and the time required to determine each eigenvalue was approximately 2 sec. The following is a summary of the primary duties of the main program: (1) To direct the reading of data and the prinout of results. (2) To initiate and implement the "golden section" minimum searching procedure. (3) To call upon the subroutine BETA. to calculate the rheological coefficients Po through P, for any prespecified dependence of these quantities on the dimensionless shear rate 7T. (4) To call upon the subroutine SOLVE. to integrate the disturbance equations from X = 0 to X = 1, using a Runge-Kutta technique.* (5) To determine the matrix of partial derivatives appearing in Eq. (4.92). (6) To call upon the subroutine INVERT.** to solve Eq. (4.92). *The computer library Runge-Kutta subroutine was employed in SOLVE. **INVERT. is a modified version of the program appearing on pp. 390391 of Applied Numerical Methods.3 157

158 (7) To iterate on the quantities G1, G2, and G3 until satisfactory convergence is achieved. A description of the important symbols employed is given below. Program Symbol Definition A a = rl/d A(I,J) Elements of the matrix of partial derivatives appearing in Eq. (4.92). A vector consisting of V, Di, and V at X = 1 is appended as an additional column. ALPHA a = (n2-n1)/Ql B aa BET(I) Pi (i = 0,...,9) COEF(I) Ci (i = 1,2,) in Eqs. (4.113) CRIT Cc DEL(I) 6Gi (i = 1,2,3) in Eq. (4.92) DELX AX, step size in Runge-Kutta integration EL El = elasticity number = Tro/pd2 EPS C GESNU(I) GUESS(I) * RATIO(I) GUESS(I) Gi (i = 1,2,3) LEFT Left boundary in the golden section interval (LEFT < cc < RIGHT) MAX Total number of applications of the golden section test NU Maximum allowable fractional change in Gi during final iteration

159 Program Symbol Definition POWER(I) Pi (i = 0,1,2) in Eqs. (4.113) PSI(I) *i (i = 1,...,;..,6); 1, = DS, *3 = V, *4 = D2*, *5 = D3~, 46 = DV RATIO(I) (1.0-RATIO(Ii)) represents the fractional change in Gi used in numerically determining the nine partial derivatives in Eq. (4.92) RIGHT Right boundary of the golden section interval (LEFT < Ec < RIGHT) R Modified Reynolds number, Rm SOLN Taylor number (eigenvalue) TAUGAM |7T| TAYLOR Tc TEST Boolean variable which specifies whether or not the critical eigenfunctions tc and Vc should be printed The MAD program listings and a typical computer output are presented in the pages which follow.

$COMPILE MAD,PRINT OBJECT,DUMP,EXECUTE 036933 03/24/671 l:37 -3:.1 A; MAD (06 JAN 1967 VERSION) PROGRAM LISTING. THIS PROGRAM SOLVES THE EIGENVALUE PROBLEM WHICH ARISES AS A RESULT OF ANALYSIS OF THE STABILITY OF A VISCOELASTIC FLUID CONTAINED BETWEEN CONCENTRIC ROTATING CYLINDERS. THE EIGENVALUE PROBLEM IS HEREIN BROKEN UP INTO A SET OF 6 FIRST ORDER DIFFERENTIAL EQUATIONS IN 6 UNKNOWNS (PSI(1)...PSI(6)) SUBJECT TO THE BOUNDARY CONDITIONS (AT X=O AND X=1), PSI(1)...PSI(3)=O. THE COEFFICIENTS IN THE DIFFERENTIAL EQUATIONS ARE ALL FUNCTIONS OF THE TAYLOR NUMBER (GUESS(3)), THE DIMENSIONLESS WAVE NUMBER IEPS), AND GEOMETRIC AND RHEOLOGICAL PARAMETERS. SINCE, IN AN EIGENVALUE PROBLEM, THE SOLUTION IS SPECIFIED UP TO AN ARBITRARY CONSTANT MULTIPLIER, PSI(4) IS TAKEN TO BE EQUAL TO 1.0 AT X=O. FOR ANY ASSUMED VALUES TO THE TAYLOR NUMBER AND MISSING INITIAL CONDITIONS (AT X=Ot PSI(5)=GUESS(1) AND PSI(6)=GUESS(2)), THE EQUATIONS ARE INTEGRATED BETWEEN X=O AND X=1 USING A RUNGE-KUTTA FOUR-POINT INTEGRATION SCHEME. OF COURSE, FOR AN ARBITRATY GUESS, THE VALUES OF PSI(1)...PSI(3) AT X=1 WOULD NOT LIKELY TURN OUT TO BE ZERO, BUT THE INITIAL GUESS IS SUBSEQUENTLY CORRECTED USING A NEWTON-RAPHSON TECHNIQUE. DERIVITIVES FOR THIS TECHNIQUE ARE HEREIN ESTIMATED BY A FINITE DIFFERENCE APPROXIMATION. ITERATION CONTINUES UNTIL THE MAGNITUDES OF THE CORRECTIONS TO THE GUESSES BECOME SMALL COMPARED TO THE GUESSES THEMSELVES. THE CRITICAL TAYLOR NUMBER (I.E. THE MINIMUM TAYLOR NUMBER ON OVER ALL DISTURBANCE WAVE NUMBERS) IS LOCATED USING THE GOLDEN SECTION SEARCHING TECHNIQUE. THE PROGRAM EXAMINES THE REGION BETWEEN EPS=LEFT AND EPS=RIGHT FOR A MINIMUM AS FOLLOWS. AFTER EVALUATING THE TAYLOR NUMBER AT THE.382 AND.618 POINTS OF THIS INTERVAL, IT ELIMINATES ONE PORTION (EITHER BETWEEN 0 AND.382 OR BETWEEN.618 AND 1) AS NOT CONTAINING THE MINIMUM. LEFT UOR RIGHT IS THEN ADJUSTED SO THAT THE INTERVAL IS SHRUNK TO.618 OF ITS ORIGINAL SIZE. THIS PROCEDURE IS REPEATED COUNT TIMES..... DIMENSIONING AND MODE DECLARATION.... DIMENSION PSI(6),GUESS(4),GESNU 4),DEL(3),RATIO(4), A4*4), *f701 1 BET(9),EPS(2),SOLN(2),CUEF(9),POWER(9)}.01 INTEGER I,N,INDEXCOUNT,K,MAX *0.02 BOOLEAN TEST *003 REFERENCES ON *004.... PRINTOUT.... START PRINT COMMENT $1SOLUTION PERFORMED FOR FOLLOWING VALUES OF PA *005 1 RAMETERSS *005 READ AND PRINT DATA *006 B=ALPHA*A *C07 PRINT COMMENT SOTHE TAYLOR NUMBER AS A FUNCTION OF THE DIMENS *008 1 IONLESS WAVE NUMBER IS GIVEN BY$ *O-08

...GOLDEN SECTION.... INDEX=O *009 COUNT=O *0 10 EPS ( 1 )=LEFT+. 382* (RIGHT-LEFT) *011 EPSt2)=LEFT+.618*(R IGHT-LEFT) *0,12 GOLDEN COUNT=COUNT+l *013 THROUGH END, FOR VALUES OF K=1,2 *014 WHENEVER K.E.INOEX, TRANSFER TO END *0 15 1... INITIALIZATION.... BEGIN THROUGH NEXT, FOR I=l,l,I.G.4 *016 Cl GESNU( I )=GUESS( I ) 0117 2 NEXT A(I,4)=0. * 018 THROUGH SETUP, FOR VALUES OF N=4,,t2,3 019 "l GESNU (N)=GUESS(N) *RAT IO(N) *020 2 WHENEVER N. GE.3 *021 2 GESNU(2 )=GUESSt2) *022 l 2 R=SQRT.(GESNU(3)/(-2.*B)) *n23 1 3 EXECUTE BETA.(BET,R,BEL,COEF,POWER) *024 I -i 2 END OF CONDITIONAL *025 I WHENEVER N.E.2, GESNU(Il=GUESS(l) *C26 ^2 THROUGH INIT, FOR I=Il,I.G.3 *027 C2 INIT PSi I(I )=O0. *028 3 PSI (4)=1. *029 2 PSI (5)=GESNU( 1) *030'2 PSI (6)=GESNU(2) *031.... INTEGRATION.... EXECUTE SOLVE.(EPS(K),ALPHA,BETGESNU(3),R,PSI,DELX,O,A) *032 n2.... SETUP OF MATRIX OF COEFFICIENTS.... THROUGH SETUP, FOR I=1,1,I.G.3 *033 r2 SETUP A(I,N)=(PSI(I)+A(I1,4))/ (GUESS(N)*(RATIO(N)-1.)) *034:... MATRIX INVERSION.... EXECUTE INVERT. (A,DEL,3) *035.... ITERATION ON GUESSES.... THROUGH CHANGE, FOR N=1I,,N.G.3 *036 91 CHANGE GUESSINI:GUESS(N)+DEL(N) *037 02 THROUGH VEST, FOR I=1:1,I.G.3 *038,1 VEST WHENEVER.ABS.(DEL(IfJ.G..ABS.(NU*GUESS(I)),TRANSFER TO BEGIN *039.... GOLDEN SECTION.... R=SQRT.(GUESS(3)/1(-2.*B)) *040 1

EXECUTE BETA. (BETRtBoEL,COEFPOWER) *041 01 END ~SOLN(K) =GUESS (3)*BET(O).P.2 *042 01 END ~~PRINT RESULTS EPS(K),# SOLN(K) *043 01....CHECK ON WHETHER INTERVAL HAS BEEN DIVIDED SUFFICIENTLY.. WHENEVER COUNT.E.MAXt TRANSFER TO SKIP *044...INTERVAL.SHRINKI.NG.... WHENEVER S0LN(2).L.SDLN(l) 04 LEFT=EPS(l1) *04650 EPStl)=EPS(2) *046 01 SOINLI1)=SOLN( 2) *048 01 EPS (2)=LEFT+.6I8* (RIGHT-LEFT) *049 01 INDEX1l *050 01 OTHERWI SE *051 01 RIGHT=EPS (2) *052 01 EPS(2)=EPS( 1) *053 01 SOLN(2)=SOLN(l1) *054 Cl EPS( 1)=LEfT+.38Z* (RIGHT-LEFT) *055 01 INDEX=2 *056 01 END OF CONDITIONAL *057 01 TRANSFER TO GOLDEN *058...FINAL RESULTS..... SKIP WHENEVER SDLN(2).L.SDLN(1) *059H CRIT=EPS( 2) *060 01 TAYLOR=SOLN( 2) *061 01 OTHERWI SE *062 Cl CRIT=EPS( 1) *063 01 TAYLOR=SOLN( 1) *064 01 END OF CONDITIONAL *065 01 PRINT COMMENT $OTHE CRITICAL TAYLOR NUMBER IS $ *066 PRINT RESULTS TAYLOR *067 PRINT COMMENTSOFOR A CRITICAL DIMENSIONLESS WAVE NUMBER OF$ *068 PRINT RESULTS CRIT *069 PRINT COMMENT SOTHE RHEOLOGICAL PARAMETERS ARE$ *070 PRINT RESULTS BET(O)...BET(9) *071 WHENEVER TEST *072 THROUGH FINAL, FOR I=1,1,I.G.3 *073 01 FINAL Psi (I).=.0, *074 01 01 PSI (4)=1. *075 01 PSI (5)=GUESS( 1) *076 01 PSI(6)=GUES $2) *077 01 WHENEVER INAEX.E.I1 *078 01 I NDE X=2 *079 02 OTHERWISE *080 02 INDEXJl *081 02 END OF CONDITIONAL *082 02 PRINT COMMENTSOTHE CRITICAL EIGENFUNCTIONS ARES *083 01 EXECUTE SOLVE.AEPS(.INDEX),ALPHA,BET,GUESS(3).RPS1,DELXl,A) *084 01 END OF CONDITIONAL *085 01 PRINT RESULTS GUESS(1)...GUESS(3) *086 TRANSFER TO START *087 END OF PROGRAM *088

TABLE OF REFERENCES ***SYMBO0L TABLE REFERENCES ALPHA 007,032,084 A 007,018,032,034,035,084 BEGIN 016,039 BET 024,032,041,042,071,084 B 007,0239024,0409041 CHANGE 036,031 C0EF 024,041 COUNT 010,013i,044 CR11 060,063,069 DEL 035,037,039 DELX 032,084 EL 024,041 END 014,015,043 EPS 011,012,0O32,043,046,047,049,05290539055,060,063,084 FINAL 073,074 GESNU 017,020,022,023,026,030,031,0O32 GOLDEN 013,058 GUESS 017,020,022,026,034,037,039,040,C42,076,077,084,086 INDEX 009,015,050,056,078,079,081,084 INIT 027,028 I 016,017,018,P027,028,033,034,038,039,073,074 K 014,015,032,042,043 LEFT 011,012,046,049,055 MAX 044 NEXT 016,018 N 019,020,021,026,034,036,037 NU 039 POWER 024,041 PSI 028,029,030,031,032,034,0749,075,076,077,084 RATIO 020,034 RIGHT 011,012,049,052,055 R 023,024,032,040,041,084 SETUP 019,033,034 SKIP 044,059 SOLN 042,043,045,048,054,05.9,06l,064 START 005,087 TAYLOR 0619,064,067 TEST 072 VEST 038,039 ***FUNCTION REFERENCES BETA 024,041 INVERT 035 SOLVE 032,084'SQRT 023,040 THE FOLLOWING NAMES HAVE OCCURRED ONLY ONCE IN THIS PROGRAM. COMPILATION WILL CONTINUE. MAX *044 NU *039 TEST *072

SCOMPILE MADPRINT OBJECTDUMPtEXECUTE 036933 03/24/67 10:37 48.9 AM MAD (06 JAN 1967 VERSION) PROGRAM LISTING........ EXTERNAL FUNCTION(EPS ALPHABETTRPSIDELXBOOLA) *001 ENTRY TO SOLVE. *002 THE SUBROUTINE SOLVE. INTEGRATES THE DIFFERENTIAL EQUATIONS ASSOCIATED WITH THE EIGENVALUE PROBLEM UNDER CONSIDERATION (IN THE INTERVAL O.L.X.L.1) SUBJECT TO ANY PRESCRIBED INITIAL CONDITIONS AND TAYLOR NUMBER. A RUNGE-KUTTA FOUR POINT INTEGRATION SCHEME IS EMPLCYED.....DIMENSIONING.... DIMENSION Q(6),F(6) *003 BOOLEAN BOOL *004....INITIALIZ ATION.... EPSQ=EPS*EPS *005 X=O. *006 EXECUTE SETRKD. (6PSI(1).F(11,Q,X,OELX) *007 WHENEVER BOOL *008 PRINT COMMENTSO X *009 01 I 1 PSI PSI' VS *009 END OF CONDITIONAL *O10 e1 WRITE WHENEVER BOOL, PRINT FORMAT RESULT,X,PSI(1).. PSI(3 *011....CHECK ON WHETHER INTEGRATION IS COMPLETE.... WHENEVER X.G..99, TRANSFER TO FINISH *012....INTEGRATION.... CALC S=RKDEQ.(O) *013 WHENEVER S.L.1.5 *014 F(1)=PSI(2) *015 01 F(2 )=PSI(4) *016 0 1 F(3)=PSI(6) *017 01 F(4)=PSI(5) *018 01 F(6)=(EPSQ*PSI(31+T*BET(O)*PSI(1)/(2.*R)+BET(2)*PSI(4) *019 01 1 -BET(3 *EPSQ*PSI(1) 1/BET( 1) *019 F(S)=EPSQ*(BET(4)*PSI(4)-EPSQ*BET(5*PSI (1) *020 01 1 -2.*R*( 1.+ALPHA*X)*PSI(3)*BET(O)+BET(6)*F(6) *020 2 +BET(7)*EPSQ*PSI(3)+(BET()81*PSI(6)-2.*BET(9)*PSI(2)/A) *020 TRANSFER TO CALC *021 01 END OF CONDITIONAL *022 01 TRANSFER TO WRITE *023 VECTOR VALUES RESULT=S1H,F36.2,3FL8.6*$ *024 FINISH FUNCTION RETURN *025 END OF FUNCTION *026

$COMPILE MAOPRINT OBJECTOUMPtEXECUTE 036933 03/24/67 10:37'5.7 3M MAD (06 JAN 1967 VERSION) PROGRAM LISTING EXTERNAL FUNCTION(BET,RBELPCOEFPOWER) *(001 ENTRY TO BETA. *002 THE SUBROUTINE BETA. CALCULATES THE RHEULOGICAL PARAMETERS BET(O)...BET(9) WHICH APPEAR IN THE ELGENVALUE PROBLEM UNDER CUNSIDERATION..... FREQUENTLY APPEARING QUANTITIES.... TAUGAM=-B*R *EL *003 SIGl=-COEF(1)*TAUGAM.P.(POWER(O)-POWER(l)) *004 SIG2=-COEF(2)*TAUGAM.P.(PCWER(O)-POWER(2)) *005....CALCULATION OF BET(O)...BET(9).... BET(O)=TAUGAM.P.POWER(O) *006 BET(1)=1.-POWER(O) *007 BET12)(SIGI+SIG2 )/2. *008 BET(3)=(SIGI+S1G2)/2. *009 BET(4)=2. *010 BET(5 )=1. *0I1 BET(6J=POWER(2)*51G2 *012 BET(7)=SIG2 *013 13 BET(8)=(1.-POWER( 1) )*SIGL *014 BET(9)0=. *015 FUNCTION RETURN *016 END OF FUNCTION *017

SCOMPILE MAD,PRINT OBJECT.DUMP,EXECUTE 999107 06/22/66 10 18 36.0 AM M-AD (29 APR 1966 VERSION).PROGRAM LISTING... EXTERKN-AL -FUNCT[QN ( -A,-ECtM) *00 1 ENTRY TO INVERT. *002 THE SUBROUTINE- I'-N-VE-R-T'. IS AN ADAPTATION OF THE PROGRAM APPEARING ON PP.390-391 OF'APPLIED NUMERICAL METHODS' BY CARNAH-AN-,;iW[TKES,'AND LUTHER. IT SOLVES M SIMULTANEOUS LINEAR EQUATIONS USING THE GAUSS-JORDAN REDUCTION SCHEME WITH THE — MAXIMUM PIVOT- CRITERION. THE MXM MATRIX OF COEFFICIENTS WITH THE VECTOR OF CONSTANTS APPENDED AS THE M+1 TH COLUMN IS STORED IN THE A ARRAY. NORMAL MODE IS INTEGER 00C3 FLOATING POINT BIGA, AJCK, A, DEL *004 DIMENSION i Ri1''10), JC( 10) *C05 MP =M+1 *006 THROUGH Li, FOR K=i1,1K.G.M *0C7 BIGA=O. *008 01 THROUGH L2, FOR I=1,i,I.G.M *009 01 THROUGH L2, FOK J=1,i,J.G.M *010 02 THRO-UGH L3, FORK I'1-i,, I1.E.K *011 03 THROUGH L3, FOR JI=1,1,J1.E.K *012 04 L3 WHENEVER I.E.IR(I).OR.J.E.JC(J1), TRANSFER TO L2 *013 05 WHENEVER.ABS.A(I,J).G. BIGA *014 03 -IGA=.ABS.A(I,J) *.015 01 03 0n IRIK)=.I *016 01 03 -JC (K) -'J *017 01 03 L2 END OF CONDITIONAL *018 01 03 BIGA=A(IR'K),JC(K)) *019 01 THROUGH L4, FOR J=1,.1J.G.MPI *020 01 L4 A(IR(K),J)=A( IR(I K), J)/BIGA *021 02 THROUGH Lit, FOR I=1,itI.G.M *022 01 WHENEVER'-' I.NSER'('K'{' *0,23 02 AJCK=A( I, JC (K) *024 01 02 THROUGH L6, FOR J=l,1,J.G.MPI *025, 01 02 L6 A(I,J)=A(I,J)-AJCK*A(IR(K),J) *026 01 03 L1 END OF CONDITIONAL *027 01 02 THROUGH L7, FOR I=1,1,I.G.M *028 L7 DEE 1 J C ( I -'R=AF- IRI MP1- *029 01 FUNCTION RETURN *030 END OF FUNCTION *031

SOLUTION PERFORMED FOR FOLLOWING VALUES OF PAKAVETERS NU=lE-3, ELX=.05,,ATI.0(1)=1.) 05,1.O1l,.01,2.,MAX=13,TEST=OB, ALPHA=-l.OA=12.O,, IGHT=4t.U,LEFT=3.0,COEF( 1= 0.00,0.00 BET(2)=C.0,BET(3)=U.0,BET(4)=2.,BET(5)=l.,B3ET(9)=0. POWER(0)=0.O,.423,.375,EL=5.90,GUESS( 1)=-7.0*,-.5, 3E3,-1.,TEST=1B-* IME tAOR1 UZ46UMBER AS. A_-UNCTION OF -THE Ol1M.ENSl0NLSS_ WAVE NUMBER _.1S___GIVENLBY EPS(I.) = 3.3820009 SOLN(1) = 3420.783661 EPS(2) - 3.618000, SGLN(2) - 3499.300537 -ESZ 1 =-.3.236076, SCLN(L) & 3395.898743_ EPS(2) = 3.382000, SOLN(2) 3420.783661. EPS(1) = 3.145924, SULNIL) = 3390.279144 -,EP.4-2) 2- 3 60 76,.-.. SLN(2 = 12 339 5..98743_. EPS(1) 3.090181, SCLN(l) = 3390.762146 EPS(2) = 3.145924, SGLN(2) = 3390.279144, -~P-4 — z -.. —- 3.i145924-t -SOLN(1) 3 3390-.279.144 -EPS(2) = 3.180344, SOLN(2) 3391.511444 EPS(l) = —.3.124623, SULNIL) = 3390.099457 r-P "-( 21 = n,~L ~ - SOLN12- =. 3390.229.144 EPSLJ = 3.111475, SOLN(1). 3390 212.036 _.EPRS(I2 2= _Z.12-462-3 SULN(2) = -3390.099457 IP S( 1 - 3 1246 23 -. — SQ-L9tiJ.. 3390-.099457 __EPS 1 s - -- 3.1327,64 SOLN(21 = 339-0.116333 -IfPSJIL4 — ___ ~- 3,~~.4139607 SOLNf11) = 3390.122253 EPS(2) = 7-3 124623 —------ -SC4.N-L2)-) — - ~310 —0 — n9-57~ L-p&SI SLI3f"24623,._. SOLNL1J, 339l.n0945_FP P L 3.127738, SOLNt2J —= - - -3a9-0O-.0a8358-&P-S12).= ~..? 3.~129654, 3.17_18 ~SQLN 1.. 339f.0- 1023 5 8-F-P-SA 2):-. —_12t 96 54 SC L N z~ 339" 1 02164

EPS.(1)..........3,12 65 45.t. SOLN(l.1 = 3390.097748 EP SL2 __ 3.12 738. SLN(2.__. 3390.098358. _EPSl -1). b.125813.SOLN ( 1 ) = 3390.098083 EPS(2) = 3.126545, SCLN(2) = 3390.097748 fHE CRITICAL TAYLQR_.NU4ER IS_ —._ TAYLOR = 3390.097748 FOR ACRITICAL DIMENSIDNLESS WAVE NUMBER OF _CRiT. = 3. 12_6545.... Jf EJHE_0L0GlCAL PARAME IEBS ARE BET (0).. BET( 9)' 1.OOOOOOE+00 1.000000E+oo.OOOOOOE+00 OOOO00 0E+00 2.OOOOOOE+00 1.OOOO00000E+00 -.OOOOOOE+00 -.OOOOOOE+00 -. OOOOOE+00.OOOOOOE+00 iHE CRITICAL EIGENFUNCTIONS ARE x _P.....PS1._ PSI'. V....00.000000.000000.000000.05.001109.041665 -.024597.10.003915.068225 -.049351.15.007718.081941 -.073886 H.20.011931.084979 -.097460. 25- *,016073 _.079422. -. 119160.30.0 19764.067270 -.138048 _, 35.022722.050421 -.153277.40.024757.030657 -. 164167'.45_.025'766.009610 -.170261.50.025722 -.011255 -.171349. _ 5_, 0.2?4666 _-.03_0655 -. 167471 ~.60.02Z699 -.047497 -.158908.65.0199'74 -.060866 -. 146146.70.Q 16683 -.070002 -.129843.75.013055 -.074259 -.110768. 80.009348 - 073054 -.089742.85., 0058 51 _. _6.5796 -_0 6..90.002881 -.051810 -.044939 __.95. 000795. _ -.3 0251_' -. 0223 57 1.00.000000.000000.000000 GUE$S. ( L,.LG E SS( 3.3 -6. 39tE+Q4 -4,90693E-0 1 33900g98E+0 3.

BIBLIOGRAPHY 1. Beard, D. W., Davies, M. H., and Walters, K., "The Stability of Elastico-Viscous Flow Between Rotating Cylinders, Part 3, Overstability in Viscous and Maxwell Fluids," J. Fluid Mech., 24, 321 (1966). 2. Bird, R. B., Stewart, W. E., and Lightfoot, E. N., Transport Phenomena, John Wiley and Sons, Inc., New York (1960). 3. Carnahan, B., Luther, H. A., and Wilkes, J. 0., Applied Numerical Methods (preliminary edition), John Wiley and Sons, Inc., New York (1964). 4. Carter, L. F., A Study of the Rheology of Suspensions of Rod-Shaped Particles in a Navier-Stokes Liquid, Ph.D. Thesis, University of Michigan, Ann Arbor, Mich. (1967). 5. Chandrasekhar, S., "The Stability of Viscous Flow Between Rotating Cylinders, " Mathematica, 1, 5 (1954). 6. Chan Man Fong, C. F., "On the Stability of Elastic-Viscous Flow Between Rotating Cylinders," Rheologica Acta, 4, 37 (1965). 7. Coleman, B. D., Markovitz, H., and Noll, W., Viscometric Flows of Non-Newtonian Fluids, Springer-Verlag, New York (1960). 8. Coleman, B. D. and Noll, W., "An Approximation Theorem for Functionals with Applications in Continuum Mechanics," Arch. Rational Mech. Anal., 6, 355 (1960). 9. Coleman, B. D. and Noll, W., "Simple Fluids with Fading Memory," in Second-Order Effects in Elasticity, Plasticity, and Fluid D2namics, Reiner, M., and Abir, D., Eds., Pergamon Press, Inc., New York, 530 (1964). 10. Crandall, S. H., Engineering Analysis, McGraw-Hill Book Co. Inc., New York (1960). 11. Datta, S. K., "Note on the Stability of an Elasticoviscous Liquid in Couette Flow," Physics of Fluids, Z, 1915 (1964). 12. Davies, M. H. "A Note on the Stability of Elastico-Viscous Liquids," Appl. Sci. Res., A15, 253 (1965). 169

170 13. Einstein, A., "Eine neue Bestimmung der Molektuldimensionen," Ann. d. Physik, 1_, 289 (1906). 14. Ericksen, J. L., "Anisotropic Fluids," Arch. Rat'l. Mech. Anal., 4, 231 (1960). 15. Ericksen, J. L., "Poiseuille Flow of Certain Anisotropic Fluids," Arch. Rat'l. Mech. Anal., 8, 1 (1961). 16. Farol Research Engineers Ltd., The Weissenberg Rheogoniometer Instruction Manuel. 17. Fredrickson, A. G., Principles and Applications of Rheology, Prentice-Hall, Inc., Englewood Cliffs, N. J. (1964). 18. Giesekus, H., "Zur Stabilitat von Stromungen viskoelastischer Flussigkeiten," Rheologica Acta, 5, 239 (1966). 19. Goddard, J. D., presented at the Society of Petroleum Engineers Symposium on Mechanics of Rheologically Complex Fluids, Houston, Texas, Dec. (1966). 20. Goddard, J. D. and Miller, C., "An Inverse for the Jaumann Derivative and Some Applications to the Rheology of Viscoelastic Fluids," Rheologica Acta, 5, 177 (1966). 21. Graebel, W. P., "Stability of a Stokesian Fluid in Couette Flow," Physics of Fluids, 4, 362 (1961). 22. Graebel, W. P., "The Hydrodynamic Stability of a Bingham Fluid in Couette Flow," in Second-Order Effects in Elasticity, Plasticity, and Fluid Dynamics, Reiner, M. and Abir, D., Eds., Pergamon Press, Inc., New York, 636 (1964). 23. Harris, D. L. and Reid, W. H., "On the Stability of Viscous Flow Between Rotating Cylinders, Part 2, Numerical Analysis," J. Fluid Mech., 20, 95 (1964). 24. Huppler, J. D., Experimental Determination of the Secondary Normal Stress Difference for Aqueous Polymer Solutions, Ph.D. Thesis, University of Wisconsin, Madison, Wisconsin (1965). 25. Leslie, F. M. "The Stability in Couette Flow of Certain Anisotropic Fluids," Proc. Camb. Phil. Soc., 60, 949 (1964). 26. Lewis, J. W., "An Experimental Study on the Motion of a Viscous Liquid Contained Between Two Coaxial Rotating Cylinders," Proc. Roy. Soc., A117, 388 (1928).

171 27. Marsh, B. D., Private communication. 28. Merrill, E. W., Mickley, H. S., and Ram, A., "Instability in Couette Flow of Solutions of Macromolecules," J. Fluid Mech., 13, 86 (1962). 29. Metzner, A. B., White, J. L., and Denn, M. M., "Behavior of Viscoelastic Materials in Short-Time Processes," CEP, 62, 12, 81 (1966). 30. Noll, W., "A Mathematical Theory of the Mechanical Behavior of Continuous Media," Arch. Rat'l. Mech. Anal., 2, 197 (1958). 31. Oldroyd, J. G., "Non-Newtonian Effects in Steady Motion of Some Idealized Elastico-Viscous Liquids," Proc. Roy. Soc., A295, 278 (1958). 32. Perry, J. H., Chemical Engineers' Handbook, 3rd ed., McGraw-Hill Book Co., Inc., New York (1950). 33. Pipkin, A. C. and Owen, D. R., "Nearly Viscometric Flows," Physics of Fluids, 10, 836 (1967). 34. Rayleigh, Lord, "On the Dynamics of Revolving Fluids," Scientific Papers, 6, 447 (1916). 35. Reiner, M., "A Mathematical Theory of Dilatancy," Am. J. Math., 67, 350 (1945). 36. Rivlin, R. S., "Hydrodynamics of Non-Newtonian Fluids," Nature, 160, 611 (1947). 37. Rivlin, R. S. and Ericksen, J. L., "Stress-Deformation Relations for Isotropic Materials," J. Rat'l. Mech. Anal., 4, 323 (1955). 38. Rubin, H. and Elata, C., "Stability of Couette Flow of Dilute Polymer Solutions," Physics of Fluids, 9, 1929 (1966). 39. Taylor, G. I., "Stability of a Viscous Liquid Contained Between Two Rotating Cylinders," Phil. Trans., A223, 289 (1923). 40. Thomas, R. H. and Walters, K., "The Stability of Elastico-Viscous Flow Between Rotating Cylinders, Part 1," J. Fluid Mech., 18, 33 (1964). 41. Thomas, R. H. and Walters, K., "The Stability of Elastico-Viscous Flow Between Rotating Cylinders, Part 2," J. Fluid Mech., 9, 557 (1964).

172 42. Truesdell, C., The Elements of Continuum Mechanics, Springer-Verlag, Inc., New York (1966). 43. Walowit, J. A., "The Stability of Couette Flow Between Rotating Cylinders in the Presence of a Radial Temperature Gradient," AIChE J., 12, 104 (1966). 44. Walters, K., "Non-Newtonian Effects in Some General Elastico-Viscous Liquids," in Second-Order Effects in Elasticity, Plasticity, and Fluid Dynamics, Reiner, M. and Abir, D., Eds., Pergammon Press, New York (1964). 45. Weissenberg, K., "Abnormal Substances and Abnormal Phenomena of Flow," Proc. of First Intl. Cong. Rheol., North Holland, Amsterdam, I29 (1948). 46. Wilde, D. J., Optimum Seeking Methods, Prentice-Hall, Inc., Englewood Cliffs, N. J. (1964). 47. Wills, A. P., Vector Analysis with an Introduction to Tensor Analysis, Dover Publications, Inc., New York (1958).