‘ ‘4... I ‘H b‘ Id .I. .1 I \ U .l.‘ .v. . Ir 4 u L K II. I I‘I‘.‘ [4“ ‘ I O I o I. I _ ._ III to uni-“H1. I hurt; fiIBII‘INIaI... A . . - ‘ I. l .I .I. 1 I _ r I J numumflc . 0 I“... . l. 1 .. OI flo‘ho-llmllm‘ Q t .01.! O . Jiswjl‘ Il‘ll' ' . 'l . -. I. ‘4%.01Qllfi b OI. I I $0“)? ,3 . .l D. Q I I1 g . - . idfihuuu‘lvlldu. . II I. ' :II. .I -I.‘ In. I. II t I. chart? I _..II..| lIII . . .II .I. ‘I. .II |1|IQII‘#IIII\III‘III]II(I.III Ilrl‘lQIII. “¥!0 ‘0'] . .1 u " IJIh'II‘I'II‘IIL I‘ll-1.0..L1IIIII u Idnj‘iil "‘ll". ‘IOI O I I I 1"..-“ It“! ‘t‘ll Ei‘l It‘ll-ll. I‘ll-Ill" {Ni-I. III- P.“% Ilwll I'lr II“ I I3 I II- III- lll‘llll Ills} +Ill'llfi" I] I! I -I'It’lliulrl .I‘.‘ I III IIO-I.|.|.I‘ II-I Ulc‘fi‘S]... nuts... 7. I... I ll- iay‘ “E. ‘IIIIvE: tl‘fi" vl‘I-l II‘i‘I’I‘II‘IIIIOIII‘I‘t l.‘l10§|o..|.l‘ . I .300. -. 4'00“” I '15.. I. 1.11} _ IILI'III! I'll-‘1 -Il‘ul-|..‘|‘l'."lnl 1““4" ‘1‘ in! a - II . I ‘.II ‘III II. I- .‘l'... ‘9 l‘llilr" II wl‘ull‘IIII-Illl‘“ II II n alv... %1l .t.“o- 0“. III.-. .I. o ‘IQP‘... .Ill 1' ‘Vlloli I: I A I. I. II LII I ‘l‘llll‘l‘l 11.1.1113 I--. ll!IlliII.I‘III~I..|1!I1\IIH\III|I‘.I|Iai .LI.-\.II‘IIIQ. - Ill-Ill. Il'..l| it’l‘oxihhl‘lth'0IOI‘I-l‘ I 1-... {$5 '..!.|‘ .n‘lJ'i‘l.‘ ll.‘ 1 , I: II.- nl.“ OII toiOd- ‘0. .HI.QII-II|II‘O“\4|‘|.J l‘rllll“.>l‘l‘ll.hh“‘:l‘l'. ill-O I o ..l . I ‘0 .al -Ilol‘il Huu_:l_‘. -IIIII I:I“I‘I'm G o I In..- ‘ .b- all in. I‘ll-[I‘d THESIS This is to certify that the dissertation entitled TRANSLATION OF TWO CONTACTING SPHERES IN A VISCOELASTIC FLUID presented by Mohammad Amin Jefri has been accepted towards fulfillment of the requirements for Ph . D. degree in Chemical Engineering Date June 15, 1983 MSU is an Affirmative Action/Equal Opportunity Institution 0-12771 IVIESI_J RETURNING MATERIALS: Place in book drop to LJBRAfiJES remove this checkout from _—c—- your record. FINES will be charged if book is returned after the date stamped below. :don us: 0555 [3.5) F‘W‘f (5% M . «‘L "Fig. -. Giuh' Ed “ 0c: "‘\ :Q‘COQ I i..— / TRANSLATION OF TWO CONTACTING SPHERES IN A VISCOELASTIC FLUID By Mohammad Amin Jefri A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Chemical Engineering 1983 TRANSLATION OF TWO CONTACTING SPHERES IN A VISCOELASTIC FLUID By Mohammad Amin Jefri This work was undertaken to study the effect of a visco- elastic fluid on the translation of two contacting spherical parti- cles of unequal size in creeping flow. The elastic effect of the medium was investigated theoretically by employing a second order fluid model and developing a numerical scheme that evaluates the elastic effect on the drag force. The accuracy of the scheme was established by carrying out a detailed error analysis at each step of this scheme. It was revealed by analyzing the numerical results, that the region near the stagnation points at the surface of the larger sphere has the major contribution to the elastic effect on the drag, while the region at the contact point has no contribution. The results obtained for several size ratios of the large to small sphere where the numerical scheme is valid showed an appreciable elastic effect on the drag. This effect increased with increasing size ratio, up to a ratio of 3. Experiments carried out on the settling of these particles in a solution of 0.2 wt.% Separan in corn syrup at particle Reynolds 4 6 numbers in the range of 10' to 10' yielded the following results. First, it is seen that a stable orientation exists in the direction Mohammad Amin Jefri of the applied force (gravity) along the line of centers with the larger sphere underneath. Secondly, the deviation from Newtonian drag for equal spheres is zero. This agrees with previous theoretical results. In the case of unequal spheres, a 10 percent reduction in the drag coefficient below the Newtonian value is observed, at a Neissenberg number of 0.1. This reduction is seen to increase with increasing Neissenberg number. Good agreement was seen when compar- ing the experimental and theoretical results. To my parents--Zain Almasri and Amin Jefri To my wife, Hend To my son, Faris ii ACKNOWLEDGMENTS I would like to thank Dr. K. Jayaraman, my major advisor, for his sincere guidance, patience, and deep involvement in this dissertation. My thanks also go to Drs. E. Grulke, P. Wood, and L. Vallejo for their valuable comments. To my mother, father, wife, and son Faris, thank you for your support and encouragement over the years. I would also like to thank the Government of Saudi Arabia for financially supporting me during the period of my study. Finally, thanks to all my friends, especially Dr. Ekong Ekong, for his friendship. TABLE OF CONTENTS Page LIST OF TABLES. . . . . . . . . . . . . . . . vi LIST OF FIGURES . . . . . . . . . . . . . . . v11 NOMENCLATURE . . . . . . . . . . . . . . . . ix Chapter I. INTRODUCTION . . . . . . . . . . . . . . 1 1.1 Introduction 1 1.2 Particle Motion in Uniform Newtonian Flows 2 1. 3 Particle Motion in Viscoelastic Fluids 5 II. MOTION OF AGGREGATES . . . . . . . . . . . 13 2.1 Particle Interaction . . . . . . . 13 2.1.2 Effect of Particle Size . . . . . . 13 2.2 Hydrodynamic Effect . . . . . . . . . 15 2.3 Bulk Stress . . . . . . . 17 2.4 Objective of Present Research . . . . . . 18 2.5 Statement of the Problem . . . . . . . . 18 III. DRAG CALCULATION FOR TWO TOUCHING SPHERES . . . . 22 3.1 Equation of Motion . . . . 22 3.1.1 Reciprocal Theorem (Background) . . 24 3.1.2 Application of the Reciprocal Theorem . 29 3.2 Use of the Newtonian Solution . . . . 31 3.3 Results for Equal Spheres . . , 42 3.4 Numerical Procedure for Evaluating the Drag on Unequal Spheres . ', 43 3.4.1 Evaluation of f(s ,g) for Unequal Spheres. 44 3. 4. 2 Evaluation of Hankel Transforms . . . 51 3.4.3 Quadrature Scheme over n and g . . . 60 3.5 Contribution from Different Regions . . . . 62 3.6 Results for Unequal Spheres . . . . . . . 67 iv Chapter IV. EXPER 4.1 4. 2 ##h 0101-9 V. CONCL 5.1 5.2 APPENDICES BIBLIOGRAPHY IMENTAL PROCEDURE AND RESULTS . Single Sphere Experiments Apparatus 4.2.1 Hall Effect Correction for Newtonian Fluids - . 4. 2. 2 Hall Effect Correction for .Viscoelastic Fluids . . . . . . . . Materials . 4.3.1 Test Fluids . 4. 3. 2 Preparation of Doublets. . Rheological Properties of Test Fluids Experimental Procedure . . Results - . 4.6.1 Previous Observations Non-Contacting Particles 4. 6. 2 Contacting Particles in Newtonian Liquid . 4. 6. 3 Contacting Spheres in Elastic Fluid USION AND RECOMMENDATION Conclusion Recommendation . Page Table LIST OF TABLES 3.1 Condition numbers K of at different values of wwwww 0101th a and S . Change in A, B, C, and D with decreasing a Behavior of f(s,E) with decreasing 6 'Absolute error in f(s,€) as a is decreased Integrand after integrating over s (a=.2, k=5) . Effect of increasing number of quadrature points on integral shown . . . . . . . Drag on pair of contacting spheres--elastic contribu- tions . . . . . . . . . Viscosity, relaxation time, and density of test fluids . Newtonian fluid results Elastic fluid results Comparison between theoretical drag coefficient ratio Xt and experimental drag coefficient ratio Xe . Experimental Raw Data vi Page 45 46 50 50 64 65 68 83 102 108 110 120 Figure wwwN ##9## 01-th h-b-b-P «:0on .10 LIST OF FIGURES Schematic diagram of spheres and coordinates Schematic diagram of the spheres and coordinates . Montonic decay of f(s,E) vs s at a = .5, g = 0.25 Monotonic decay of f(s,g) vs. s for a = 0.5, E = 0.25 . . . . . . Monotonic decay of M vs s for a=0.5, g=0.25 . BE Monotonic decay of $13 vs s for a: .5, 5=o.25 Wall correction factor vs. d/D of Sutterby (1973) Schematic diagram of drop cylinder Photograph of centering device . Photograph of several size doublets Calibration of Weissenberg Rheogoniometer by ASTM fluid . . . . . . . . . . . Corn syrup viscosity vs. shear rate at 25°C Viscosity vs. shear rate for corn syrup at 36°C Viscosity vs. temperature for corn syrup Viscosity vs. shear rate of 0.2% Separan in corn syrup at 25°C Normal stress coefficient vs. shear rate of 0.2% separan in corn syrup . . . \I‘i‘I Page 20 33 48 56 57 58 73 77 78 81 84 85 86 87 88 89 Figure 4.13 First normal stress difference vs. shear rate for 0.2% Separan in corn syrup . Photograph of a doublet in a Newtonian fluid . Photograph of a doublets in a Viscoelastic fluid Photograph of particles separated by a distance in Newtonian fluid . . . . . . . . Photograph of same particles at the bottom of the test cylinder . . . . . . . Photograph of two particles in Viscoelastic fluid at the top of the test cylinder Photograph of same two particles after some dis- tance down the tube . . . . . Photograph of same particles as they converge Photograph of particles separated by large critical distance at the top of the cylinder in a visco- elastic fluid . . . . . . . Photograph of same particles as they diverge down the tube . . . . . . . . . . . . Correction factor for drag coefficient deviation from Newtonian Value due to presence of elasticity vs Ne* . . . . . . . . . Correction factor for drag coefficient deviation from Newtonian values due Unpresence of elasticity vs We* . . . . . viii Page 90 93 94 95 96 97 98 99 100 101 106 107 11> 21> > H O D! m H 11% m U . 112,133 NOMENCLATURE unknown function of s in Equation (3.37) plate area of viscometer in cm dimensionless Rivlin-Ericksen tensors of exterior fluid small sphere radius (cm) coefficient of exponential in f1 in Equation (3.63) unknown function of s in Equation (3.37) coefficient matrix for solving for A(s)--D(s), see Equation (3.54) multiplier of the exponential function f1 function of a in Equation (3.62) unknown function of s in Equation (3.37) unknown column vector in Equation (3.54) drag coefficient elastic drag coefficient Newtonian drag coefficient characteristic particle dimension unknown function of s in Equation (3.37) drop cylinder diameter reference field deformation gradient dimensional deformation gradient dimensionless deformation gradient in Equation (3.29) Sphere diameter ix I'T'I I N '71 VT! '11 HTl-HT6 diameter of viscometer plate unit vector along z-direction drag force vector drag force on contacting particles surface drag force vector due to Newtonian contribution Newtonian drag force first non-Newtonian contribution to the total drag second non-Newtonian contribution to the total drag of order drag force due to elasticity on surface of doublets Newtonian drag in the z-direction twice differentiable function. See Equation (3.34) function of x, n, and 5 given in Equation (3.51) deformation gradient in Equation (1.8), (1.10) function of x, n, and g at small 5 values in the Hankel transforms in Equation (3.66) exponential decaying function in the Hankel transforms 0 to w and 0 to xu given by Equation (3.63) correction factor to the Newtonian fluid for contacting particles absolute error in evaluating f(s,g)--See Table 3.4 hydrodynamic.couple gravitational acceleartion function of s in Equation (3.78) function of g and n in Equation (3.50) transfer function of load cell analytical solution of the Hankel transforms 0 to w See Equations (3.73), (3.74) 11 (2) o ’ Ko Re Hankel transforms given by Equation (3.51) unit tensor complex number = 7:1 integrand of the Hankel transforms at small value of x arbitrary unit vector Bessel function of the first kind and order 0 Bessel function of the first kind and order 1 Bessel function of the first kind and order m wall correction factor raw sum norm condition number in Equation (3.59) wall effect tensor for single sphere settling ratio of large sphere to small sphere radius time constants of Equation (1.8) distance required to attain terminal velocity See Equation (4.5) characteristic container length particle mass fluid mass displaced by particle order of Bessel function first normal stress difference (Newton/meterz) number of simultaneous equations outward unit normal pressure integrand function of E and n of a double integral over a and n in Equation (3.47) Reynolds number r-direction of a cylindrical polar coordinates system xi (I) U) M l1 U) 1 to T l—i l—l -—t —| —1 O l—l [—5 "—1 "—1 H ll—‘I O (‘1' position vector fluid surface far from the particle fluid surface just about the particle particle surface Hankel transform variable. See Equation (3.35) function of a given by Equation (3.68) torque on particle surface function of w in Equation (3.48) torque vector in Equation (1.9) Newtonian torque non-Newtonian torque vector in Equation (1.9) stress tensor of Rivlin-Ericksen type. See Equation 1.3) stress tensor for Viscoelastic contribution stress tensor for reference fluid time (sec) mantissa of a single precision floating point number uniform translation velocity (cm/sec) reference uniform translation velocity (cm/sec) bounded uniform translation velocity uniform translation velocity vector dimensional velocity vector in Equation (1.3) dimensional transpose velocity vector in Equation (1.3) z-direction of a cylindrical polar coordinates system fluid volume surrounding doublets fluid velocity field vector xii IIO') - (9(1), where w(1) is the angular velocity vector .1. 9(1), 9(1) is the velocity vector and its transpose K811) and K82) are two time constants related to the viscosity and normal stress of the fluid. 10 It is by the contribution of these constants to the total hydrodynamic force F and couple G, information was obtained on those preferred orientation. It was shown thatatransversely isotropic particle in a quiescent field will have a terminal orientation in the direction of the external force. In this case, the particle would not rotate once it reaches the terminal state. A particle with its symmetry axis in the plane of the shear will not leave that plane while a particle with its axis parallel to the vorticity axis will always maintain that orientation. For long transversely isotropic particle, Leal (1975) used the slender-body approximation to calculate the hydrodynamic force and torque for simple translation. It was shown that the particle will acquire a terminal orientation that is parallel to the axis of symmetry in the direction of the external force. This was also seen to be the case experimentally. In simple shear flow the results were identical to those of Brunn (1977a). Brunn (1979) investigated the effect of particle shape on the orientation. He considered a near sphere particle and included the particle shape in the analysis. The medium was taken to be represented by the second order fluid. The result of particle sedimentation gave the same conclusion as that of a perfect sphere. This is a terminal orienta- tion in the direction of minimum resistance. The results showed no such agreement for shear flow. In this case, it is shown that the particle migrates in the direction of its axis provided that this is the vorticity axis. In elongational flow, the behavior is quali- tatively the same as in a Newtonian fluid. In a review by Brunn 11 (1980) the motion of rigid particles in Viscoelastic fluid was sur- veyed. The second order fluid model was used to describe the medium. A general formulation for arbitrary rigid particle in a steady motion of negligible inertial effect was considered. A regular perturbation expansion around the Newtonian solution in power of small Neissenberg number was assumed to obtain an expression for the drag force and torque on the particle surface. F - £0 + we E1 + ' I IO + We 11 + . . . (1.9) where the subscript 0 is for Newtonian contribution and the sub- script 1 is for the non-Newtonian (normal stress) contribution. In pure translation the contribution from E1 was obtained via application of the reciprocal theorem which is given as a volume integral around the total fluid volume surrounding the particle. (1.10) I—l o I '11 I--' I N t -——\ Q. < [—1 l-h O o "-11 ’3 L4 "-01 > 0 where i is an arbitrary vector, u is the viscosity are the Newtonian and reference field deformation 0 gradient :0’ ll—h > 12 Brunn (1976) had stressed the difficulty of evaluating the integral of Equation (1.10). He clarified that the difficulty is not due to the volume integral itself, but to the tedious repeated tensorial product :0 . :0: i. The same method of analysis was also outlined by Leal (1975). The merit in using the reciprocal theorem is in that there is no need for obtaining the velocity field in order to obtain the force and torque on the particle surface. Leal also found that a slender or symmetric rod-like particle in a simple shear flow has a correction in the drag force at order one rather than at order two in Neissenberg number as for a sphere. In the next chapter particle interaction is going to be dis- cussed and the problem statement‘is to be presented. Chapter III and IV are devoted to the theoretical and experimental analysis. Finally, Chapter V would present the conclusion and recommendation. CHAPTER II MOTION OF AGGREGATES 2.1 Particle Interaction 2.1.1 Effect of Particle Size Increasing the concentration of particles and shear rate in suspensions leads to the formation of particle doublets due to the hydrodynamic interaction. In general, suspensions of single particles will behave differently from suspension of agglomerates. This is caused by the difference in shape between agglomerates and single particles and by the fluid entrapped in the interstices between the particles. When the particle size is of lum or less, other factors which are known as colloidal forces effect the interaction and sta— bility of doublets in addition to the hydrodynamic force. The first of the colloidal forces is the Brownian. These have been reviewed by Russell (1980). This force is due to random collision between particles due to fluctuation in thermal conditions. The dominance of this force over the electrostatic force is for particle sizes of nano- meters. Electrostatic effect is dominant for particle sizes in the range of .1 to 10pm in diameter. This has been treated by Overbeek (1948). At this size, two electrostatic effects are encountered between the particles; an attractive field as a result of the Van der Naals forces and a repulsive field due to particle 13 14 surface charge. The repulsion effect acts as an opposing force against the attraction. Thus, the superiority of either the attrac- tion or repulsion effect is the limiting factor for forming an aggre- gated or a floctuated system. The stability of the doublets formed was investigated by several workers. Papenhyijzen.(1972) investigated the effect of high and low deformation on aggregated suspension. A network model describing the particle arrangement as a random chain formation was developed. The relative order of magnitude of the effects of hydrodynamic and non-hydrodynamic forces on breaking these chains was calculated for two different suspensions. It was found that the hydrodynamic force is of the same order of magnitude as the non-hydrodynamic. Hoffman (1974) performed a theoretical and experi- mental study of the role of interactive forces on the dilatant vis- cosity behavior in concentrated suspensions of polymer resins in shear flow. The phenomona was explained by the effect of the shear rate on the layer of ordered chains of particles which pass one over another in the direction of flow. It is the disorder of such an arrangement, caused by shear rate increase, that results in dilatancy. A mathematical model was postulated to describe this behavior. Experimental results gave strong evidence to the importance of the repulsive force and shear stress effect. The mathematical model has not been conclusive yet. Zeichner and Schowalter (1977) carried out a similar study of the interparticles forces in a shear flow on the stability of colloidal systems. 15 2.2 Hydrodynamic Effect When the particle size is above those of colloidal size only hydrodynamic force effects exist. Michele et al . (1977) observed that rigid spheres, as well as air bubbles of 60 to 70 pm diameter suspended in visoelastic polymer solutions subjected to laminar shear flow, aligned themselves to form finite chains. They reported also that when two spheres come into contact in such shear flows, no rotation was observed. Riddle et al. (1977) observed pairs of identical rigid spheres (of diameter 0.3 to 0.6 cm) falling along their line of centers in Viscoelastic fluids and found that for initial sparations less than a critical value, the spheres come in contact. All observations indicate the formation of chains of particles in suspensions of both Newtonian and Viscoelastic fluids. The study of hydrodynamic effects on these systems started investi- gating the effect of the interaction between two particles. A theoretical analysis of O'Neill (1969b) studied the slow viscous flow caused by the motion of two equal spheres almost in contact. The spheres were perpendicular to their line of centers. Two cases were considered: one of a translation with uniform equal velocities and the other is of rotation with equal and opposite angular velocities. The drag force for each sphere was obtained for the first problem and the value of each was less than that of a single sphere in the same fluid. The method of solution involved the use of the contacting spheres coordinates. The solution of the creeping flow equation was obtained in terms of several Hankel trans- forms expressed in terms of these coordinates. In the second problem 16 a singular purturbation expansion around the limit of zero separation between the spheres was carried out. In the same year Cooley and O'Neill (1969a) were able to solve the problem of two arbitrary contacting spheres translating slowly in a viscous incompressible fluid. This problem is the same problem we are going to solve in a Viscoelastic fluid. These workers expressed the axisymmetric stream function in the same coordinates of O'Neill, mentioned earlier, as a Hankel transform too. The drag force exerted by the Newtonian fluid on either sphere is less than the drag ona.single sphere, over a range of ratios of sphere diameter. The usame problem was solved later for simple shear flow by Simon and Goren (1971) and Nir and Acrivos (1972). In a Viscoelastic medium, which is the main concern of this work, the study of the hydrodynamic effect on interacting particles has just started. Brunn (1977b) found that equal spherical particles in contact in a second order fluid would yield no correction to Stokes' drag if the solution was considered at order We. For the case where the particles are separated from each other in such a way that the distance between the spheres divided by the spheres radius is much greater than one, the spheres seem to converge and they orient themselves along their center line. No data exist-on the effect of changing size ratio for contacting particles in these fluid. The experimental work that was carried by Riddle et al. (1977) had been for fluids which possess considered shear thinning behavior and so they could not be compared to theories which assume constant shear viscosity as the second order fluid model does. 17 2.3 Bulk Stress The previous factors which influence particle interaction have been studied by a number of investigators in an attempt to deter- mine an effective viscosity for concentrated suspensions of Newtonian medium. Adler (1978) used the cell model to get a concentration dependent effective viscosity. The defect of the model is its dependence on the shape of the cell. Frankel and Acrivos (1967) used the classical hydrodynamic lubrication theory to study the same systems by calculating dissipated energy in the gap between the spherical particles. Both the cell model and the lubrication method suffers from neglecting particle interaction in the analysis. Polymeric suspensions have been investigated experimentally by Highgate and Nhorlew (1970). Three different Viscoelastic systems of various types of rigid spherical particles have been studied. The size of the particle is 100 um in diameter and at concentration of 10% by volume or less. The relative viscosity, defined by ratio of viscosity of suspension to viscosity of suspending fluid, was measured along with the first normal stress difference. This had been done for several solid concentrations within the above range. The results showed that comparing suspension properties to the suspending medium at the same shear stress was a function of concentration only; whereas if the comparison is made at the same shear rate both concentration and shear rate dependence was noticed. The same observations were reported by Kataoka et al. (1978). Until recently, no theoretical explanation was available. 18 Sun and Jayaraman (1982) have derived the bulk stress for sus- pensions of neutrally buoyant spherical particles in a second order fluid medium. They showed that the bulk viscosity of the suspension has a shear thinning factor which is directly related to the elastic- ity of the medium. Their expressions are borne out by the data of Highgate and Nhorlow for systems of concentration up to 7%. These results suggest that if an understanding of systems of higher con- centrations (moderate) is to be achieved, a basic understanding of the role of elasticity in the motion of doublets must be pursued. Thus we exclude colloidal systems and dispersion forces from any future consideration within the scope of this work. 2.4 Objective of Present Research Understanding the behavior of aggregates in flowing polymer liquids is still in a very early stage. This in part is due to the fact that most theoretical analysis has been centered around single particle; while most practical systems are composed of doublets and chains. It is hoped that studying the effect of medium elasticity on the drag experienced by two rigid, contacting spheres in uniform translation would add some light to the subject of filled polymer systems. 2.5 Statement of the Problem Two rigid spheres in contact, one of radius a and another of radius ka with their line of centers along the z-axis in a 19 cylindrical polar coordinates as shown in Figure 2.1 were considered. The pair translates along the z-axis with a constant velocity U in an incompressible vi'scoellastiic. fluid at negligible particle Reynolds numer, Re ~ a Upf/uo where of = fluid density and “0 = zero shear rate viscosity of the fluid. If the fluid relaxation time is small (but finite) compared to the time scale of motion g, the fluid motion will be "rheologically slow" so that the second order stress constitutive equation may be used: V1902 g = ”P9 + 2U [Q “'7rgjf + (91 + 2v2)(Q ' 2)] (2-1) ) where g = is a unit tensor Q = rate of deformation gradient tensor given by g= —%- (v! + vf) (2.2) 11\u_andinm s are the primary and secondary normal stress coefficient. 3? / gzx 5 denotes the corotational derivative given by @Q 32 1 QT=E+{Q-vg} +.2_(_U,.g}-{g-w}) A (2.3) This model will allow us to isolate the elastic effect on the trans- lation of these contacting spheres, since it has a non-shear dependent 20 2'33, "=0 A T Figure 2.1.--Schematic diagram of spheres and coordinates. 21 viscosity p. Thus in the next chapter we will derive an expression for the drag force on the surface of the spheres due to the elasticity of the medium. CHAPTER III DRAG CALCULATION FOR TWO TOUCHING SPHERES The drag on two touching spheres of arbitrary sizes translat- ing in a second order fluid is evaluated in this chapter. First, _the governing equation of motion and the solution procedure are laid down. Next, a volume integral is developed for the elastic effect with the integrand in the most useful form. This volume integral is then evaluated by a sequence of numerical steps. The errors occurring in each numerical step are discussed. 3.1 Equation of Motion The axisymmetrical fluid motion described by the problem statement of Section 2.3 may be represented by the equations: V ° g = O, V o v = O (3.1) where g is the second order stress constitutive equation given earlier by equation (2.1) as: g = 'P§ + ZUEQ ' 2 j t + (V1 + 2V2)(Q ' 2)] C (3.2) The boundary conditions require that Uez on either sphere < ll 9 far from the spheres (3.3) l< II 22 23 Since only uniform translation is considered, in Equation (3.2) the term containing SEQ given by equation (2.3) and the term containing g? Q . Q are significant for this analysis. Thus a modified Neissenberg number, We is defined as: _ U We - (V1 + 2V2) 3' (3.4) while the conventional Heissenberg number is, We = v1U/a. For small values of We (Re << Ne << 1), which is the case in the problem con- sidered here, the velocity field v and the hydrodynamic force on the two spheres 5 may be expressed as a regular perturbation expansion in powers of He. - 2 v - v0 + We v1 + We v2 + . . . 2 F = E + We F + We 52+ . . . (3-5) 0 -1 where v0, f0 denote solutions for a Newtonian fluid, with the same boundary conditions,using only the first two terms of Equation (3.2) and v1, El the correction obtained with the other two terms of Equation (3.2). No attempt is going to be made to evaluate v], since it is possible to solve for E1 without it as shown by Brunn (1980). This is done by applying the reciprocal theorem which was also used by Leal (1975). The background of this theorm, as well as its application to determine the elastic contribution F1 are presented below. 24 3.1.1 Reciprocal Theorem(Background) The reciprocal theorem is a useful device with regard to problems involving the resistance of particles and pressure drops due to fluid moving with respect to particles in creeping flow. Many of the developments and uses of this theorem stem from the work of Lorentz (1906). The theorem can be stated as follows. Let there be a closed surface which is bounding a volume of fluid where we know the velocity and stress fields for a certain steady, incompressible creeping flow in a certain geometry; the theorem says that the force and torque on any surface witin that volume for a different fluid and a different creeping flow but the same geometry may be obtained with- out solving for the velocity and stress fields in the latter situa- tion. The details of this statement can be best explained by show- ing its use for a specific case as presented below. We assume that we know the solution to an incompressible New- tonian fluid in creeping flow for a certain geometry with the equa- tion of motion and continuity. V - £0 = 0 and v - yo = 0 (3.6) where A A A “+ 90 = ' PDQ + 2U(V!o'+ V20) - (3-7) Next we consider an incompressible Viscoelastic fluid within the same geometry as the first fluid and also in creeping flow situation, the relevant equations are 25 Q o I: _ = O and V - y = 0 (3,8) + - H - poé + 2u(Vyo+ Vyo)- P16 + D 2U[2 ' ;% + (V1 + 2V2) 2 ’ 2] (3.9) is the stress tensor of the second order fluid model. Furthermore, let us say that we are interested in getting the contribution at the first perturbation in He for the second fluid such as the force on any surface within the total volume that is enclosing the fluid. Equation (3.8) may be rewritten with 0(Ne) terms as v - (g, + weij1 + Negl) = 0 V - (v0 + Nevl) = 0 (3.10) where 'ld. H II - B1§ + 2pm!1 + vyi) iuo = - Poé + 2u ":1 > Note here that the area integral is evaluated over the entire surface bounding the fluid. This comprises a fluid surface just about the particle, So, as well as a fluid surface S0° far from the particle. Also, the normal 9 is directed away from the fluid at the fluid- 1 particle boundary. Furthermore, as shown in Appendix A [El : v§o - g, : v21] dV = o (3.20) Vf so that we may write = J ;1 : V go dV (3°21) V 29 Note here that the volume integral on the right hand side involves only velocity fields in a Newtonian fluid, v and go. In order to 0 apply the reciprocal theorem in this form successfully, the comple- mentary or known fields (go, go) must be chosen to leave only one unknown on the left hand side, such as the o(Ne) contribution to force on the particle in the other problem. 3.1.2 Application of the Reciprocal Theorem Pursuing the objective outlined in the last section, of leav- ing only one unknown on the left hand side of Equation (3.21), let us choose the complementary problem such that Q0 = 0 far from the particle on Sco ; that is the fluid is quiescent far from the particle. Further, to obtain the z - component of the force contribution at 0(Ne) on the particle, let us choose a uniform translation of the particle along the z — direction for the complementary problem, i.e., l<> o = Oez on Sp the particle surface (3-22) The other problem of interest with unknown velocity and stress fields in a second order fluid is also one of a uniformly translating particle in a quiescent fluid since the velocity is specified on the particle surface v = Uez on Sp (3.23) This is met by the zeroth order term v0 and y1= 0 on Sp (3.24) 30 Finally, the fluid here too is quiescent far from the particle v1 = 0 on Sm (3.25) and we obtain from Equation (3.21) - U F12 = g1 ; v20 dV (3.26) Vf where F12 is the z- component of the 0(We) contribution to force on the particle surface. If we choose, U = U, we obtain F12 = - .2011 (V,+2\),) (20 - go) : V‘lo W (3.27) Vf Of‘ F12 = %E (v. + 2v.) I (2,, - 2,) = 2, W (MB) v _ f Working with nondimensional quantities on the right hand side (length scale a, velocity scale U) 20 the same as go here, * *- = '2U(V1 + 2V2) U2 ( * 0 F dV* (3.29) "O "O O 0 “U 12 O 01" 31 Flz = -2(uUa)We [ ) : Q dV* (3.30) V with nondimensional quantities denoted by an asterisk. This relation given by equation (3.30) is a direct conse- quence of the result of the reciprocal theorem as given by equa- tion (3.21). It is seen that we could avoid evaluating v1 since £1 at order one in Weissenberg number is associated with terms involving 30' At this point we would like to stress that all the work in the remaining part of this chapter is centered around eval- uating this volume integral as given by equation (3.30). 3.2. Use of the Newtonian Solution * The dimensionless, axisymmetric velocity field, Yo needed for determining the volume integral of Equation (3.30) may be obtained from a stream function 0 in cylindrical coordinates as y; = (v*, o, v*); v* = _1_$p_. v* = -'—1- 9p- (3.31) Several investigators have solved for this stream function in vari- ous coordinates. Cooley and O'Neill (1969a) approached this problem using tangent sphere coordinates (mentioned briefly in Chapter I) which are related to cylindrical coordinate by the relations Z*=‘2_§.___;r*=——J_2 ’¢=¢ £2 + n2 £2 + n2 or n +1: = (2*?1”) - 1 = /-‘I (3.32) Figure 3.1 is a schematic diagram of the spheres and the coordinates. In terms of these coordinates, the surfaces of the two spheres are given by g = 1 and g = -a s - l/k. The region occupied by the fluid is given by -a < E < 1 and 0 < n < m. The point of contact of the two spheres is given by n = m; and the region far away from the spheres by E = n = 0. . Cooley and O'Neill showed that the equation of motion for creeping flow of an incompressible Newtonian fluid about the pairs of spheres may be written in terms of the axisymmetric stream function as 4 AW = 0 (3.33) where the operator4/<,is defined in tangent sphere coordinates by ./\E.[(€2 + n2)-& £315. n)] =% (£2 + n2)3/2 [fl- $1.???)- 332] (3.34) 352 ° 802 whereé?(g,n) is any twice differentiable function. The solution to Equation (3.33) waS'Uuyiobtained by them as Hankel transform involving 33 =0 ft‘ztxér fl Figure 3.1. Schematic diagram of the spheres and coordinates. 34 01, a Bessel function of first kind and order 1. oo w = ( 2 n 2)3/2 {(A + EC) sinh 5E + (B + £0) cosh 5E} E + n 0 J1(sn) ds (3.35) where A, B, C, D are functions of 5 found from the no-slip boundary conditions on the moving sphere surfaces,€ = (51 = 1, £2 = -a) which can be written as w - -2n2(€2 + n2)'2 at (E = £13 £2) = 8€n2(€2 + n2)'3 010) (“£3 (3.36) These equations, along with Equation (3.35) yield (as shown in detail by Cooley and O'Neill (1969a) a set of four linear equations in the unknowns A, B, C, and D for arbitrary size Spheres. A sinh $51 + B cosh 5E1 + 51C sinh 551 + E1 0 cosh $51 = ‘Slill - -2. (1511+ s 1) A sinh sgz + B cosh $52 + ngsinh sgz + 52 D cosh $52 = ’SIEZI -1 '29 (lag! + S ) 35 A cosh sgl + B sinh 551 + C(s'1 sinh sgl + £1 cosh 551) + 0(5-1 COSh $61 + £1 sinh sgl) = 2519'51511 1 A cosh $52 + B sinh 552 + C(s' sinh sag-Fazcosh 552) 1 ‘51521 0(5- cosh 5E2 + £2 sinh 5&2) = 2E2 e (3.37) Solution of these equations for the unknowns A, B, C and D furnishes the solution to the stream function 0. The force 50 on the pair of - Spheres has been shown by Copley and O'Neill (1969a) to be given by E0: Zfiuaez sBds . (3.38) It is to be noted that this is obtained by adding the forces on the individual spheres given in Equation (4.3) of their paper. In order to evaluate the elastic contribution to drag, the integrand of Equation (3.30), must be expressed in terms of the con- tacting sphere coordinates. This is because the Newtonian solution of which use is going to be made of here, is given in these coordi- nates, and the volume of fluid around the spheres is most easily expressed in these coordinates. 'Huaprocedure involved performing the ~k * * tensorial operation Q0 oQ : D in terms of cylindrical coordinates, o -0 followed by coordinates transformation with the aid of Equation (3.32) 36 —3 * * * ‘l .15 0 l( Bvr 3 z ) * * * 3r 2 32 Sr * * * Q0 - O Vr/r 0 a * a * * l( v" v2 ) o 3’72,— * 2 32 ar* 32 * 'k and the dot product Q0 - 90 yield _ 3 * 3v* v* v* 3v* 3v* 1' v z <—£—)2+1/4(—4;+—é>2 o -—";< 3: + .) Br 32 3r Zr 32 Gr * 1: * Vr QOOQO- 0 :5 0 * a* a* 3* 3* -v v v v v 2 1 -E($+ 3.) 0 (—§- +z(—:+ Zr 32 ar 32 az * 3V2 2 "1; ) 3r (3.39) The above quantities yield for 37 v*2 3 * ** 3 'k 'k * * * v v v 2 v 2 ° 9 = Q = - 3E " -—%é—-+-5;—-( -—%é—) + —5ir ° 0 o I;;3.32 r 32 4r 3 * a * v * v H:— + -i—>21 (3.40) 32 Sr The velocity components in Equation (3.40) can be expressed in terms of the stream function using Equations (3.31) and some rearrangement to give. 2 2 * * * 3 2 20-20 Qo=-*4[-33. (a—i—afl’r) +317(§“’:)2;-2—3-4’: Y‘ Z Z l‘ I‘ Z Z l‘ 2 * +r 31.3122 3* <-1/r*i“’7)1 (3.41) 32 32 Br Br At this point coordinate transformation is to be performed. To start with, the volume element of Equation (3.30) in tangent sphere coordinates can be written as: dV* = dgdnd¢ = Bdgdnda, (3.42) "a“nh¢ (n2 + a2)5/2 where hg, hn’ h¢ are scale factors which may be obtained with the procedure discussed by Happel and Brenner (1965). .The trans- formation of Equation (3.41) into the contacting spheres coordinate system is a rather lengthy but straight forward process. It involves handling a large number of terms which resulted from repeated use of the chain rule as required for each expression 38 in the above equation. To clarify the previous statement further, 3!, is going to be transformed below. With w = w(n, g) given 32 aw ,. 34 35 9;) an 3—1' ‘3: _T+3n ‘1' (3-43) 2 32 32 The partial derivativeség;, 91*- and others like in? SZ 82 3r and §§;-can be obtained from Equations (3.32) which yield, 3r 2 2 3 3 - ._D.;=.. gn,_§?=n__2£_ 32 32 2 2 Ln; .-. - (n 25 ), 3g = - an (3.44) 3r 3r With the aid of Equation (3.44) the individual terms in Equations (3.41) can be written in terms of the contacting coordinates as: :3 _(n2_€2) 34_ am 31W — 33% ‘ ‘—-?r- 53- US 2 2_ 2 31111) 3 W = LII—2L2 33%(3119) ' NE 1(314’) 32* 312W (31¢) Ill 0) N G- I I m :5 file) A O.) H '6 V + m N I J . N O.) 39 = * 3 = 2n _ ii_ 3221(1- 7' 8:; (32¢) _—n2.+§2 [5” 36, (32‘?) + gzznz %(32w):1 (3.45) Now we can express Equation (3.30) in terms of contacting spheres by utilizing the results,(fi’ Equation (3.42) and (3.45) which yield F12 _ 1 EEU" -324 dz dn Q(€,n) (3.4s) -, 0 where 2 23/2 Q(€.n) = ilL—{531—l [11(4) + 12(4) + T3(w)1 (3.47) 11(4) =-§% (alw)925. 3 _35,’ 1 II 0 31 (1.] _§ -_- 3f .1 dx 3%; F? (x dx) 7" O (3.51) 42 -It is to be emphasized that I1 through I6 have to be obtained before any other numerical calculation can be carried out. 3.3 Results for Equal Spheres For equal spheres, a = 1 and the boundary conditions given by Equation (3.36) lead to nonzero values only for B and C in Equations (3.37) so the stream function 0 is an even function of g l .- _ n V (52 + n2)3/2 ): {EC sinh $5 + B cosh sElJ1 (sn)ds Furthermore, Cooley and O'Neill have noted the explicit expressions ' -[2 + 25 + s-1(1-e'25)]/[s + sinh s cosh s] w I [1 + Zs-e'ZSJ/[s + sinh s cosh s] (3.52) ('5 II This feature is useful here in obtaining an analytical answer for the integral in Equation (3.46). A quick look at the individual terms of Equations (3.45) in terms of the even stream function would show that 310 = even 311w = odd 312w = even _ 82211) = Odd (3.53) and so it is readily seen from Equation (3.48) that for this case, T1(w), T2(W), T3(w)-—all turn out to be odd functions of g. Inte- gration over a from -1 to +1 should yield zero. Hence, the drag 43 exerted by a second order fluid on a pair of identical, touching Spheres is the same up to 0(We) as the drag exerted by a Newtonian fluid with the same viscosity. This result will be used later in the case of arbitrary spheres to check the numerical procedures developed to evaluate the integral in Equation (3.46). It is worth noting here thata pair of identical, touching spheres is a body of revolution with fore-after symmetry and thus belongs to the class of transveresly isotropic particles. For such particles Brunn (1977a) has shown that the 0(We) contribution to the drag is zero. 3.4 Numerical Procedure for Evaluating the 0rag,on Unequél Spheres In this section we are going to have three subsections. Section 3.4.1 is devoted to evaluating the function f(s,g) for unequal spheres; Section 3.4.2 is the detailed evaluation of the Hankel transforms; and Section 3.4.3 to discussing the quadrature scheme over n and g. The importance of the first two subsections is because for uneuqal spheres (a f 1), the functions A, B, C, and D of 5 (see Equation (3.37)) must be evaluated numerically; so the integral of Equation (3.46) has to be evaluated numerically; managing such an expression is no trivial matter. This is not because of the triple integration that had to be carried out, but rather because of the great need of very accurate numerical evaluations of the function f(S,€) and subsequently the Hankel transforms. 44 3.4.1 Evaluation of f(s,€) for Unequal Spheres The function f(s,E) is evaluated thousands of times in the product quadrature scheme used for the multiple integration of Equation (3.46). Hence, the accurate evaluation of A(s), 8(5), 0(5), and 0(5) is at the heart of the lengthy sequence of the numeri- cal steps in this work. These functions of s are obtained as the solution to the four linear equations (3.37). Detailed error esti- mates for different valuse of S and a may be obtained by writing (3.37) as «299 = 9 (3.54) where 9T = [A(s), 3(5), c1 (3.55) Defining the vector norm of c by ||§||0° = max |ci| (3.56) 1 . and the matrix (row sum) norm of gby 1. .|0° = max (2 | ) (3.57) 24 J. , L4 _ we obtain an upper bound on the relative error Ilégll/llcll according to Goult et al. (1974), where R is a condition number defined by 1841,0810- 45 llésll/llcll _<_ Rowe-52‘? (3.58) (3.59) n is the number of Equations (4) and t, the number of bits in the mantissa of a Single precision floating point number is 48 on the cyber 750 at Michigan State University. 10'5 say, K may be as large as 5 x 10 . 7 So for a relative error of Table 3.1 shows the condition numbers of the matrngiB for different values of s at several values of a between .05 and 5. a marked increase at both very low and large values of s. It is to be noticed that the values of K show Furthermore, TABLE 3.1.--Condition numbers K ofgat different values of 01 and S. S a 0.1 1.0 6.0 10.0 20.0 0.05 0.1217x106 0.9038x102 0.4732x104 0.3226x106 0.805 x1010 0.1 0.1062x106 0.8075x102 0.3650x104 10.2042x106 0.3095x1010 0.2 0.8221x105 0.6595x102 0.2164x104 0.8139x105 0.4554x109 0.5 0.4281x105 0.4184x102 0.4469x103 0.5004x104 0.1400x107 1 0.1864x104 0.2822x102 0.5633x102 0.882 x102 0.1681x103 5 0.1073x104 0.1813x104 0.3501x1013 0.4992x1020 9.2283x1038 46 the actual magnitude of A, B, C, and D are also larger with decreas- ing a at large values of s as can be seen in Table 3.2. For TABLE 3.2.--Change in A, B, C, and D with decreasing a), s=10 a=.5 s=10 a=.2 A 5.538 x 10'4 4.762 x 10"2 B -5.540 x 10‘4 -4.762 x 10'2 c 9.989 x 10'4 0.1832 0 -9.987 x 10‘4 -0.1832 both the small values and large values of s, Cooley and O'Neill (1969a) have asymptotic estimates for A, B, C, and D. = -2(a-1)(a2 + 4a +1) (a+1)3s for s_: 0.1 (3.60) and for large S. A, B, C, and 0 all are 0(56-285) (3.61) where B = min {1, 0} Thus in our evaluationcfiithe function f(s,g) to overcome the error in evaluating A, B, C, and D at both small s 5_.1 and large 5 3_10 we do the following. First at small S the expressions of Equations 47 (3.60) are used and the function f(s.€) at this range of s < .1 will be referred to as fO from now on. On the other hand for large 5 we obtain with Equation (3.61) for A, B, C, D, the following asymptotic expression for f(s,g) f - se'5(c"E ) (b2 + b3g) (3.62) where b2 and b3 are functions of a. A plot of f(s.§) versus 5 shows a monotonic decay for the function above S = 7.5 as seen in Figure 3.2. It was seen possible to fit f(s,E) with a Single exponential function f = b1 9-915 (3.63) 1 where b1, a1 are obtained from fitting f(S,g) keeping 0,5 fixed. More will be said about the fitting procedure in the next section. For convenience, the form of Equation (3.63) is used for values of 5.: 10. Thus for 5.: 10 the function f(S,E) is referred to as f1. Another look at Table 3.1-would show that the relative error in most cases is less than 10'7; but for a = 5 the error at larger values of 5 becomes enormous. This large error is typical for a greater than 1. However, of all possible values for a, it is enough to consider the range of 0 < d < 1, because the results for the remaining values may be found by reversing the sign of the right-hand side of Equation (3.37). The last point that we would like to stress is the magnitude of error in calculating f(s,g) as a is decreased. In Table 3.2 we 48 (It? 7%“) -/.Cl - 2.0 - /.0 6.0 j /0.0 5 Figure 3.2. Monotonic decay of f(s,€) VS S at a = .5, E = 0-25 49 have shown that the values of A, B, C, and 0 increase in magnitude_ as a is lowered from .5 to .2 at S = 10. This increase in magnitude has also caused the value of f(S,g) to increase markedly for this . change in a at s = 10. Thus, the accuracy of evaluating f(s,g) is decreased for values of a less than 0.5. The behavior of f(s,g) as a is decreased is given in Table 3.3. It is to be noticed that the values of f(s,§) increase by one order of magnitude as a is decreased from .5 to .2 at large values of 5 above s = 1. Finally in Table 3.4 we present the absolute error lAfl in calculating f(s,g). |Af| is obtained by finding the error in A, B, C, and D, and multiplying by f(s,€). The table shows us the change in the magnitude of the absolute error as a is decreased._ It is seen that as a is decreased, the absolute error is higher in two particular situations. The first is that for all values of g as a is decreased, the error is largest 10. The second situation is seen to be associated with two at 5 distinct regions on the surface of the contacting spheres. These are at the stagnation point on the larger sphere at g = -a and in the region far away from the Spheres at E = 0. These absolute errors would accumulate each time the function f(s, E) is evaluated with the numerical scheme. In that scheme the function f(s,E) is evaluated 27,000 times on the average. Multiplying this number by Af would give an upper bound on the total error involved in evaluating f(s, E). Using the maximum value of |8f| over 4 s and a, we obtain an upper bound of 10' on total error due to function evaluation at a = 0.2. TABLE 3.3.--Behavior of f(s,€) with decreasing a 50 a=.5 g=-.5 a==.5 g=0.0 a=.5 g=1.0 s f(5.€) f(S.€) f(S.€l 1 -1.8196 -1.9108 -1.471 6 -0.066383 -0.020822 - .00578 10 -0.008086 -0.000554 - .00009988 a= 2 g=- 2 a=.2 g=o.0 a=.2 g=1.0 1 -1.965 -1.9881 -1.462 6 - .22088 - .18998 - .00578 10 - .0812 - .047621 - .00009988 TABLE 3.4.--Absolute error in f(s,€) as a is decreased 8:.5 g=-.5 a=.5 2=0.0 a=.5 g=1.0 S IAII. lAfl lAfl 1 1.5 x 10'11 1.6 x 10‘11 1.21 x 10'11 6 5 9 x 10‘12 1.85x10‘12 5.13 x 10'13 10 8 1 x 10'12 5.55;<10’13 1.0 x 10"13 a =.2 g=-.2 6:.2 5:0.0 6:.2 g=1.0 1 2.5 x 10’11 2.63 x 10‘11 1:95 x 10'11 6 9.5 x 10‘11 8.17 x 10‘11 2.5 x 10'12 10 1 3 x 10‘9 7.62 x 10'10 1.59 x 10‘12 51 ‘§,4.2 Evaluation of Hankel Transforms There are several problems to be addressed as we proceed in describing the procedure of evaluating the Hankel transforms I1 . . . I6 given in Equation (3.51). For convenience, only one Hankel transform is going to be used to describe the method of solution. The others are evaluated in a similar fashion. Also, we will use the variable x = Sn which is the argument of the Bessel function. In particular, we choose the transform 11 =-%-I:f(x/066)J1(x)dx . (3.64) As we discussed in Section 3.4.1 that asymptotic expressions f0 and f1 will be used over a range of small 5 and a range of large 3 values repsectively. Hence "2 xu f(x/n,g)Jl(x)dx = fodl(x)dx + le(x)dx O 0 xx. + f101(x)dx (3-55) xu We note that in Equation (3.65) XH= Sn and as we mentioned earlier it is more convenient for the behavior of the Bessel functions to use x rather than Sn. This is primarily to keep track of the number 52 of cycles over the limits of integration and to provide enough quad- rature points in each cycle. The quadrature points which we are referring to here are those of a Gauss Chebyshev quadrature scheme which is used for the numerical integration as will be discussed shortly. In Equation (3.65) xz = .1n; however, x2 should be less than or equal to 0.3 so that the power series representation of the Bessel function (to be discussed below) is valid. Thus, x = min [.ln,.3]; 9. and Xu = 10n as established earlier in Section 3.4.1. In the sub- sequent paragraphs each of those integrals on theright hand side of Equation C365) is going to be discussed separately. In the first integral on the right hand side of Equation (3.65) x9 0 f0 Jl(x)dx (3.66) fo has been defined in Section 3.4.1 and because of the unbounded functions A(s), B(S) as s + 0 (see Equation (3.60)), for 01(x), we considered the first few terms of a power series representation of the Bessel function 3 _x__ x<'O.3 ~ X J1m" 2‘ 16 — When this is done, the integrand of Equation (3.66) iS no more unbounded and can be written as: 2 11 = f001(x) = (Tn/2)sinh(£n§) - (Di-63) sinh (X—f) 2 - ncosh (é?) + (Lgn) cosh (é?) (3.67) 53 T = -2 (01-1)(012l+4a + 1) where 3 (a + 1) (3.68) The integration of Equation (3.66) was carried out numerically using Gauss-Chebyshev quadrature scheme. The integrands for the rest of the Hankel transforms I2 . . . 16 for this limit 04+ xi were obtained by a Similar procedure and are listed below. 2 (1 - %$r)((T/2)sinh(é§) - cosh (é§)) 1'2 = - - 1i , . 25. 13 - n( )((T/2)51nh(XE/n) - cosh( )) 4 n2 9 i4 = (Tcosh(xE/n) + 2Sinh 0%)) 01(x) 15 =%(1 cosh(’%)—2 sinh (%))(% Jo(x)-Jl(x)) i6 = (T—rf s1 Ming) 4% cosh(£n§))dl(x) (3.69) The second integral on the right hand side of Equation (3.65) over the range x2 to x was obtained with the Gauss Chebyshev u quadrature scheme. Here the function f is the actual function as given with the original Hankel transform in Equation (3.51). A(s), B(S), C(s), and B(S) are evaluated by the Gauss-elimination routine from Equation (3.37), and no truncation is used for the Bessel func- tion. The number of quadrature points supplied was based on deter- mining the number of cycles for that range of integration, and at least three points per cycle were supplied. 54 Finally the third integral on the right hand side of Equation (3.65) is to be discussed. In Section 3.4.1 we mentioned that the errors in evaluating A(S) -- D(s) is also large at large 5. We also mentioned that f1 = ble‘alS (3.70) This form has been obtained by observing the behavior of f(s,E) at large values of s. Figure 3.3 shows this behavior where a monotonic decay iS observed after about 5 = %-= 7.5. It should be mentioned that for the otherHankel transforms I2 . . . 16, their function of s are born by derivatives of f(S,E) with respect to E as can be seen in Equation (3.51). The behavior of these functions at larges is shown in Figures 3.4 and 3.5 where a similar monotonic decay is observed at about the same value of S =-%-= 7.5. Thus, fitting each of f(s,E), i3féélél-andfig-g-Eiflby an exponential near s = 10 would pre- vent the error in evaluating A(S)--B(s) at larger value of s. In f1, b1 and a1 are determined by the fitting procedure of fig; 27:1;- for each set of a and E separately. For each a, several values of E between -a and 1 were picked and a range of s values (between 7 and 12),that was seen to yield a uniform decay, was chosen. Even though we have been able to overcome the problem of evaluating those functions of s at large s by the exponential fit, we have noticed that for a less than 0.33 there was a considerable difference between the values of the functions f(s,EL-§££§&§l 9 3E 2 I I I . 9.;éggél and their exponential fit f1, f1. f1 . The difference 55 I 2 is larger and more significant between 3fg§,E)’ 3 géi'g ) and f f1 . 3f 36 10 2, while the error with f is of order 10'6. However, in the case 19 For example, when a = 0. 2, the fitting error associated with 2 is of order 10 4 and the error associated with $33-1s of order of a = 0.5, the fitting errors are much smaller. Thus, once more we see that this range of a values would have one additional error to that discussed in Section 3.4.1 in evaluating the function because of the fitting procedure. It should be stressed that this error especially in fitting Efééiéi-and 322355 is going to be large as we repeat their evaluation in the numerical scheme for these a values less than 0.33. The Hankel transforms of a decaying exponential is available as DI H ‘ 4. ~ .23 N . I e‘als 01(sn)ds = 1 (3.71) 91 I—' + J N Hence, we may write the last integral of Equation (3.65) as xu=10n f1 J1 (x)dx = HTl- f1J1(x)dx J Xu 0 W— XU=10n =b a1 + nZ-a1_ Talx 1 Na: + 02 ble—f‘] Jl(x)dx (3.72) 1 0 And the expressions for the remaining Hankel transforms in the limit 0 +'w are 56 0.0 775) '-/.0 - 2.0 — /.o 5.0 f /0.0 5 Figure 3.3. Monotonic decay of f(s,§) vs. 5 for a = 0.5, g = 0.25. 57 3f(s,€) 35 .17 .— .16 "' .14— .12 .08 __ .06 .04 -- 7 7.5 8 8.5 9 9.5 10.0 Figure 3.4. Monotonic decay of _8f_g_2_,€_l vs. 5 for 0L=0.5, £70.25- 58 32f(s,€) agz . -.07 '— -.04 — -.10 l l l I L I 6.5 7 7.5 8 8.5 9 9-5 Figure 3.5.--Monotonic decay of 3;- vs s for OL = .5, €= 0.25. - a + 2a n - (a2 + 2)3/2 H12 - aggl = b1 ( 1 21 2 313/2 ) (3.73) n (al + n ) H11 2 a1 a1 H13 - .,.. - b ( ¢ (1- . . \- 1 :3 w. .1' .152. .9312 Baln - ) (3.74) (ng + a§)5/2 It is to be noticed that HT4 and HT6 are given by expression similar to that for HT1 in Equation (3.72), but with different coefficients. Their coefficients are a4, b4, and a6, b6 respectively and obtained from fitting-gg-and-ggg. H15 is similar to H12 with as and b5 obtained from fitting g2. The indices associated with HT's are for the purpose of identifying the specific Hankel transform I1 . . . 16 in Equation (3.51). The other integral on the right hand side of Equation (3.72) is (x)dx (3.75) This integral was obtained numerically using the same quadrature scheme. Now that we have discussed the individual integrals and explained the motivation behind that procedure, we rewrite Equation (3.65) 60 I1 = rf(x/n,€(dl(><)dx .0 x1 = .In or .3 2 2 = [TM-é- - §g)sinh(’-‘n§) + 11(58- -1)cosh(5n§)]dx ‘0 xu=10n t J f(X/n.€)J1(X)dX x =.lnor .3 - x =10n VU 'X‘ ‘1“? + “2 " at1 "3‘1"; b ( ) - b e J1(x)dx (3.76) 1 nJai +112 0 This form is the one that was used in the numerical evaluation. An + additional check on Unanumerical values of these Hankel transforms is available from the work of Soni and Soni (1973) on asymptotic esti- mates at large n. These estimates are discussed in the next section. 3.4.3 ,Quadrature Scheme over n and g It was mentioned earlier that in terms of the contacting spheres coordinates employed here, the region at which n = w is the point of contact of the spheres and the region of n = g ='0 is the region far from the spheres. The contribution to the elastic effect at the contact point was investigated through the use of asymptotic estimates of the six Hankel transforms given in Equation (3.51) as n + w and by writing theintegral over n as follows. 61 m 10 [0(n.€.11 . . -15)dn = 01mg. Il . . ~15) dn 0 0 + Q(n,g, asymptotic estimates of 1's) dn (3.77) 10 where the asymptotic estimates of the Hankel transforms were obtained by a theorem of Soni and Soni (1973). Soni and Soni have related the behavior of a function of s that is unbounded at‘s + 0 of the form s'Yg(s) where y > 0 and 0 < y < m + 3/2 with the limit of its Hankel transform at n‘+ m. lim "Y n+m /sn Jm(sn) s g(s)ds 0 -Y + "2 iifl‘ééilfliifil H 9w”) (318) where the function g(s) is bounded over the entire range of s. It is to be noticed that in Equation (3.78) si'Y g(s) is equal to the specific function of s in the various Hankel transforms. Further the form of equation (3.67) has n% which is not present in our defini- tion of the Hankel transforms, so we have to multiply the results of Equation (3.73) by n'%. The coefficient y is chosen in such a way as to make f(s,g) bounded at s = 0, m is the order of Bessel function. Application of Equation (3.78) to the six Hankel transforms showed that as n +'w the following estimates are obtained: 62 11 ' -2 12 = 0 I3 = 0 I4 - 0 15 = 0 16 = 0 (3.79) In the last integral of Equation (3.77) when infinity was replaced by 15 or 20 and using the results of Equation (3.79) along with other functions of g and n yielded values of the order of 10'5 which were practically zero. This result had established two things. First, the point of contact of the spheres yield no contribution to whatever result we get for the elastic effect on the drag. Second the value of n = 10 is a reasonable upper limt for n. 0n the other hand, at low n the value of the integrand is large, particularly at E §_0. It is useful to note here that the range of low n with a near 0 describes the region far away from the spheres while the range of low n with 5 near -a describes the stag- nation points region about the larger sphere. More will be said in Section 3.6 about the contributions from different regions to the elastic effect on the drag. 3.5 Contribution from Different Regions We proceed to look at the integrand Q(g,n) of the double integral over n and g in Equation (3.46). This integrand involves large sequence of algebric expressions. It is useful to look at the 63 behavior of Q(g,n) in different regions of fluid around the spheres. In Section 3.4.3 we have already shown that the region near the equatorial section of the spheres and in particular the point of con- tact has no contribution to the drag. This region is that of n + infinity. 0n the other hand, in the regions of low values of n and nonzero values of g we observed a different behavior. The numerical values of the integrand Q(g,n) for a = .2 are listed in Table 3.5. The values of n and g are those at which the integrand started to have an appreciable magnitude. Let us first look at low values of n, with g approaching zero, which represent the region far away from spheres as can be seen from The values of Q(g,n) are significant in this region as shown in Table 3.5 Next we looked at the region near E_:_-a and at low n which describe the stagnation points on the surface of the rear sphere. It was seen that this region contribute the most to the elastic effect on the drag. The range of low n values over which the integrand Q(g,n) takes on appreciable magnitude is larger around the larger sphere (E < 0) than around the small sphere (g-> 0). The large contribution from the stagnation point region is understandable in the light of previous work. Of particular importance to us is the work of Leal (1975) where he treated the 64 Table 3.5.--Integrand after integrating over x (a=.2, k=5) E n Q(€.n) 0.79 0.007 0.29 0.713 0.007 0.29 0.63 0.007 0.26 0.54 0.007 0.23 0.447 0.007 -0.21 0.087 0 007 0.42 0.01 0.062 -O.60 0.01 0.007 -950.7 -0.056 0.062 253.8 -0.056 0.007 -2.08 ~0.112 0.332 -1.01 -0.112 0.17 5.74 -0.112 0 062 165.8 -0.112 0.007 2.71 -0.154 0.17 18.38 -0 154 0.062 93.83 -0 154 0.007 -641.3 40.183 0.062 64.74 -O.183 0 007 -1.31 -0.198 0.062 1.56 -0.198 0.007 0.21 65 case of a long, slender rod-like cone without fore-aft symmetry. In that work, the drag force due to the elasticity of the medium was represented by two integrals. The first is a surface integral which is evaluated at the surface of the particle using asymptotic expres- sion for the integrand. These expressions are not valid for the region close to the ends of the particle. Further in the volume integral it was seen that the region of fluid close to the particle is dominant; and because of the slender body approximation the stagna- tion points at both ends of the particles where omitted. However, the fact that we have an appreciable contribution even near a = 0 (i.e., far from the spheres) is at variance with the result of Leal. It is to be mentioned that the shape of our particles is different from Leal's and this may explain this variance in the results for this region. We would like to conclude this section by presenting an estimate of the error associated with the integration over g and n. It was just seen that the region of small n and E < 0 has the largest contribution to the drag on the spheres surfaces. Thus, when we evaluated the numerical value of the integral of Equation (3.46) for equal spheres (a = 1) where a high degree of accuracy was available in the Hankel transforms evaluation, we obtained a final answer of 0.12. This value was obtained by integrating over a from -1 to 0 and over n from 0 to 1 which are the regions of major contribution to the drag. The integration from-*1to 0 over a and 1 to 10 over n 4 has a value of 10' which is practically zero. Before we can say anything about the value 0.12 for the integral of equal spheres, we 66 further investigated the effect of changing the number of quadrature points over the range -a to 0 over E and 0 to 1 over n for several 0 values. As can be seen in Table 3.6 that the effect of increasing TABLE 3.6.--Effect of increasing number oquuadrature points on integral shown 30 points each 50 Points each . 1°25 I: (15“ l: d“ 0.33 0.98 0.95 0.50 0.48 0.46 0.67 0.28 0.27 the number from 30 to 50 has an effect on the second decimal point of the integral value. The other‘integral -01 to 0 over a and 1 to 10 ‘4). Thus, we conclude that even over n is also very small 0 (10 with this increase in the number of quadrature points, the first decimal point is unchanged. We also know that in the integration over 8 and n the integrand does not involve 6. Thus, we say that since for equal spheres the result of Equation (3.46) should be zero, there- fore, the value of 0.12 is error due to the quadrature over 5 and n. This magnitude is also associated with other values of a's # 1 which we will present their final results in the next section. 67 3.6 Results for Unequal Spheres Before we present the results of the drag correction factor Flz/GHuUa due to the elastic medium for unequal spheres (a f 1), we would like to recall attention to Sections 3.4.1 and 3.4.2. In Section 3.4.1 we aimed at establishing confidence in the numerical evaluation of the function f(s,g) of the Hankel transform. We have seen that as a is decreased from 0.5 to 0.2 the absolute error lAfl increased by an order of magnitude. Furthermore, we showed that as we repeat the evaluation of f(s,£) the error adds up to reach an upper bound of 10"4 for a = .2. In Section 3.5.2 we discussed the error involved in fitting éflélél-and Ejjlgiél-at a less than 0.33 and said that the BE BE order of magnitude of the difference between these functions and the fitting function f1 is large. We also mentioned that the difference will add up as we repeat evaluating these functions. Thus, in light of these points, it was decided that values of 8': .33 where the upper bound of the total error in evaluating all the functions of s and also the difference between the functions and 8 4 their fitting values is of order 10' to 10' , is an acceptable degree of accuracy. In Table 3.7 we present the result for unequal spheres for k =-§-between 1.5 and 3 or a between .33 and .67. The Foz , _ 6nuUa the correspond1ng Newton1an drag correction factor for those particle size ratios. These results F - 12 table includes both wuua and show that there is an appreciable elastic effect on the drag at order one in Weissenberg number. Also, it is seen that this effect increases as k is increased. 68 TABLE 3.7.--Drag on pair of contacting spheres--elastic contributions K -FOE/600Ua -F12/610Ua 1 1.29 .64 1.5 1.66 1.23 2 2.09 1.89 3 3.04 3.42 Note: F2 = FOZ + We F12 Before we conclude this chapter we would like to recall that the problem treated here is that of two contacting spherical particles translating along their center line in the positive z-direction with the small sphere leading. He mentioned earlier in Section 3.5.1 that the case where the spheres translates in the negative 2- direction with the larger sphere leading can be considered by changing the signs of the right hand sides of Equation (3.37). The sign of the quantity Flz is not going to be affected by this sign change of Equation (3.37). since it is associated with the quadratic quantity QO- 90' .FOF the case of spheres translating in the positive 2 direction FOz acts in the negative z-direction; the net drag F = F 2 + WeF 0 12 If we consider We = .1 and for a = 1.5 F = -6nuUa(l.66) - .60uU6(1.233) 69 which is greater than the Newtonian value. 0n the other hand, the case where the large sphere is leading in the negative z-direction, Foz acts up in the positive z-direction. Thus F = 6nuUa(1.66) - .6nuUa(1.233) which is less than the Newtonian value. In Chapter IV we shall dis— cuss the results of Table 3.7 in light of the experimental results. CHAPTER IV EXPERIMENTAL PROCEDURE AND RESULTS 4.1 Single Sphere Experiments Sedimentation is a common method for measuring the drag coef- ficient of particles translating in a fluid. It is very accurate when care is exercised in both the design and procedure. Several investigators--Sutterby (1973), Sakai et al. (1977/78), Sigli and Coutanceau (1977), Chhabra et al. (1980), Broad and Mena (1974), and Acharya et al. (1976)--have performed such experiments using single- rigid spheres with both Newtonian and Viscoelastic fluids. The only measurements needed are settling velocity, particle size, and density and the fluid density. In terms of isolating the effect of fluid elasticity on drag, Boger and coworkers (1980) have obtained the most accurate results using so-called Boger fluids. Most of these experiments were done at very low particle Reynolds numbers. The work of Sigli and Coutanceau (1977) has addressed inertial effects. It was concluded that inertial effects generally opposed the elastic effect. This means that as the Reynolds number is increased, the elastic effect is decreased. Creeping flow conditions are commonly observed since most of the theoretical work applies only to such flow conditions. The purpose of the present work is to investigate the effect of elasticity on the drag experienced by rigid contacting spheres 70 71 translating in a non-shear thinning fluid along their line of centers. This is to be performed over a range of Weissenberg number and at low Reynolds number. The remainder of this chapter will include apparatus design, material properties and preparation, experimental procedure, and finally, results. 4.2 Apparatus 4.2.1 Hall Effect Correction for Néwtonian Fluids A drop cylinder made of Plexiglass was designed for the experi- ment. In this design, the following factors have been taken into consideration, which are very critical from the standpoint of the degree of accuracy. In selecting the cylinder diameter, the wall effect is the major concern. For a single sphere in a Newtonian fluid, Faxen (1932) has made a theoretical analysis of the correc- tion to Stokes law due to the presence of a boundary. The resulting expression for drag is F = GnuUaK (4.1) where 1 K: 1-2.104(d/D) + 2.09(d/0)3 - .95(d/0)5 d/D is the ratio of sphere to tube diameters The above expression was verified experimentally by Bacon (1936) for d/D up to 0.32. The linearity of the constitutive relation of Newtonian fluids enabled Brenner (1964) to derive expressions for the 72 first order effect of wall proximity as a correction formula for the terminal velocity given as _'E g = 0” + . + 0(C0/2) (4.2) 6002 and a correction force formula obtained by a simple inversion of the above expression. In Equation (4.2), g is the velocity of the particle when settling in an infinite medium, 0” is the velocity when the particle is settling under the influence of an outside force 5 at a distance from a boundary whose wall effect tensor is 5 and C0 is a characteristic particle dimension, usually the radius. The second order tensor 5 was obtained by the method of reflection. Sutterby (1973) studied the wall and inertial effect experimentally overa range of d/D between 0.0025 to 0.125 for a range of Reynolds number between 0.00001 to 3.78. The falling sphere data were corre- lated as a relationship between 05/0 5 K, d/D and Re = pUd/u. Here “5 is the fluid viscosity obtained from Stokes law and u is corrected fluid viscosity. The value of K for several values Re was given in a graphical form and reproduced here as Figure 4.1. Agreement with Faxen results is up to Re = 0.2; where beyond this value inertial effect is appreciable. 4.2.2 Hall Effect Correction for Viscoelastic Fluids In the case of a Viscoelastic material, the wall effect is not fully determined. Unlike the case of viscous flow, the non- linearities of the fluid's constitutive relations for a Viscoelastic 73 0-5 0-4 l k\\\\ 0.3 2.5 / N W C12 15 ,rr”’/A 1.9/ 0-8 , 1 0! 0-6 ’ o x’ o. FAxéN (5) FOR Re-O.2 o 00 e - d/g C1CXD cycxs (DJC) Figure 4.1. Nall correction factor vs. d/0 of Sutterby (1973)- 74 medium doesn't in general allow the wall effect to be expressed as a force correction formula. Only a velocity correction factor is possible. Caswell (1970) presented the expression for the first order effect of the wall proximity on particle settling as: g = Uco + 5 - f/6nu02 + 0(2'2) (4.3) The only restriction on the constitutive equation for the validity of the above relation is that it must describe an isotropic fluid which has a lower Newtonian regime with zero shear viscosity. This general relation was examined for the case where the stress expres- sion for the medium is represented by the third order fluid model. Translation induced by a force alone and rotation induced by torque alone was considered. The solution involved velocity perturbation expansions and Green's function method was utilized for obtaining the wall effect tensor 5. The final expression for the unbounded velocity was expressed in terms of the zero shear viscosity for a sphere settling in a Viscoelastic medium as: 6naU F F 'E “0 U03 6na2 611a2 where A is a combination of material constants given in Caswell (1970) Various experiments were carried out to estimate the critical ratio of particle to tube diameter above which wall effect is significant. Sigli (1977) observed that for a ratio greater than 0.25 the wall 75 proximity increased the effect of fluid elasticity. This effect is seen by a decrease in the particle velocity. Boger and coworkers (1980) had established experimental conditions for negligible wall effects. They investigated situations for d/D between 0.04 and 0.2. It was found that the terminal velocity decreased as the sphere to tube diameter ratio increased, but this reduction in the terminal velocity as a result of the proximity of the tube wall to the falling sphere was less than 2% of the unbounded terminal velocity for the case of d/D = 0.2. In the present work the diameter of the tube is 200mm. The maximum particle diameter is less than 25.4 mm; so the maximum parti- cle to tube diameter ratio is about 0.125. This selection will provide negligible wall effect. The method of Thomas and Walter (1965) was used to estimate the distance (Le) required for the sphere to attain their terminal velocity using the following equa- tion: 17opgafiF '- =——— (4.5) e M 2 “0 where OP is the particle density pr+ pF 0‘6 ‘ —96'f'—F M' = 4/3 flpf a' 76 pf = fluid density a' = doublet radious = a3(1 + k3) k = ratio of large to small sphere diameter a = small particle radius F = (M-M')g = net gravational force M = 4/3 h a' pp and g = gravitational acceleration The maximum entrance length needed was thus estimated to be less than 200 mm. The designed tube height is 1200 mm, allowing for at least 800 mm of test section. Two sections of 300 mm each were marked along the tube length to improve precision. First, a section of 50 mm from the top was left empty to be able to draw a vacuum on the solution without drawing it out. This section is then followed by a 200 mm liquid filled section for attaining terminal velocity. Following this there are the two sections of 300 mm which are sepa- rated by 150 mm. A schematic diagram of the details of the drop cylinder is shown in Figure 4.2. Only measurement considered are those reproduced in the two equal sections. Finally, a section of 200 mm is allowed for the bottom edge effect. The accuracy of this arrangement is seen in the results presented for a Newtonian fluid. The spheres are supposed uibe free from any attachment and fall under gravity along their line of centers. Extreme care is required to have the spheres dr0pped axially in the fall tube at the start of each experiment. A special centering device, shown in Figure 4.3, was designed for each pair. This device, a Plexiglass funnel positioned in the center of the tube cover. The lower end 77 t. Vaccum section E5. 5 Terminal Velocity Section First Measurement Section 1200 mm Measurement Separation Distance Second Measurement Section Bottom Edge Effect Section .1. J N a S a a a: l_ —"'l 200nm L— Tube Diameter Figure 4.2. Schematic diagram of drop cylinder. 78 Figure 4.3. Photograph of centering device. 79 of the funnel is immersed in the fluid and is only 0.001 inch larger than the large sphere. This design yielded good degree of alignment of the pairs along the axis before they enter the fluid. Friction in the funnel was eliminated by making the passage very smooth. The tube wall thickness was 1/4" so as to withstand the pressure of the highly viscous fluid. 4.3 Materials 4.3.1 Test Fluids For the purpose of isolating the elastic effect on the trans- lation of the doublets, a Viscoelastic fluid that has a constant shear viscosity is needed. This fluid has the characteristics of what is referred to as a second order fluid with constant shear viscosity and constant normal stress coefficients. lim N 1 -> o (u 6)) =60 and 01= $5 = 20010 (4.6) where 01 is the first normal stress coefficient N1 is the first normal stress difference, 7 is the shear rate 10 is a time constant In the present work two media are used, a Newtonian corn syrup made by A. E. Staley and a Viscoelastic solution of Separan® in this corn syrup. The syrup was chosen due to its clarity and high vis- cosity at the experiment conditions which was 1580 poise at 25°C. 80 Terminal velocity was obtained over a short distance down the tube and low Reynolds number flow around the spheres was obtained because 4 to 10‘5). of such fluid viscosity. The range of Re is between (10' The Newtonian medium was used for the purpose of confirming the accuracy of the drag measurements. The Viscoelastic medium was pre- pared by dissolving small quantities (0.2 wt.%) of Separan AP30 synthetic polyacrylamide, manufactured by Dow Chemical, in the corn syrup. The polymer was sprinkled in the syrup at different depth in the preparation tank and left for two days to swell. Then a slight rotation of the solution to achieve uniformity of the polymer con- centration. The solution obtained was fairly clear and homogeneous. Separan AP-30 has good thermal stability below 210°C, and good resistance to shear degradation. 4.3.2 Preparation of Doublets Spheres choosen for the experiment were steel ball bearings having extremely small tolerance on diameter and sphericity- The sphere diameter was measured carefully at several points. The density was obtained from the measurements of weight and volume of sphere. Doublets were formed by joining two spheres together over a size ratio of 1 to 7. The joining process was performed very care~ fully. A commercially available adhesive (Super glue) was used, a small drop was enough to bond the spheres together at a minimum contact point. The surfaces of the spheres were wiped clean and no foreign material such as glue was on them as they were dropped in the cylinder. Figure 4.4 is a photograph of several size doublets. 81 .889 1‘8, “:9 :t “A" 4' Figure 4.4. Photograph of several size doublets. 82 4.4 Rheological Properties of Test Fluids An R-16 Weissenberg Rheogoniometer was used to measure the flow properties of all test fluids. Steady shear measurements were carried out to obtain the torques and forces on the cone and plate and angular velocity of the rotating plate. A cone angle of 0.5522° and a plate diameter of 7.5 cm were used in the measurements. To measure normal stress difference, piezoelectric load cell (922E) connected through an amplifier to a storage oscilloscope were used here. The cell was calibrated before the measurements by placing a certain weight and recording the voltage output. The first normal stress difference is then calculated by: 2HAnlg __f_ "1‘4 (mm c (6/4)d§ where An1 is the oscilloscope steady-state reading in volts 9 is the gravity force H is the transfer function of the cell in mass/volt (11 is the plate diameter and the factor (2) is for the force on the lower and upper platens. The first normal stress coefficient is then obtained by Equation (4.6). A standard oil supplied by ASTM was used for viscosity cali- bration. Temperature control was also used to obtain viscosity- temperature data for comparison with the manufacturer's data. This 83 comparison is shown in Figure 4.5 to be excellent. The Newtonian corn syrup was also characterized at various temperatures. The Viscoelastic Separan solution in corn syrup was tested at the experi- ment condition. The effect of inertia on normal stress measurements was studied very carefully. No normal force was observed with the pure corn syrup. Since the Separan solution has the same density almost as the corn syrup, inertial effects are absents in these measurements.. Corn syrup viscosity vs. shear rate is shown in Figure 4.6 and 4.7. The temperature dependence of the viscosity is given in Figure 4.8. Figure 4.9 shows the viscosity vs. shear rate of the elastic fluid. Constant viscosity is observed up to a 1 shear rate of 5 sec- . For the first normal stress coefficient 1 Figure 4.10 show that 01 is constant again up to y z 2 sec- . In any case our settling experiments involve shear rates less than or 1 as will be seen in Table 4.3. The first normal equal to_1 sec- stress difference is shown in Figure 4.11. The flow properties of all test fluids are summarized in Table 4.1., A new sample was used TABLE 4.1.--Viscosity, relaxation time, and density of test fluids Test Fluid Temp C° 00(Poise) 10(sec) of(gm/cm3) Standard oil 25 740 -- Corn syrup 200 25 1580 -- 1.4288 .2% Separan 25 1760.57 0.26 1.3501 84 IO 0 Experimental data - Line of manufacturer data I'll u ( PULSE) S Line of manufacturer data [<91 1 1 1 1 l 2.2 2.4 2.6 2.8 3.0 3.2 3.4 100V T (’k') Figure 4.5 Calibration of Weissenberg Rheogoniometer by ASTM Fluid. 85 [0 I 1 [III i 2, (5” Fl IIITII u( poise) [0 W l lllll I [[11111] I I IIJIIII IO [0 IO 7' (586") Figure 4.6. Corn syrup viscosity vs. shear rate at 25°C. 86 I03:- _\ : ——t + + 4.. 0.1 .83 " a 2 \_ IO :— 2 : ’0’ J ILIIJHII l [JILHII J 10" 10° 10 o -I 7(sec ) Figure 4.7. Viscosity vs. shear rate for corn syrup at 36°C. 87 IO ll'III 1 1 I0 11 (pOISC) IIIIIIII .3 IO r4 I llIllll Rheogoniometer data ' Manufacturer value 1 l 1 L J I 0 2.4 Figure 4.8. 26 2.8 3.0 3.2 3.4 100/ 7‘ (“k") Viscosity vs. temperature for corn syrup. . 3.6 3.8 88 10“ C ”44—9—1 .— a”‘s ‘3 «1’0 5' .3 : o L Q- . V b .11 102:- .- J l lllllll lillllflj l 70" IO 10' 7" ( sec") Figure 4.9.--Viscosity vs. shear rate of 0.2% separan in corn syrup at 25°C. 89 [III I I [VIII Imiill ,0 J l lJ-Jlll.‘ I illlllll -’ o . l0 70 lo - -I 7’ ( sec ) Figure 4.10. Normal stress coefficient vs. shear rate of 0.2% separan in corn syrup. 90 IIIII I0 I l ITII'I' NI(N/m‘3) .2 I0 I I IIIIII ’0 J l LIIJIII J l I 1 Hill I I 11111]! Io°- IO 10" 7" ( sec") 8 First normal stress difference vs. shear rate Figure 4.11. . for 0.2% separan 1n corn syrup. 91 at each shear rate. This was necessary to avoid any error due to material degradation or insufficient relaxation time after shearing the sample. Each time a sample is loaded twenty minutes were allowed for the material to relax before a measurement is taken. Evaporation of the sample was cut to a minimum by applying a thin film of silicon oil of comparable viscosity to the exposed sample in the gap between the cone and platen. 4. 5 Experimental Procedure The Newtonian fluid was used first for reproducing available theoretical results and confirming the suitability of the experimental arrangement. Pure corn syrup filled the drop cylinder up to a level that is enough to have the centering device immersed. The cylinder was left in a constant temperature room for few days to allow the entrapped air to escape and thermal equilibrium to be reached. Appli- cation of a vaccum on the cylinder helped further speed up the process of getting rid of the air bubbles. Each pair of spheres was released carefully in the centering device to avoid any eccentricity. The terminal velocity was recorded with two electronic timers, were each section is timed by a separate timer. A period of at least twenty minutes was allowed after each test to avoid any error due to disturbances. This was netessary due to the highly viscous fluids used. This time was also higher for larger doublets. The same procedure was followed for the visco- elastic medium. It is worth noting that it took a longer time to clear the solution from air bubbles. 92 The orientation of non-spherical particles is an important question to be addressed in sedimentation experiments. Photographs in Figures 4.12, 4.13 show that for both Newtonian and Viscoelastic medium the orientation in which the line of centers coincides with the axis of the cylinder and the larger sphere is underneath, is stable and maintained throughout the fall. The same behavior was predicted by Boun1(1977a) for transversely isotropic particle in a quiescent field, as discussed earlier in Chapter I. 4.6 Results 4.6.1 Previous Observations Non-Contactinngarticles For the purpose of completeness non-contacting spheres were also dropped with two different initial separation distances. This was done in both the Newtonian syrup and the elastic liquid. As expected, in the Newtonian syrup the initial separation has no effect on the spheres as they settle in the tube. On the other hand, in the elastic liquid we observed convergence of the spheres for small initial separation distance and divergence for large separation distance as they settle in the tube as can be seen in Figures 4.14 to 4.20. These phenomena were also reported by Riddle et al. (1977). The fluids used by Riddle were of considerable shear thinning. In our case the elastic fluid is non-shear thinning and thus we attribute these observations to the elasticity of the medium. 93 Figure 4.12. Photograph of a doublet in a Newtonian fluid. 94 Figure 4.13. Photograph of a doublets in a Viscoelastic fluid. 95 Figure 4.14. Photograph of particles separated by a distance in Newtonian fluid. Top section of the test column. 96 Figure 4.15. Photograph of same particles at the bottom 6f the test cylinder. 97 Figure 4.16. Photograph of two particles in Viscoelastic fluid at the top of the test cylinder. ' 98 Figure 4.17. Photograph of same two particles after some distance down the tube. 99 Figure 4.18. Photograph of same particles as they converge. 100 Figure 4.19. Photograph of particles separated by large critical distance at the top of the cylinder in a Viscoelastic fluid. 101 Figure 4.20. Photograph of same particles as they diverge down the tube. 102 4.6.2 Contacting Particles in Newtonian Liquid The theoretical and experimental values of the correction factor fN to STokes' law for doublets of various size ratios in the Newtonian corn syrup are summarized in Table 4.2. The TABLE 4.2.-~Newtonian fluid results. a = small sphere radius (cm). k = ratio of radii, fex = exp. drag correction factor, ftheo = theo-drag correction factor a I k fexp. ftheo .1587 cm 1 1.216 1.29 .1587 cm 2 2.08 2.09 .1587 cm 4 3.97 4.02 .1587 cm 5 5.03 5.01 .1587 cm 7 6.99 7.01 experimental values have been obtained from the following equation F G (4.8) fN = 6nuUa where Fe = (473)6(op-6f)(1+k3) a3g where g is the acceleration a is the small sphere radius k is the ratio of radii (k.: 1) pp and of are the particle and fluid densities 103 Theoretical values of fN were given by O'Neill (1969a). The maximum discrepancy between the experimental and theoretical values of fN is seen to be at k = 1; and amounts to 6%. This is believed to be due to experimental and wall effect error. Nevertheless, it should be mentioned that wall effect contribution to this value is very negligible. This conclusion is reached on the basis of the results of larger k values. So the error is mostly due to small errors in measuring the terminal velocities. It is to be noted that creeping flow was ensured since the values of Reynolds numbers are between lO'4 to 10's. The results obtained here confirmed the suitability of the cylinder design for carrying out the experiment for the elastic medium. 4.6.3 Contacting Spheres in Elastic Fluid The results of this section are to be presented in terms of a correction factor which accounts for the deviation in drag coefficient from Newtonian value due to the presence of fluid elasticity. The correction factor was obtained for several values to the product of 10 and the average shear rate which is in the case of doublets of spheres in creeping flow can be arbitrarily defined as XoYa v = IOU/a(1 + k) (4.9) where 20 is a relaxation time in seconds and all other quantities are as defined before. 104 The various values of the product Aoiav were obtained by varying the doublets size ratio. In general, for creeping flow, the drag coeffi- cient can be written as: cd = 2F (4.10) of 6202(1 + k2)n where F is the drag experienced by the doublets. The correction factor can then be written as: C Xe = 599- . (4.11) ds where Cde is the drag coefficient in the elastic fluid C is the drag coefficient for a Newtonian fluid of the ds same viscosity as that of the elastic medium It is worth noting that in calculating Cde’ we set F = FG (4.12) while in the calculation of Cds’ we set F = 6nuUa fN,theo (4.13) with the viscosity 0 and the observed terminal velocity U are taken from the experiment with the elastic fluid. Flow around a doublet of spherical shape in a non-viscometric and the shear rate at the 105 surface varies from point to point over the surface., This case is the same as for flow around a single sphere. An average shear rate defined by Yav = U/a is usually used there. In our case 7 is to be defined as av yav = U/a(1+k) (4.14) At this point we will define YavAo 5 Ne* ‘ (4.15) In Table 4.3 the results for the elastic fluid are presented. These results include the size ratio k, the small sphere radius a, Ne*, the average surface shear rate Iav’the correction factor due to elasticity Xefor the doublets and the correction factor for the single spheres XeS given by Boger and coworkers (1980). These results cover a range 4 to 10'6 as shown in the Table of raw of Reynolds numbers between 10' data in Appendix B. The results of Table 4.3 are also shown in Figure 4.21 and 4.22. The value of Ne* was varied by changing the size ratio of the spheres. Thus in Figure 4.21 Ne* values were obtained from doublets whose small sphere radius is 0.0794 cm and k was increased from 1 to 4. This arrangement covered a range of Ne* from 0.0129 to 0.064. In Figure 4.22 the small sphere radius is * 0.1587 cm and k is increased from 1 to 7 to cover a range of We from 106 1.0 o .9 '<_=4 XES, .8 '— L Xe 1 O Xes (single spheres, Boger, 1980) 4 Xe (contacting spheres) 7 T' I l 1, I I .01 .02 .03 .04 .05 .06 We* Figure 4.21. Correction factor for drag coefficient deviation from Newtonian Value due to presence of elasticity VS Ne*. 107 Xes Xe k 2 .8 "‘ o Xes (single spheres Boger 1980) AIXe (contacting spheres) 6 k=7 7 1 1 I 1 1 1 11 1 1 1 l 1 L 1 111 10'2 10'1 10'0 He* Figure 4.22. Correction factor for drag coefficient deviation from Newtonian values due to presence of elasticity vs Ne*. 108 Table 4.3.—-Elastic fluid results k a(cm) yav(sec'1) Iav Ao=we Xe Xes 1 0.0794 .05 0.0128 .99 1 i 1 0.1587 .1 0.0261 .98 1 ‘ 1.5 0.0794 .07 0.018 .91 1 2 0.0794 .104 0.027 .87 1 ‘ 2 0.1587 .22 0.056 .84 1 4 0.0794 .25 0.064 .82 1 4 0.1587 .51 0.133 .80 1 5 0.1587 .695 0.178 .76 .92 6 0.1587 .879 .228 .74 .87 7 0.1587 1.05 .272 .74 .85 0.02611x10.272. These figures show that there is a significant reduc- tion in drag coefficient of the elastic medium below that of Newtonian liquid for contacting spheres. It was further seen that this reduc- tion is dependent on Ne*. We observed that there is a linear reduc- tion over Ne*=n0-0128 to 0.035 where a reduction of 15% was seen. For Ne* between 0.035 and 0.064 a leveling off is observed and a reduction of 18% is reached. As Ne* is further increased-above 0.064 we observed further reduction in the drag below the Newtonian value. The maximum reduction was seen to be 26% at Ne* = 0.272. The case of single spheres doesn't show deviation except at Ne* = .13. A linear reduction is seen there which reaches a maximum of 15% at * We = 0.272. 109 At this point we would like to see how the elastic effect predicted by the numerical scheme of Chapter III compares to the results of this chapter. Itharticular, we will compare the results for the value of k = 1.5. In Chapter III we defined we = (0, + 20,) g- (4.16) where V1 = 01/00 V2 = Wz/UO Ipland¢m_ are the first and second normal stress coefficients In this work we could measure- only 01. It is a common practice to assume that 02 = - .101, so We ((1)1 ' 2X "11:71) 3%- : '8w1-a_%_ (4.17) O ‘ 0 According to Equation (4.6) W1 =‘2UOAO ' (4.18) Therefore Ne 1.6 00 U/a (4.19) So for k 1.5 and a = 0.0794 cm, U = 0.0145 cm/sec. (from Appendix B) we get Ne = 0.076. Thus using F02 and F12 for k = 1.5 in Table 3.6, we can evaluate a corresponding theoretical value to Xe as Xt using F = F02 + weFlz (4.20) 110 to get Xt = 0.944 while the experimental value Xe = 0.91. In Table 4.4 we present the results of k = 1.5 and 2. It is noticed that experimental values are in close agreement with those predicted by theory. TABLE 4.4.--Comparison between theoretical drag coefficient ratio Xt and experimental drag coefficient ratio Xe: k He He* 1-Xt 1-Xe 1.5 0.076 0.018 .06 1 .03a .09 2 .13 0.027 .1 1 .02a .13 aThis uncertainty is due to error in numerical scheme, already discussed. CHAPTER V CONCLUSION AND RECOMMENDATION 5.1 Conclusion This investigation has focused on determining the elastic effect of a non-shear-thinning medium on the sedimentation of con- tacting particles. The flow conditions are those of creeping flow with Reynolds number in the range of 10'4 to 10'6. As far as we know, we are the first to observe the elastic effect on the sedimen- tation of contacting particles with a non-shear-thinning (Boger) fluid. The suitability of our apparatus was checked with a viscous, Newtonian corn syrup in which our observations on orientation as well as sedimentation speed were consistent with previous reports in the literature. The experimentally observed settling speed on these particle pairs is reproducible to within 2 percent. In the elastic fluid medium, this study provided new infor- mation about a stable orientation for unequal spheres. This orien- tation is that of the line of centers along gravity with the larger sphere leading. It was further observed that even if the initial orientation was not coinciding with the terminal orientation, the particles will reorient themselves to the terminal and stable orientation within a length of at most 15 cm. For completeness, we repeated the experiments which were carried out by Riddle et al. 111 112 (1977) on spheres separated by several distances. Convergence was seen for particles with small initial separation and divergence was seen with large initial separation. It is to be noted that while the fluid used by Riddle was considerably shear thinning, we observed similar behavior with the Boger fluid, so the above behavior might be attributed to elasticity of the medium alone. Measured settling velocities for unequal spheres showed a significant reduction in the drag coefficient below the Newtonian values. The reduction was seen to be a function of Ne* = Iavxo’ where a deviation of 26 percent is observed over a range of Ne* between 0.0128 and 0272, based on a dimension (1 + k)a of the two spheres. A theoretical analysis was carried out to investigate the translation of these contacting particles along their line of centers, in an effort to obtain the elastic effect of a second order fluid medium of a constant viscosity on the drag force on the con- tacting spheres surface. A volume integral was developed for the drag contribution at 0(We) involving the zeroth order, Newtonian velocity field. Steps were taken to insure accuracy at each stage of the numerical solution by checking condition of matrix in solving a set of linear equations repeatedly; comparing asymptotic limits of all six Hankel transforms at both ends (n + O and n +1w) with numeri- cal estimates, taking care to obtain simple analytical approxima- tions to the behavior of the integrand at large 5 as well as small 5. Since the function f(s,€) is evaluated thousands of times in the 113 numerical scheme, accurate evaluation of f was the heart of the numerical steps in this work. Thus we carried out a detail error analysis and obtained an upper bound on the total error in evaluating f(s,F,). Our analysis revealed that the error increased with decreasing a over 0 < a < 1 and in the increasing 5 values above 5 = 6. Using the worst case error an upper bound on the error in the total result including 27,000 function evaluations, would be 10'4 for a = 0.2 and 10'7 for a = 0.5. In addition, there is the 2 - error in fitting exponentials to the large 5 behavior of f,%§, ggi. This fit was necessary to evaluate the complete Hankel transform of these functions. The error of the fit was worst for azflag2 and this error too increased with decreasing a over 0 < a < 1. We have tabulated results in Table 3.7 above a = 0.33 for which these errors add up to only 20 percent of the final result for reduction in drag. These results indicate a definite nonzero elastic contribu- tion to drag for the a values listed. The contribution of the integrand involving n and 5 alone was investigated over the entire range of g and n. This investi- gation revealed that the region near the stagnation point on the larger sphere is the most important region for estimating the elastic effect. The region near the contact point of the two spheres has no contribution at all. These results are in accordance with findings of Leal (1975) in his analysis of the drag on slender cone in a second order fluid, Leal found that the major contribution to a similar volume integral was obtained from the fluid close to the particles; however, because of the slender body approximation, he 114 was able to omit the stagnation points at both ends of the particle for very long rods. The numerical results obtained for the integral of Equation (3.46) for a 3 0.33 show an appreciable elastic effect on the drag experienced by the unequal contacting spheres. This effect was larger for larger size ratio lye Comparing these theoretical predic- tions with the experimental values for k = 1.5 and 2 showed that experiments are in good agreement with theory. 5.2 Recommendation The theoretical and experimental analysis carried out here has established the existence of an elastic effect due to the medium on the drag force of two contacting spheres of unequal size translating along their center line in creeping flow. The major problem and difficulty in the numerical solution is in evaluating the function f(s,g) at large values of s and small values of a. The fitting procedure also involved an additional error with decreasing a over 0 < a < 1. Therefore, in order to obtain results for low a values some other methods for evaluating the function f(s,§) and fitting 2 f(s,g),fl%§flandfi%§§l should be pursued in future work. APPENDICES 115 APPENDIX A PROOF OF VANISHING INTEGRAND 0F EQUATION (3.20) 116 APPENDIX A PROOF OF VANISHING INTEGRAND OF EQUATION (3.20) Proving that (21’ vgo- fi : 131) = 0 (A.1) Vf Tensorial notations are going to be adapted for this proof. We consider the integrand term by term as follows HI: Vvo = - P12 + Zugl: Vyo = [-P1 (Sijeiej + “(3:2 eiej + 13:13— 8131.)] 23:" - (-PlalJ :izm) (euej° emen) + “ (2:;1 :::m )(eieJ emen) 117 118 3v . 3v 1; om , + u(3x1. 5xn )(ejei' amen) - p 5 . 3:99. 5 5_ + u(.3!li.i!9fl. 5 5_ 1 ij 3x in jm 3x. 3x in am n j n + avij avom 6. 5. ) 3x1 3xn jn 1m Since m = j and n = i av . 3v . 3; . 8v . 3; . =-p 5,,_2.1+u 1103 _11 _9.1_) 1 13 3x1 axj 8x1 3x1 3x1 and i f j + 1 - v30) (A-2) o A+ A .- 11 (V!1°VVO + W by the same procedure it can be proven that would yield the same quantity as in (A-2) thus the integrand of (A-l) vanishes. APPENDIX 8 EXPERIMENTAL RAN DATA 119 TABLE B.1.--Experimental Raw Data a(cm) t(sec) U(cm/sec) Egynes) (gynes) Re 0 0794 1 3854.8 0.0079 26.44 26.72 1.9 x 10'6 0.1587 1 949.3 0.0316 208.96 213.74 7.7 x 10‘5 0.316 1 247.8 0.121 1692.66 1659.59 5.9 x 10"5 0.0794 1. 2067.7 0.0145 57.83 63.9 3.5 x 10"6 0.0794 2 _ 1214.7 0.0247 118.96 136.43 6.1 x 10"6 0.1587 2 293.6 0.1022 952.1 1133.4 2.5 x 10‘5 0.0794 4 306.7 0.0978 857.2 1046.5 2.4 x 10‘5 0.1587 4 74.3 0.403 6871.79 8579.62 9.8 x 10‘5 0.1587 5 45.3 0.662 13329.11 17538.3 1.6 x 10'4 0.1587 6 30.7 0.976 22955.69 31021 212 2.4 x 10'4 0.1587 7 22.5 1.33 36390.6 49309.7 3.2 x 10‘4 120 BIBLIOGARPHY 121 BIBLIOGRAPHY Acharya, A., R. M. Mashelkar, and Ulbrcht, J. Rheol. Acta, 15, 454 (1976). Adler, P. M. Rheol. Acta, 11, 288-295 (1978). Bacon, L. R. J. Franklin Inst., 221, 251 (1936). Brenner, H. Chem. Eng. Sci., 15, 242-251 (1961). Brenner, H., and Happel, J. Chem. Eng. Sci., 12, 599 (1964). Brenner, H. Low Reynolds Number Hydrodynamics, Noordhoff (1965). BroadbenE, J.)M., and Mena, B. J. Non-Newtonian Fluid Mech., 2, 1 1974 . Brunn, P. Rheol. Acta, 15, 104-119 (1976). Brunn, P. Rheol. Acta, 15, 163 (1976). Brunn, P. J. Fluid Mech., 82, 529-547 (1977a). Brunn, P. Rheol. Acta, 15, 461-475 (1977b). Brunn, P. Rheol. Acta, 18., 229-243 (1979). Brunn, P. J. Non-Newtonian Fluid Mech., 1, 271-288 (1980). Caswell, B., and N. Shwarz. J. Fluid Mech., 15, 417-426 (1962). Caswell, B. Chem. Eng. Sci., 25, 1167 (1970). Caswell, B. Appl. Mech. Div., 22_ASME (1977). Cooley, M. D. A., and M. E. O'Neill. J. Camb. Phil. Soc., 55, 407-415 (1969a). Cooley, M. D. A., and M. E. O'Neill. Mathematika, 16, 37-49 (1969b). Cox, R. G., and H. Brenner. Chem. Eng. Sci., 22, 1753-1777 (1967). 122 123 Chhabara, R. P., P. H. T. Uhlherr, and Boger, D. V. J. Non- Newtonian Fluid Mech., 5, 187-199 (1980). Chan, P. C. H., and L. G. Leal. J. Fluid Mech., 21, Part 1, 131-170 (1979). ‘ Dean, N. R., and M. E. O'Neill. Mathematika, 19, 13-25 (1963). Faxen, H. Ark. Mathematik Astron. Fys., 11, 1 (1923). Frankel, N. A., and A. A. Crivos. Chem. Eng. Sci., 21, 847-953. (1967). Gauthier, F., H. L. Goldsmith, and Mason, S. G. Rheol. Acta, 19, 344-364 (1971a). Gauthier, F., H. L. Goldsmith, and Mason, S. G. Trans. Soc. Rheol. 15:2, 197-330 (1971b). Gieskus, V. H. Rheol._Acta, 3, 59-71 (1962). Goldman, A. J., R. G. Cox, and Brenner, H. Chem. Eng. Sci., 11, 1151-1170 (1966). Goldman, A. J., R. G. Cox, and Brenner, H. Chem. Eng. Sci. 21, 637-659 (1967). Goren, S. L. J. Colloid Interf. Sci., 15, 94-86 (1971). Goult, R. J., and Hoskins. Computational Methods in Linear Algebra. New York: John Wiley (1974). ’ Highgate, D. J., and Whorlow, R. w. Rheol. Acta, 2, 569-576 (1970). Hoffman, R. L. J. Colloid Interf. Sci., 45, 481-506 (1974). Kataoka, T., T. Kitano, M. Sasahara, and Nishijima, K. Rheol. Acta, _11, 149-155 (1978). Leal, L. G. J. Fluid Mech., 59, 305-337 (1975). Leslie, F. M. and R. Z. Tanner. 11, Pt. 1, 36-48 (1961). Lin, C. J., K. L. Lee, and Sather, N. F. J. Fluid Mech, 41, 35-47 (1970). Lorentz, H. A., Abhand. L. The Gret. Phys. 1, 23 (1906). Michele,(J., R. Patzold, and Donis, R. Rheol, Acta, 15, 317-321 1977). 124 Nir, A., and A. Crivos. A. J. Fluid Mech., 52, 209-223 (1973). Oldroyed, J. G. Proc. Roy. Soc. A, 215, 278 (1958). O'Neill, M. E. Mathematika, 11, 67-74 (1964). Papenhuijzen, J. M. P. Rheol. Acta, 11, 73-88 (1972). Riddel, M. J., C. Narvaez, and Bird, B. J. Non-Newtonian Fluid Mech., 2, 23-35 (1977).' Russell, W. B. J. Rheol., 24, 287-317 (1980). Sakai, K., Adachi, K., and N. Yoshioka. J. Non-Newtonian fluid Mech., 1, 107-125 (1977/1978). Sigli, D., and Countanceau, J. J. Non-Newtonian Fluid Mech., g, 1-21 (1977). Soni, K., and Soni, R. P. SIAM J. Math. Anal., 4, 466-481 (1973). Stimson, M., and Jeffery, G. B. Proc. R. Soc., A111, 116-119 (1926). Sun, K., Jayaraman, K. J. Chem. Eng. Commun., 29, 137-155, (1982). Sutterby, J. L. Soc. Rheol, 1114, 559-573 (1973). Thomas, R. H., and Halters, K. Rheol. Acta, 5, 23-27 (1966). Zeichner, G. R., and Schowalter, H. R. AICHE J., 11, 243-254 (1977).