THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING ON THE STABILITY OF PARALLEL FLOWS OF ViSCOELASTIC LIQUIDS Dean T. Mook A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Engineering Mechanics 1966 October, 1966 IP-751

Doctoral Committee: Associate Professor William P, Graebel, Chairman Professor Jack L. Goldberg Professor James D. Murray Professor Chia-Shun Yih

AC KNOW LEDGMENTS It is a pleasure to express my appreciation and gratitude to my Chairman, Professor William P. Graebel, for his excellent guidance, timely suggestions and encouragement throughout the course of this work. I also wish to thank the other members of my committee for their time and interest. I also wish to express my appreciation to the Computing Center of The University of Michigan for the use of their facilities, the Industry Program of the College of Engineering for the reproduction of the final manuscript, the COMCOMP project sponsored by the Advanced Research Project Agency of the Department of Defense for the use of the on line interactive computing system, and the Fluid Dynamics Branch of the Office of Naval Research which partially supported this work under Contract Nonr-1224(49) with The University of Michigan. Last but not least, I wish to thank the members of my family whose sacrifices made it possible for me to complete my graduate studies.

TABLE OF CONTENTS Page AC KNOW LEDGMENTS........... ii LIST OF FIGURES. o..... i CHAPT ER I. INTRODUC TION........... 1 CHAPTER II. IMPLICATIONS OF THE CONSTITUTIVE EQUATIONS.......... 4 CHAPTER III. THE STABILITY PROBLEM... 15 CHAPTER IV. AN ASYMPTOTIC SOLUTION TO THE STABILITY EQUATION VALID FOR SLIGHTLY VISC OELASTIC LIQUIDS. e o... 31 A. The Determination of the Characteristic Equation for a Navier-Stokes Liquid.... 31 B. The Determination of the Characteristic Equation for a Viscoelastic Liquid... 39 C. Solution of the Characteristic Equation.. 43 CHAPTER V. CONCLUDING REMARKS AND SUMMARY.. 52 BIB LIOGRAPHY.............. 54 APPENDIX I...... 56 APPENDIX II. O. O......... 60 111

LIST OF FIGURES Figur e 1. Neutral Stability Curves for Various Values of X, Comparing with Results of Others.. 49 2. Critical Reynolds Number versus the Parameter R. 50 3. Neutral Stability Curves for Various Values of X and Curves of Constant c.... 51 iv

ABSTRACT SA. general linear hereditary stress-strain law is used in studying the stability of plane parallel flows of viscoelastic liquids. In I:eepiny/ with the usual procedure, a stability equation valid for liquids of arbitrary memory is derived by superposing an infinitesimal, two-dimensional, periodic disturbance on the primary flow. It is shown that in the case of liquids with short memory this stability equation reduces to the one given by Walters in 1962. An asymptotic solution to the stability equation is obtained using the method of inner and outer expansions, as presented by Graebel in 1966. Corrections to his expansions are noted, and the results are compared to those obtained by Chan Man Fong and Walters in 1965. The results indicate that the elastic property of the liquid has a destabilizing effect on the flow. In addition the sum of tensor quantities associated with a given particle at various times is discussed. A comparison of the sum obtained by adding corresponding components when they are referred to the same set of material axes and the sum obtained by adding corresponding components when they are referred to the same set of spatial axes is given.

CHAPTER I INTRODUCTION The interest in developing a theory for viscoelastic materials has grown with the widespread use of viscoelastic solids, such as vinyls, and viscoelastic liquids, such as polymer solutions. The latter has been considered for reducing ship drag by forming a layer of viscoelastic liquid along the hull. Granting this layer can be formed, one can then focus attention on the study of parallel or nearly parallel flows of such fluids, and a problem of basic interest would be that of the stability of these parallel flows. Any such study presumes the knowledge of a constitutive equation which adequately describes the material properties. From the many proposed equations appearing in the literature the one selected for this study is,DCg.,, (I. l) and -, "t' 2.C q 4 Z X tV,, x to X C" t. — ) g (I. 2) where the Pik are the components of the stress, p is a scalar pressure, gik the components of the metric tensor, Pik the components of the partial stress. The e(1) are the components of the strain rate. A -1

-2particle of the fluid occupies the position mrn at time t, and the position m x at time t. The relaxation function 4, (t - t) is given by ( _ -t) ('4t -7 T (I. 3) where N(T) is the distribution function of relaxation times. An equation of the type (I. 2) is a form of the Boltzmann linear superposition principle. Forms of (I. 2) have been discussed by Oldroyd [1], Walters [2], and Gurtin and Sternberg [3], to name a few. The integral form has the appeal of being general, but it also has the disadvantage of requiring a knowledge of the displacement field. (Indeed, in the only earlier studies of stability problems using this constitutive equation known to the author, Walters [2] as well as Thomas and Walters [ 13], working with (I. 2), were only able to determine the displacements for a short time interval and, consequently, obtained stability equations which only apply to liquids with short memory.) For this study (I. 2) was considered general enough to yield some meaningful results and yet simple enough to keep the problem tractable. In this thesis a stability equation which is valid for liquids with arbitrary memory is derived, and it is shown that for materials with short memory this equation reduces to one obtained by Walters earlier. An approximate solution to the stability equation is then obtained using essentially the method of inner and outer expansions.

-3The use of this method to solve the Orr-Sommerfeld equation has already been discussed by Graebel (1966). Chan Man Fong and Walters [4] published two solutions to the stability equation given in Walters' first paper using both a perturbation technique and an extension of the asymptotic method Lin used to solve the Orr-Sommerfeld equation. They found it possible to make only a qualitative comparison of the two results. In the present work, comparison in made with one set of their results, and reasons are given for believing the present results to be more accurate. In the paper byChan Man Fong and Walters as well as in this thesis Poiseuille flow is the only primary flow for which the stability equation is solved, although the derived stability equation is valid for all parallel primary flows.

CHAPTER II IMPLICATIONS OF THE CONSTITUTIVE EQUATIONS In formulating constitutive equations the following assumptions are generally made: First, the dynamic and kinematic variables involved in the constitutive equations must be independent of the rigid body motion of the material. This assumption is sometimes referred to as the principle of material objectivity or the principle of material frameindifference. Secondly, it is usually assumed that constitutive equations do not relate the dynamic and kinematic quantities associated with one particle to those quantities associated with a different particle, or, equivalently, the behavior of a particle is completely determined by the motion of a small neighborhood of that particle. Thirdly, the constitutive equations must be invariant under coordinate transformations. This means that the quantities involved must be absolute tensors. Moreover, whenever the histories of the quantities are involved, their components must always be referred to the same material axes (a material axis is one made up of the same particles of the material for all time) and the present properties of the material cannot be influenced by happenings that will occur in the future. -4

-5To aid in the formulation of constitutive equations consistent with the assumptions above, Oldroyd [1] proposed using a convected coordinate system which moves through space and deforms in such a way that the coordinate curves are also material axes. At any instant a one-to-one relationship is assumed to exist between the convected coordinates, i, and a system of spatial coordinates, x or x; namely, (II.la) or XL ~~~~~- s~ (~iSt ) ~(II. Ilb) where t > t. The spatial coordinate system used at time t may differ from the one used at time t, but the same arbitrary convected coordinate system is used at all times. At each instant (II. 1) may be regarded as a coordinate transformation. The ~ may be eliminated from (II. 1) to give (II. 2a) and -' -'(II. 2b),i Thus, the x represent the position at time t of the particle which is in the position represented by xi at time t. Whenever the velocity field has been prescribed, the relationships indicated in (II. 2), collectively referred to as the motion, are obtained by integrating

-6a-~'(II. 3a) and A=t tra (II. 3b) where D/Dt denotes differentiation with respect to time holding the xj (or equivalently, the i,) constant and 8/at denotes differentiation holding the x constant. One approach which can be used to develop constitutive equations is to represent the properties of the material symbolically through the use of one-dimensional rheological models and then to deduce, in a manner consistent with the assumptions stated above, the constitutive equations from the equations describing the model behavior. (The latter deduction is of course not necessarily unique. ) For example, there are two classical models associated with viscouselastic materials. One, the Voight model, has an elastic element in parallel with a viscous element. The other, the Maxwell model, has a viscous element in series with an elastic element. The former characterizes what has been called a viscoelastic solid and the latter characterizes what has been called an elastico-viscous liquid. For the latter one would have in mind a material that flows under any stress, that shows some strain recovery upon rapidly removing the stress, and that exhibits stress relaxation. (Stress relaxation occurs

-7when the stress decreases while the strain is held constant. See, for example, the photographs on page 129 in [5].) For the Maxwell model, the forces in each element are equal and the total elongation is the sum of the elongations of each element. Hence the differential equation describing the model is,F F C where F denotes force, E elongation, k the spring constant, c the damping constant, and a dot differentiation with respect to time. F may be expressed as a function of E as follows: ____ t t-t T(= tTe- e T ~ T""d where T = k/c. By associating stress with force and strain with elongation, a possible constitutive equation generalizing this is tot t t- t I T' (,TI - eF i~tl Kp i~tXe e 7d (Ct)dt (0 4) where T (called the relaxation time) = rj/L, with rn denoting the elastic (1) modulus, and pL the viscosity, rr' and E the components of the partial stress and strain rate tensors, respectively, referred to a convected coordinate system. (Henceforth, the components of any tensor will be denoted by Greek letters with Greek indices when

8referred to a convected coordinate system and by Latin letters with Latin indices when referred to a spatial system. ) It follows from (II. 4) that the stress relaxes exponentially with time. By performing the integration indicated in (II 4) with the t constant and with E(1) always referred to the same material axes, the assumptions of the initial paragraph are satisfied. For M Maxwell models connected in parallel, the stress corresponding to the composite model is the sum of the stresses corresponding to each of the Maxwell components, namely, 2 \ ~ t t:-i FTp be 2X Tn ) n -\ - oO For a continuous distribution of Maxwell components connected in parallel with a spectrum of relaxation times, the limit of this constitutive equation is t FFO<Fb Dirt)=d 2 )t-)CAi(,ta (11.5) where c0 t ee itr, (Il. 6) has been introduced, N(r) and 4, (t - t) being as defined previously.

9The Maxwell elements corresponding to the limit values of the spectrum have degenerate forms. At T = 0 there corresponds a purely viscous element and as T approaches oo a purely elastic element. Letting IL or nT approach infinity is equivalent to replacing the damper or spring by a rigid member, and setting L or T equal to zero is equivalent to removing the element. Hence (II. 5) may be associated with a combination of (either'or both) Maxwell and Voight models connected in parallel; that is, (II. 5) is a general form capable of describing the properties of viscoelastic solids as well as liquids. Gurtin and Sternberg [3] have shown that (II. 5) is the most general linear, hereditary stressstrain law for small strains. Oldroyd [ 1] has defined a convected differentiation, here denoted by 6/6t, as the inverse of the integration indicated in (II. 5). That is, the convected derivative is the tensor whose components are the derivatives of the convected components with respect to time as measured with respect to the convected axes. He has also shown th how to obtain the components of the convected derivative of an nthorder tensor in a spatial coordinate system. The results for a second order tensor are Ti rt n (tr. 7) This result will be used in Chapter IV. (Oldroyd has published a comprehensive discussion of the various time derivatives which are

-10used in the theory of continuous media in the appendix to [6]. There are other convected derivatives which could be used in constitutive equations (e. g., Jaumann's definition) which would result in quasilinear equations similar to (I. 2). ) (II. 5) implies a definition for the sum of tensor quantities associated with different points in space, but with the same particle. To illustrate this, the following example is considered: suppose the strain increases from zero in (say) two finite steps occurring at times tl and t2. These strain increments are labelled AE c (tl) and / E (tz). According to (II. 5), the resulting partial stress is /\\, I- L%(' t~)LCcpl\ r /i(+t)Z\C.1r W) ) The term inside the brackets is the sum (apart from the scalar factors 4 (t - ti)) of the corresponding components of the strain rate at the two times tl and t2. Since a convected coordinate system is being used, the strain rate is associated with the same particle and the components are referred to the same set of material axes. In general, however, the particle occupies different positions in space and the material axes have different spatial orientations at the two times. This suggests that the sum of tensor quantities associated with different points' in space, but with the same particle, is obtained by first referring all the components to the same conrvected coordinate system, where the tensors to be summed are associated with the same point, and then

- 11- 1 adding corresponding components. These sums are then the components of the sum of the tensor quantities. This procedure was apparently first proposed by Oldroyd [ 1] and is not a generalization of the usual procedure for summing vector quantities. The usual procedure for adding vectors is to "shift, " or parallel transport, the vectors to a common point and then to add corresponding components: these sums are then the components of the sum. More generally, then, the addition of tensor quantities requires shifting the tensors to a common point in space and then adding the corresponding components. The following example is considered in order to illustrate the difference in the two procedures. A second order tensor is used, but the extension to a tensor of any order is immediate. Suppose B represents a second order tensor associated with a given particle. At time t the position of this particle is denoted by xi and at time t by x1. The components of the tensor are denoted by b ij(x, t) and bij (x, t) at times t and t, respectively, when they are referred to a spatial reference frame and by C (g, t) and PIp( t) when they are referred to a convected reference frame. Accordingly, the first procedure yields for the components of the sum and the second yields

-12Sj = a,t 2 i j b(k x,') where k the G are the contravariant base vectors at x andg are the covariant base vectors at x. The g.k are often referred to as "shifters" and were apparently first discussed by Doyle and Ericksen [7]. They are related to the displacements through the following equation where the U are the components of the displacement associated with the G. and I represents covariant differentiation. In order to show that the 6.. and the CP are components of different tensors, it is convenient to transfer both sets of components to the same spatial reference frame; they then become Sji i x,t. ) bqn lt, t)u/ bne:._ X, -. tt N a+ 4n ( -u\ dv(Rl IJ It - and A j, =4 JXX

-13The two expressions are quite different, the difference being that the sum as proposed by Oldroyd is formed by adding corresponding components when they are measured with respect to the same set of (deforming) material axes and base vectors, and the other sum is formed by adding components when they are referred to the same set of (undeforming) spatial axes and base vectors. Oldroyd's definition is consistent with the idea of material objectivity, and the other is consistent with the familiar procedure for the addition of vectors. The author has found no other discussion comparing these two procedures for forming a sum in the literature. Since in fluid mechanics the boundary conditions and equations of motion are generally referred to a spatial reference frame (Lodge [8] has worked out some examples in which the boundary conditions, etc. are referred to a convected coordinate system), it is desirable to refer the constitutive equations to a spatial coordinate system also. The result is t -- t61ij - 2 i 9 ( t -t' ) J e x, c))t (II. 7) or, if the contravariant form of (II. 5) is used, t,' X = z f 1 ( t - t ) -h _ * - ) (II. 8)

-14(II. 7) and (II. 8) are not equivalent. Even in a cartesian coordinate system the two stress components Pij and pij are generally not equal, as will be seen in the following chapter where a steady flow example is considered. (II. 7) and (II. 8) must be regarded as two different constitutive equations and, as Oldroyd has stated, an experiment is necessary to determine which one should be used for a given liquid.

CHAPTER III THE STABILITY PROBLEM The solution for steady flow between parallel plates is next worked out using both (II. 7) and (II. 8). The stability problem for parallel flows is then formulated for each. The stability of a steady flow is studied by superposing a wavy infinitesimal disturbance on it and then determining under what conditions this disturbance will grow. Cartesian coordinates are used, the z-axis and the y-axis being chosen parallel and perpendicular to the plates, respectively. For the steady flow the velocity components are assumed in the form L= V=o, w=W(y (11,. 1) By inspection the motion is i<= X, (III. 2a) -- = Ad J(III. 2b) and (HI. zc) i - w(We~- t ) (III.Zc) The non-zero component of the strain rate is e 2Dw (III. 3) -15

-16where D represents differentiation with respect to y. The non-zero displacement gradients are (III. 4a) ~-~ _ I (III. 4b) 3D - and 45 \ E DWvB -o (III. 5b) where B \T h_(T) dT (III. 5c) Substitution of (III. 5) into the equations of motion and assuming W (4h) = O result in W= (\- (I[1.6) The only non- zero total stress components are'ZO - \^ 5(III.7a) _w IT- (I.5a

-17(III. 7b) and - (III. 7c) where 2p WCVW i + + crV\s CA (III. 7d) corresponds to the stress found in a Navier-Stokes liquid. Using (II. 8), the only non- zero components of the total stress are 6X VI I' (III. 8a);~.'2.u, y (nI. 8b)...- (III. 8c) Hence, (II. 7) and (II. 8) predict different normal stresses, the forms given in (III. 7) and in (III. 8) differing from that corresponding to a Navier-Stokes liquid in the terms containing B1. The velocity fields, however, are the same for all three.. These normal stress differences are consistent with the sudden expansion or contraction of the stream when some non-Navier-Stokes liquids suddenly emerge from a pipe into the atmosphere (the Merrington effect), as well as with the differences in the shape of the free surface for such different liquids

-18undergoing Couette flow (the Weissenberg effect). In the development of the Orr-Sommerfeld equation consideration is limited to a disturbance that corresponds to a velocity field which is both temporally and spatially periodic. Subject to this limitation, Squire [9] has shown that it is sufficient to consider a disturbance that corresponds to a two-dimensional velocity field. In the viscoelastic case, only such disturbances will be considered also, although no proof of the sufficiency of this exists. Accordingly, the disturbance velocity components are taken in the form a* = c, C * = steB, avd ue= -rlyB (In. 9) where E = exp i A (z - Ct); A is the (real) wave number, and C the (complex) wave speed. The components of the total velocity are tk = =u ot*, vd X = W + J (I1I. 10) The total displacement is assumed to be the sum of the primary flow displacement and the displacement resulting from the disturbance; namely, X = X, (III. 1la) t b + ~((~t ~it, ~ - ~ ~i)~, ) (HI. lb) and wh=ei-re (ta-n) + rersntepstioa time, ) (III.o lc) where y and z represent the position at time t of the particle in the

-19position represented by y and z at time t; subject to the conditions ojt = tb - O. (III.12a) and It (3 l t = )=o. (III. 12b) Substituting (III. 11) into (II. 2) leads to the following equations for a and p tr w-Ai) - 0i (III. 13) and ve-w-w~~ t(WtSE)~= 0. (11.14)..E - w- ~ - ). (I. 14) (Whenever the argument of a function is not stated explicitly, it is understood to have its value at time t. ) Assuming the various derivatives of a and p are of the same order as a and p and neglecting terms containing products of the small quantities, (III. 13) reduces to -vE.&~~ a3t +W~?~ b (III. 15) and (III. 14) reduces to -He W teA-W am ~W; (HII. 16) Assuming the form a = f (y) E + h (y, z,F) and considering (III. 12a) lead to O — ~E ~,- ex AtA-Ct)+ (III. 17)

20It follows from (III. 11) that __..' )__. (=,- )- *' iLO) - WV +,j *.t and C) ) )V t -.) t \ )I (III. 17) can be rewritten as -}\.I-. t_ IF.. (III. 18)'LA(w'- A where F = expiA(C - W) (t- t), and (III. 16) becomes -o -i-',:g - smL la r- r J (11. 19) Follwing anc miart -A apo n -c)o'~.r v,/.x,- -YV) + K; w 1_ -- K ) (11. 20) l inearization exact solution is obtained for a which justifies the linearization used in the present analysis. If consideration is restricted only to a short time interval, then (III. 18) and (III. 20) can be approximated by expanding F in a power series in (t - t) and retaining only the first order terms. Hence,

-21and - e- ( W+ E)(t-t). (1I.22) These are the expressions used by Walters [2]. The more general forms (III. 18) and (III. 20) will, however, be used here. Finally, (III. 18) and (III. 20) are well-behaved at the critical point where W and C have the same value. Applying the limiting process to (III. 18) and (III. 20) results in (III. 23) and 3'i (B - ->jt" L (~DVV) 1 r-'ai;) (III. 24) in agreement with forms (III. 21) and (III. 22), in the linear terms in t - t The non-zero strain rate components at time t are e t )L, (III. 25a) (III. 25c) I I. E 2- TDW (III 2-5b)

-22To the same approximation used for the displacements, the non-zero strain rate components at time t are,)JU= DT EJ F (III. 26a) et = 24AE +BD) EF - -FA w_ C) 2' EID (III. 26b) and e') LsEF(III. 26c) Making use of the continuity equation, the non-zero derivatives of the displacements can be written as ~ — = | - C E- ( \- ) (III. Z7 a). a~ w-c +Di \ s tS IDW E (t - t (f__w) EHI), (III. Z7c) and LD )W (t - t ) T zW E F _(III. Z7d)'A~ W W- A,()-C) _D an'W-:W-c)T W,'

-23Putting (III. 26) and (III. 27) into (I. 2) leads to t t )t 2 - W-.C This can be simplified since t t,c t -t I Tlsc, lrt J-~70 d i,'F -cc -o ~ N dT I +A F- C) T

- 24Introducing Hi (y) =(B - Ko)/iA (W - C) and P33 = Pzz E results in - - ZKoD - iADW, c(III.28a) In a similar manner the other non-zero partial stress components are -=< - D'T Y 2, D (III. 28b) and ZttDW iV (III. 28c) +2 DWDZ W 42 +H)- (ZDW))A(l34)/A - h WI\KV ~ii ~~ Ko -r LDW) iw -( )(H3 + 3 C H The functions K (y), H (y) and Gn(Y) are generalizations of Walters' constants and are defined by o0 ____ T_ __ __) JT (III. 29a) WX tE-= I -c~ir~ i o Ilr\s A( W-C) T (HI. 29b) j._ -

-25and r ~ N lT O[tA C) (III. 29c) From (I. 1), the non-zero total stress components are (III, 30a) U3Y _s- 73 _ 2,~)2x = $_?D)E (I". 30b) (IlI. 30b) y2)8 = wt + q E(III. 30c) and tD; z =: (:3;t - ~)E (III. 30d) where p is the pressure associated with the primary flow and p (y) E the pressure associated with the disturbance, It is possible to obtain a partial check on the stresses as given above by considering a specific N(T) for which the integral equation can be replaced by an equivalent differential equation. The determination of the stresses in such a case does not require knowledge of the displacement field. For example, suppose whee =he \s rpren he Da (II. b31i where the 6's represent the Dirac delta function. It can be readily

-26verified by substituting (III. 31) into (II. 5) that stress and strain rate are related by the following differential equation'";* \' +,C, g C %?'r+ - (III. 32) where 6/,t represents the convected derivative of Oldroyd discussed in Chapter II. This equation often appears in the literature. Splitting the stress into -(. -- Pij -1- ~Y (III. 33) where the P..ij result from the primary flow and the EPij from the disturbance and then substituting (III. 33) and (III. 10) into (III. 32) yields, after neglecting terms O (v2), P, -t- t7 tl t + -i Ps (t - C) ) \ 2 i -F )\ t > 1:, t-1 P (IIIo 34a) -lAW),. I k-l ji) X (l/? ) Y\N -?'+ 3 3319 1-1 Mv-c) xi - E D >JII.DP (III. 34b) i 2 i+ /\ X. { l + —- C..) i, k, -A, k- c X, (,;,-:~ + X.q zL. /\-,'~ ~..O,

-27and'Al3 C,,63[\4-iT' ) 3 >\I)P (III. 34c) +' /\ +?V /., a + 2,,\z-] F- - ~,3 A-j.. +,DW -7__ + A mgT_ +w. + IiA(\'- A + )E E From (IlIo 5c), (IIIo 30) and (III 31) it follows for such an N that (III. 35a) __ __ _ 2j >9 _ ) i(III. 35b) 8\c. I = A_2 \ t, L, (III. 35c) (III. 35d) and _____ 1)( ) (' cosd__g f X _' (II. 35d) Then considering terms of the same order and using (IlL. 35), the resulting P,. agree with (IIT. 5) and the p O with (ELI Z3E o)8).

-28Substituting (III. 30) into the equations of motion, one obtains, after considering the equations satisfied by the velocity components of the primary flow, i~ \1\- C)- v -" -3C U T' and - p tw- C _D ) Laz /A -...\..t.A.i Eliminating p by cross-differentiation leads to D3lt 0)(~~ bV4 \ )_\fr\N/j\5' = ~t +AJIC ) (III. 36) In order to write (III. 36) in terms of non-dimensional variables the following dimensionless functions of y are introduced (III. 37 a) -4, (III. 37b) \.\. (III. 37c) \=) _ gL<V~3 ~(III. 37d) where V is a characteristic velocity and L is a characteristic length. Putting (III. 28) and (III. 37) into (III. 36) results in

-29-'i(\J-c)(D'i)-t~i) 2 | (izv 0l)% VT (III 38) dimensional wave speed, C/V; U is now the non-dimensional primary flow velocity, W/V and D represents differentiation with respect to +' n + o{(~-c)T' i ~3y i~+)) \' - where a is the non-dimensional wave number, AL; c is the nondimensional wave speed, C/V; U is now the non-dimensional primary flow velocity, W/V; and D represents differentiation with respect to the non-dimensional y. For viscoelastic liquids whose material properties can be adequately described by a linear, hereditary stress-strain law, (III. 38) is the stability equation for parallel flows. (For a Navier-Stokes liquid all the Pn are zero except P and (III. 38) reduces to the Orr-Sommerfeld equation. ) Walters [2] has suggested that for some viscoelastic liquids, called "slightly viscoelastic, " it is reasonable to expect the N(T) to vanish rapidly as T increases. (This is the justification he gave for using (III. 21) and (III. 22).) In this case the upper limit in the expressions for the Pn can be replaced by a finite limit (say T) which may even be quite small. Assuming I iA(W - C) T << 1/T everywhere in the flow field and neglecting all Bn for n > 2, (III. 38) reduces to

- 30= il(U - C) ( D 2_ K z )- D a u 2 ff (III. 39) I- oDR,(u-c) _ --- L) + - Lz C where VL (III. 40a) and Rn - -', -/ -Y\ (III. 40b) (III. 39) is the equation given by Walters and was derived by using (III. 21) and (III. 22) in place of (III 19) and (III. 20). Two final remarks are appropriate. (III. 38) was derived without specifying the form of W (y) and hence is not restricted to the primary flow considered at the beginning of this chapter. Also if (II. 8) had been used in place of (II. 7), the resulting equation for v would have been exactly the same as (III. 38) - a rather surprising result, considering that primary flows involving the two constitutive equations give quite different normal stress results. An approximate solution to (III. 38) is given next.

CHAPTER IV AN ASYMPTOTIC SOLUTION TO THE STABILITY EQAUTION VALID FOR SLIGHTLY VISCOELASTIC LIQUIDS An approximate solution to the stability equation (III. 38) for the primary flow U = 1 - y2 is next obtained. The procedure used to determine the characteristic equation is due to Graebel [10], with only slight modifications. The determination of the characteristic equation for a NavierStokes liquid is discussed in Part A, the determination of the characteristic equation for a viscoelastic liquid in Part B, and in Part C the procedure used to actually solve the characteristic equation and determine the points on the neutral stability curve is presented. In the interest of brevity, some of the results for Part A are merely stated; in these instances the reader is referred to the paper by Graebel for the details. A. The Determination of the Characteristic Equation for a NavierStokes Liquid It is anticipated that both c and l/oR will be small for the case of interest. This suggests an asymptotic solution. Hence, after Graebel, the flow region is divided into an inner and an outer region. The inner region is a strip that includes the rigid bottom boundary at y = -1 and the "critical point" at y = Yc, where U(Yc) = c. The outer region extends from the edge of the inner region to the centerline between the plates at y = 0 which, as a result of the symmetry of the geometry and equations, -31

- 32serves as the other boundary. It is assumed that the two regions overlap. The procedure now is to obtain a solution valid in the inner region and a solution valid in the outer region and then to merge the two solutions. The inner solution is used in satisfying the boundary conditions at y = -1 and the outer solution is used at y = O. Satisfaction of all boundary conditions and the merging results in a characteristic equation which defines a as a function of R, the plot of a vs R. for c. = 0 being the neutral stability curve. In the inner region v is obtained by introducing the change in variable (coordinate stretching) ____~ - ~ = - (IV. 1) and by putting -(_z = =-l() (Iiv"?. 2) where [L is a function of aR, expected to be small, but unknown at this point. Substituting (IV. 1) and IV. 2) into (III. 38) with all of the Rn for n' 1 set equal to zero,the proper choice for ~i is = (ca R)~ (IV. 3) and (III. 38) becomes di X Xd (2&) d X (IV. 4) A ~ -+ a t- ~i q _d 1+'''= -0 d~l dr\c Y7

- 33X(0) and X(1) are the solutions of -- uc- d = ~ j (IV. 5) and dd - (L:- izc dX=t(Y dd X 2 > t~} ) q(IV. 6) Hence, ~X"=r \ = Y(IV. 7a) and ) I {d S t~ 00 jDVtc)i/3 z -\ (IV. 7d) where -RAll(DU)C) i- 3, (DD-ccin ) (1) (2 \3 )i'i ) 1- t ) H(I) and H(2) are Hankel functions of order one-third~ 3 3

- 34Since H() increases exponentially with large positive r at a much ~~3 ~~~~~(0) faster rate than does the outer solution, X4 cannot be merged with the outer solution and is therefore discarded. Using (IV. 7) in (IV. 6), it follows that _-X _z (IV. 8a)'D U' Graebel gave the exact solution for X1) x2 z (,.tX J x- (IV 8b) where the path of integration is chosen so that the integral does not grow exponentially with'o - oo For large rB, (tDc) - IA 4( \) * t (IV. 8c) where _ cI A.. In addition Graebel gave the recursion relation L -tn ~t \nt a')) -- D n ----- dl —Cn (IV. 8d) which leads directly to X =,, o 9 3 A ~) (IV. 8e) 3' Doc(O) v0 c'> ~ (-1,. OV l(-2 "............ __~ -k ___ -- IV e

- 35where _ 0.6430 i, 39 t 0 In the outer region v is obtained by introducing the expansion i' hA,y)- &. (XL(~U.(OI ~y) 8)+,(bl f'''li~' (IV, 9) Substituting (IV. 9) into (III. 38) gives i(UC)( D- ol)'C) LO) C (o) (IV. 10) and v(0) v(n) for all n such that E (p) > 0 (L3). The Tollmien solutions n [11] are lO-'~~ (1/- ~', ~1 — F A(~ I,"- (IV. lla) and 00oo o a i)= F; ()-,) t C \ (IV. llb) where Ao' I, A \ -- 0(, A - 4 _ and3I( 4) -A,)C and Ro —c,, o - -? _;**

- 36Bn )(v3) Btz- DUr (zn+5)8%z t (;%t3)An+ 1 + D c34tl - Q+;n for n > 0. In the inner region Graebel used C, YY\`~ CZ +~ C3~ 3 (IV. 12) where the ci are arbitrary constants. Then after considering (IV. 9) and (IV. 11), it appears that for the inner and outer solutions to merge for large positive r1 the proper choice is 1/[1 for:E,ln pi for El, and 1 for Ez. In the outer region V is given by - L) _ ~CC,()_, aC+ F-Z (IV. - 3) where A is a constant that cannot be determined until higher order (0) terms are considered. X3 decreases exponentially with large positive r1 and, hence, no terms are needed to merge with it. The boundary conditions are ]- _ = asdo -= O (I'V 1 4) and 3 — O ~ ~ ~ v o(IV.Ia) if V is an even function of y, or Wt~- =:'^- ~ it <O (IV. J 5b)

- 37~ if v is an odd function of y, at y = 0. Since V satisfied the OrrSommerfeld equation and U is an even function, V can be separated into independent even and odd parts; and since v in the outer region satisfies (IV. 10), only one condition of either (IV. 15a) or (IV. 15b) is sufficient. Thenbyusing (IV. 14) and (IV. 15a) since it has been shown that symmetrical disturbances are less stable than asymmetrical disturbances, Ct DF, )Xta c i t ) (IV. 16a) C, _______ _l -) C 28 A ft (IrV. 16b) and C + C 3 a - ~(IV. 16c) to the lowest non-trivial order in Ai, where ll = -(1 +yc)/ zl/ll and Zo - yc. For a non-trivial solution of (IV. 16) to exist, \i D XDFt 1 )AIV. 17) which is the characteristic equation obtained by Graebel. One finds, however, that the plot of a vs R for ci = 0 obtained from (IV. 17) does not contain the characteristic loop of a neutral stability curve. This shortcoming of (IV. 17) can be explained by considering the magnitudes (1) of the X terms in the expansion (IV. 2), which in turn requires some

- 38knowledge of rl. With Lin's results [12] (rl 2. 5, l 1/20) as a guide it appears that X(0) may not adequately describe the solution near (1) (0) y = -1. From (IV. 8), Im(BLX( )). 3 as compared with the ImX = 0; (1) (1) (0) ) the real parts of X X1) and ILX ) are small compared to X1) and X(0) respectively; and a check shows that LX(1) is small compared to X (0) Thus by including the X() terms only the imaginary part of X2 is altered significantly. This additional term apparently is responsible for generating the loop in the stability curve. It seems now that a better choice for the solution in the inner region is c > = C - (IV. 18) + Q (,A k c,v) 1 C3X Although the X1() term does not significantly alter X1 as discussed above, including this term aids in the merger with the outer solution. For merging (IV. 18) with the outer solution, it appears that the proper choice for Eo is 1/i, for el In p, and for E 2 1, as before. Then V(y) in the outer region is given by Subs =tC, iin F (IV. 19) Substituting (IV. 18) in (IV. 14) and (IV. 19) into (IV. 15a) leads to, for the lowest non-trivial order,

-39C,)K () +i CZX2 (i>) ~+ Cz X3 (U,> O0 (IV. 20a) c~/,(v/.,) C1(7 (~,] ( O, (IV. 20b) CI 2~ C, CF. Vito) - aC DF,,C; (IV. 20c) Ik C~~D Uc. and to the characteristic equation: (0) X3 (,) _L (IV. 21) since now F K -:),. %, (j and fz (). ( ) are reasonable approximations. The plot of a vs R resulting from (IV. 21) does have a loop. (IV. 21) is essentially the equation used by Lin, although he elected to express the outer solution as an expansion in powers of a. With modern day computers the expansion in terms of the coordinate rather than az seems to be much simpler and more accurate. B. The Determination of the Characteristic Equation for a Slightly Viscoelastic Liquid For a slightly viscoelastic liquid the solution of (III 38) can be carried out in a manner analogous to the solution for a Navier-Stokes liquid. Specifically, it is anticipated that both c and IPn/al will be small. Hence, the flow region can be divided into an inner and an

- 40outer region as in Part A. In the inner region the Pn can be expressed as series in i by expanding the denominators and then integrating term by term. The results are P. = F {\ - )iDRft'c - t _ ( ARo 2 c -(IV. 22a) PR R _R_ _i R>U0c t - D ucLK (IV. 22b) r, (h- iER3ju-iR D)c. (IV.22c) and I3 = R {RW - 3ibR5DsCt 3 D 5 32c (IV. 22d) 3o-D ) where DnU = DnU(y =c) (IV. 22) are valid in the neighborhood of the critical point for any distribution function that vanishes as T becomes

- 41large, specifically for N(T-) negligible when T > T where jiA (W-C)T <1. As in Part A, the change in variable (IV. 1) and the expansion (IV. 2) are introduced. When (IV. 22), (IV 1) and (IV. 2) are substituted into (III. 40), the proper choice for pi is, as in the case of a NavierStoke s liquid, A -/L (o~R)-14 ~ (IV. 23) and (IV. 4) and (IV. 5) are modified to read +,"'t 0 and d ~ X o(IV. 25) _ A.__C.z. M)- -C- O M; v where E = LaR1. (IV. 6) is unchanged. Since the term containing R1 multiplies the fourth derivative, X1 and Xz remain as in Part A. Taking E to be small, X3 is approximated by oX3(9) c/i7, A- ~) (IV. 26)

-42where d(9O i Jm Q{C8 O(IV. 27) and d 41~~~~~4 _ ~ ~ ~ ~ __ Hid(IV. 28) Then co is the same as X3 in Part A, (IV. 7c). After using (IV. 8d) ag ain, -, = 5 tn d"~bo s j ) ~(IV. 29) Near y = -1, Ec and X1)X3 are of the same size for the values of E used in this study; hence (LX3 could produce changes in the neutral stability curve equal to those produced by E+. However, for this study it was considered sufficient merely to determine whether or not the elastic property of the material (manifested in the term containing RI) has a stabilizing effect on the flow. Thus X1 and X2 as given in Part A and X3 as given in (IV. 26) are used for the solution in the inner region, and, consequently, the present results serve only to indicate a trend. For larger values of E a perturbation scheme such as (IV. 26) may not be adequate. In this case one would resort to an exact solution of (IV. 25). See Appendix II. Since the outer solution remains unchanged and X3 is not directly involved in the merger, the characteristic equation retains the same

- 43form as (IV. 21), the only difference being X3 is now given by (IV. 26) instead of by (IV. 7c). Further, (IV. 21) was obtained using (IV. 15). Since the approximation in the viscoelastic case is based on the equation (i U-c)L (a) z A W =,0 8 -R~iL jOt RU-c) v can still be separated into independent even and odd parts. It is assumed here again that symmetrical disturbances are still less stable than asymmetrical disturbances. C. Solution of the Characteristic Equation For calculation purposes, it is convenient to introduce the change in variable (IV. 30) where:= (D 5)c) Then (IV. 21) becomes |8X~d- i,(i )+\l4, dr g0 00 (IV. 31) _ _ _ t,,' ( (DF o(- F- F-. WD,' o 7 = - - I-

- 44where cO c DwMw/3 IS-r INcT)d iC sph: N~o t/ )) The functions of hi and h2 are discussed to some extent, and tables of hi and hz and of their derivatives are given in [14] In order to make use of these tables, it is convenient to express hi as 4, t) =<7 t+ \ ( Ca-~5; ) (IV. 32) where + = Oo - aX3 +a- xO + OX3v " and and where Y3 _= -a~'-,agd b =Z v4/-) (oM-\).W) 3yv\ ~,,'+1)

-45Then for x - i;, it follows that *e ( * ) = (z oo)*,* i+ Xm l 2_w2,: ( l5\ A 2v i3 v I-o,, + I ioo o)' and -;-o<4 @ v t l,) 2 t ta 00) (io 34-'',,3 L, L oZ,) (zoo-,. where Am = am (-200)m and Bm= bm (-200)m. The Am and Bm are given in Table I in [14]. The derivatives of h1 are determined from g(IV. 32) by term by term differentiation. For the integration indicated on the left hand side of (IV. 31) the series (IV. 33) and (IV. 34) were integrated term by term for <' 5, For t. > 5, the integration was continued numerically using Simpson's rule and the asymptotic expansion h1, until there were no noticeable changes in the value of the integral in the seventh significant figure. From [14] the asymptotic expansions are

-462 )?e) TY e Cos -\I 3(IV. Lt (:v.35) and _A +i x) -_. C e; tt (IV. 36) where O (zY)( 3 /') LTrA) These are valid for The results of the integration are 3 dPRe`Riti4)= 0.504361 (IV. 37) 00 W~ \[ ( 1)Zmv I BiiT ( o 41 *I z (-,) RZO,,

* |' V ~,,, ) (r, -I + ~' X -0 (, g +9 ) 0=D LAr Zo A'Z \d \$~~\pue SIZI44Lz. " *)'V tU )' z (6~ AI) ((e) >), 5 g( o-) 2 c8~ Ai) Zbd oG6 O 6i - = ( A 81'\- al 5 o 5.

- 48 For given values of X and c the solution to (IV. 31) was obtained by plotting the real part against the imaginary part for each side. The intersections of these two families of curves gave a and (>l. R was computed from _ SUc i-' c ), (IV. 41) The results are shown in Figure 1. The graph for X -0, which corresponds to a Navier-Stokes liquid, is seen to be in close agreement with Lin's results, and the preceeding statements regarding the anticipated size of the various parameters are seemingly consistent with the final results. The results are qualitatively in agreement with those of Chan Man Fong and Walters, their X being defined as five times the value of the present one. A comparison of the two sets of results is shown in Figure 1, the X being that of the present paper rather than that of Chan Man Fong and Walters. The disagreement of the two results is not understood. However it is noted that Chan Man Fong and Walters results for x = 0 depart considerably from the accepted results of Lin. On this basis, it is believed that the present results are the more accurate.

-492.2 1 ~ I 1' 2.2 I I I ~~~~~ ~ ~~~~~~I I! "'I' I " 2.1 2.0 / - / " / 1.9 / 1.8 / 1.7 1.6 1.5 0.02 1.4 ~ =~~~~0.025 a 1.2~ ~ ~ 00 0 0 0 00 0 0 00 0=.0 1.2 1.1~~~~~~ 0.005 XZo 1.0 \ 0.9 004 NEUTRAL STABILITY CURVES 0.8 - CHAN MAN FONG WALTERS ooooo LIN 0.7 - PRESENT RESULTS 0.6, 0.5 ~,~ ~~~ ~~~~~~~~~~~~~~~~~~~ z 0.012 0.4 0.3 - 12 13 14 IS 16 17 is 19 20 21 22 23 24 25 26 27 28 29 30 31 32 R 1/3 Figure 1. Neutral Stability Curves for Various Values of X, Comparing with Results of Others

-505500 I I 5000 4500 Rcr 4000 3500 3000 2500.. 0 i 2 3 4 R,xe 2. Critica l Reynolds Number versus the Parameter R Figure 2. Critical Reynolds Number versus the Parameter R

1.8 I I! I I _i. 7 - ~~C:0.30 C:0.29 C =0.31 C=0.29 C=0.32. C=0.26 1.6 - C=0.25 1.5X-0.02 1.4 1.3 1.2- X-o.oo 1.0 0.9 0.8 0.7 0.6 0.5 12 13 14 15 iS 17 18 19 20 21 22 23 24 25 26 27 R /3 Figure 3. Neutral Stability Curves for Various Values of X and Curves of Constant c

CHAPTER V CONCLUDING REMARKS AND SUMMARY Using a general linear hereditary stress-strain law to describe the material properties of a viscoelastic material, the stability of parallel flows was studied. Following the usual procedure, a stability equation valid for liquids with arbitrary memory was derived by superposing an infinitesimal, two-dimensional, periodic disturbance on the primary flow. In the case of liquids with very short memory this equation reduces to the one given by Walters. The equations are also similar in the region of the critical point, where the disturbance speed is the same as the flow speed. An asymptotic solution to the stability equation was obtained using the method of inner and outer expansions and the results compared with those of Chan Man Fong and Walters in Figure 1. Since their paper includes only a brief description of their calculations, it is difficult to explain the rather considerable discrepancy in the two results, although the present results are felt to be the more accurate since they agree with accepted results for the Navier-Stokes liquid. Keeping only the first two terms, Lin as well as Chan Man Fong and Walters used a series expansion in ao2 for the solution in the outer region, or the outer (inviscid) solution. For a> 1 two terms may not approximate the series well, but the third term involves considerable - 52

- 53computation and was neglected by these authors. In this study a series expansion in the coordinate was used for the outer solution. This is more accurate since an arbitrary number of terms can be conveniently retained. The differences in the expressions used may explain the slight discrepancy for a > 1 between Lin's results (also shown in Figure 1) and the results obtained in this study for X = 0, but not the difference between Lin's results and those of Chan Man Fong and Walters for X = 0. All results for these viscoelastic liquids indicate that the elastic property of the material has a destabilizing effect on the flow.

BIBLIOGRAPHY [ 1] Oldroyd, J. G., "On the Formulation of Rheological Equations of State, " Proceedings of the Royal Society, Volume 200, Series A, 1950. [2] Walters, K., "The Solution of Flow Problems in the Case of Materials with Memory Part I, " Journal de Mechanique, Volume 2, Number 4, December, 1962. [3] Gurtin, M. E. and Sternberg, E., "On the Linear Theory of Viscoelasticity," Archive for Rational Mechanics and Analy.sis, Volume 1, Number 4, 1962. [4] Chan Man Fong, C. F. and Walters, K., "The Solution of Flow Problems in the Case of Materials with Memory Part II, " Journal de Mechanique, Volume 4, Number 4, December, 1965. [5] Fredrickson, A.-C., Principles and Applications of Rheology. Englewood Cliffs, N. J.: Prentice Hall, 1964. [6] Oldroyd, J G., "Non-Newtonian Effects in Steady Motion of Some Idealized Elastio-Viscous Liquids," Proceedings of the Royal Society, Volume 245, Series A, 1958. [7] Doyle, T.C. and Ericksen, J.L., "Nonlinear Elasticity'' Advances in Applied Mechanics, Volume IV. New York: Academic Press, 1956. [8] Lodge, A. S., "On the Use of Convected Coordinate Systems in the Mechanics of Continuous Media, " Proceedings of the Cambridge Philosophical Society, Volume 47, Part 3, 1951. [9] Squire, H. B., "On the Stability for Three-Dimensional Disturbances of Viscous Fluid Flow Between Parallel Walls, " Proceedings of the Royal Soceity, Volume 142, Series A, 1933. [10] Graebel, W. P. "On the Determination of the Characteristic Equations for the Stability of Parallel Flow, " Journal of Fluid Mechanics, Volume 24, Part 3, 1966, [11] Tollmien, W., "General Instability Criterion of Laminar Velocity Distributions," Technical Memorandums, NACA, Number 792, 1936. - 54

55[12] Lin, C. C., "On the Stability of Two-Dimensional Parallel Flows, " Quarterly Of Applied Mathematics, Volume 3, 1945. [13] Thomas, R. H. and Walters, K., "The Stability of Elastico-Viscous Flow Between Rotating Cylinders Part I, " Journal of Fluid Mechanics, Volume 18, 1964. [14] The Annals of the Computation Laboratory of Harvard University, Volume II. Cambridge, Mass.: Harvard University Press, 1945.

APPENDIC ES - 56

APPENDIX I AN EXACT SOLUTION FOR a For a linear primary velocity profile, an exact solution for a is developed, and it is shown that this solution agrees with (III. 18) and (III. 23). T ake =- r>(, t | iE It)~? (AI. 1) where 04 (t l t ~)=o -(AI. 2) defines y. Then t and After differentiating with respect to time, o( c. 8@(\1- C) (AI. 3) since d _ = \!/' T) and - ontinuity. v-v and __- Al from c5ontinuity. -57

- 58Equation (AI. 3) is valid for any velocity profile. For a linear primary velocity profile, \N = - =-~ + -\ - w (5) Q 0 (AI. 4) Putting (AI. 4) into (AI. 3) and integrating once gives ~' = i~l t Wt34 tA a + -L Ct | where 0o = O(t- t ). Then K o \ -S+ + \N%9)- C -_ i \ where S =\]w -w - ZOo Integrating again gives J=S ( ) ( - ) (Al. 5) where anLS(t-)t)-C and L=l~cs

- 59For W(y) - C i 0, S= W-c = e L= ( and Q = O. to the lowest order in a. (AI. 5) thus reduces to (III. 18) to this order. In a similar manner (AI. 5) reduces to (III. 23) for W (y) - C = 0. Other velocity profiles would give results much like (AI. 5) and the integration could be carried out in a similar fashion provided that it is practical to find the roots of W - C considered as a polynomial in a.r

APPENDIX II An exact solution of equation (IV. 25) is possible and has, in fact, been given by Chan Man Fong and Walters [ 4]. A modified version of their results is presented here to show the application of the present method. Equation (IV. 25) is of the form Writing and e the equation for f is the confluent hypergeometric form V I>) 0S a' J,.o (AII. I) and )= 2 i- e e ( -6 (AI0. 2) -60

-61= As z approaches ioo, U(z) approaches z and V(z) approaches z e. Thus it is suggested that U must be the solution corresponding to X3 and V to X4. (The difference in the asymptotic behavior of U and the Hankel funetions is somewhat indicative of the computational difficulties involved in computing U. ) To verify the first (and most important) of these, replace v in (All. 1) by -1/2 + i sT/2. Then in the limit taking argz = - Tr and P - 2ar/3 co e \rxt 5l- e- - Ut )2 ieex 3(A3) Expressing the right hand side of (AII. 3) in terms of Hankel functions gives the r equired e su3 (AII. 4) the required result.