..‘. . . .‘ THESIS LIBRARY Michigan State University 1., This is to certify that the thesis entitled A NUMERICAL STUDY OF THE STOKES FLOW FIELD AROUND CONTACTING SPHERES TRANSLATING AXISYMMETRICALLY presented by Kevin L. Nichols has been accepted towards fulfillment of the requirements for M. S. Jegee in Chemical Engineering Date F'Oé' 024,1/7’5 0-7639 MS U is an Affirmative Action/Equal Opportunity Institution MSU LIBRARIES “ \— RETURNING MATERIALS: Place in book drop to remove this checkout from your record. FINES will be charged if book is returned after the date stamped below. A NUMERICAL STUDY OF THE STOKES FLOH FIELD AROUND CONTACTING SPHERES TRANSLATING AXISYMMETRICALLY By Kevin Louis Nichols ‘ A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Chemical Engineering 1985 ABSTRACT A NUMERICAL STUDY OF THE STOKES FLOW FIELD AROUND CONTACTING SPHERES TRANSLATING AXISYMMETRICALLY By Kevin Louis Nichols A scheme was developed for the numerical evaluation of the stream function, velocities, rates of deformation, and vorticity around two contacting spheres of equal or unequal radii translating, axisymmetrically, in creeping flow, through a Newtonian fluid. Special procedures were developed to evaluate the various Hankel integrals, occurring in the calculations, accurately enough to obtain estimates of quantities such as shear rate with an error of no more than 2 per- cent. The detailed Newtonian flow field calculations were then used to estimate a first order effect of elasticity in terms of small Heissenberg number (relaxation time times average surface shear rate) on the drag around the particles. This effect turned out to be insig- nificant for arbitrary radius ratios. In contrast, the results of settling experiments done by Jefri (1983) in a viscoelastic fluid showed that while there is no reduction in drag for equal radii pairs there is a significant (10% or higher) reduction in drag for unequal radii pairs. Another mechanism for the observed drag reduction is suggested, namely, changing the no slip boundary condition to a slip boundary condition. To my wife, Joanne ii ACKNOWLEDGEMENTS I would like to thank Dr. K. Jayaraman, my advisor, for his instruction throughout the development of this thesis. I would also like to thank Michigan State University for their financial support during this period of my study. TABLE OF CONTENTS LIST OF TABLES . LIST OF FIGURES Chapter I. INTRODUCTION II. BACKGROUND . III. THE FLOW FIELD--(ANALYTICAL) . IV. CALCULATION SCHEME V. ACCURACY OF CALCULATIONS Scheme Used to Check Accuracy of Calculations VI. RESULTS . Velocity . . . . . . . . . . . Rate of Deformation Tensor and Vorticity Surface Shear Rate . . . . . . Attempt to MOdel Drag VII. CONCLUSIONS AND RECOMMENDATIONS . BIBLIOGRAPHY iv Page Vi 22 23 37 39 48 54 68 71 82 84 Table 5.1 5.2 5.3 6.1 6.2 6.3 6.4 LIST OF TABLES Change in integrals Ial’ a = 1, . . ., 6 using 50 quadrature points as compared to 40 quadrature points . . . . . . . . . . . . . I1 throuth 16, a = .2, n = 15 Location of w = O streamlines for a = 1 case along r axis (a = O) . . . . . . Criteria for identifying stagnant, stretch dominated, rotation dominated, and shear dominated regions of fluid . . . . . . . Average surface shear rates Quadrature points for a = .2 case F0 and F for a = 1, .5, and .2 1 Page 30 33 36 66 69 79 81 Figure 1. 01010010101 #0)“) mmmmmmm 1 .10 .11 .12 .13 .14 LIST OF FIGURES First normal stress difference and Xe vs. shear rate for 0.2% separan in corn syrup Tangent sphere coordinate system relative to cylindrical coordinates . . w and pi vs. r at z = 1 for a - 1 case 1 for a .2 case . w and Yi vs. r at z Streamlines, w, for a 1 case . Streamlines, w, for a .2 case Streamlines, w, for single sphere case Streamlines, w, (x104) in the contact region for = 1 case . . . . . . . . . . . Streamlines, w, (x104) in contact region for a - 2 case . . . . . Streamlines, w, (x104 ) in polar regions for 1 case . . . . . . . Streamlines, w, (x104) in polar regions for = 2 case . . . . . . . . Contours of vr (x102) for a = 1 case . Contours of v} (x102) for a = .2 case Contours of V2 for a = 1 case = .2 case . Contours of V2 for a v2 vs. r at z = 1 for a = .2 case . Contours at Drr (x103) for a = 1 case Contours of Drr (x103) for a==.2 case vi Page 10 33 34 4O 41 42 44 45 46 47 49 50 51 52 53 55 56 Figure mmmmmmm .15 .16 .17 .18 .19 .20 .21 .22 .23 .24 .25 .26 .27 Contours of DZZ (x103) for a = 1 case Contours of DZZ (x103) for a = .2 case Contours of Drz (x102) for a = 1 case Contours of Drz (x102) for a = .2 case Contours of vorticity (x102) for a = 1 case Contours of vorticity (x102) for a = .2 case Contours of vorticity for singls sphere case Stagnant (I), stretch dominated (II), rotation dominated (III), and shear dominated (IV) regions of fluid for a = 1 and a = .2 cases . . Dimensionless tangential surface shear rate for a = 1 case . . . . . . . . . Dimensionless tangential surface shear rate for leading sphere in a = .2 case . . . Dimensionless tangential surface shear rate for trailing sphere in a = .2 case . Dimensionless tangential surface shear rate for single sphere case . . . . . Values and locations of maximum dimensionless surface shear rates (ymax/va/R) for single sphere case, a = 1 case, and a = .2 case . vii Page 57 58 59 6O 62 63 64 67 72 73 74 76 CHAPTER I INTRODUCTION The objective of this work is to numerically evaluate the detailed flow field around two contacting spheres translating, in creeping motion, through a Newtonial fluid. As will be shown here, this may help in understanding the role elasticity plays on the drag around aggregates in viscoelastic fluids. Chhabra et al. (1980) conducted experiments with spherical particles settling in model polymer solutions to examine the influence of fluid elasticity on the drag coefficient for creeping flow around a sphere. Chhabra et al. point out that all experimenters before them did experiments in shear thinning fluids so no distinction could be made between the effects of shear thinning and the effects of elasticity. Chhabra et al. used second order fluids that, at average surface shear rates, i. involved in the experiment, were not shear thinning, and which displayed primary normal stresses that were proportional to i2. They carried out measurements at sufficiently small values of sphere Reynolds number so that inertia could be neglected. The idea then was that any difference between Newtonian drag and experimental drag would be due exclusively to elasticity. Chhambra et al. defined an average surface shear rate in their experiments as: Y=voo/R 1.1 where v is terminal velocity of sphere G) R is sphere radius The results of Chhabra et al.'s experiments show a reduction of drag occurred when the average surface shear rate was high enough that the shear modulus became shear rate dependent, and no change in drag occurred when the average surface shear rate was in the second order region for the model fluids used. Chhabra. et al. concluded from this that the drag on a sphere in a viscoelastic fluid was dependent on the elasticity of the fluid only when the fluid shear modulus became shear rate dependent. Jefri (1983) conducted settling experiments similar to Chhambra et al.'s experiments with two contacting spheres instead of a single sphere. The results of these experiments are shown in Figure 1.1. Xe on this figure is defined as Xe = CD (experimental) / CD (Newtonian) 1.2 where CD is the drag coefficient The average surface shear rate for two contacting spheres in the settling experiments is defined as: i=vw/(1+k)R 1.3- where the leading sphere--as the spheres translate axisymmetrically-- has a radius of R and the trailing sphere has a: radius of kR, legend A 5 data point from Jefri's experiments 104 First ‘ 1.0 normal _ stress .9 Xe difference (N/mZ) P .8 for two 103 -- contacting E - .7 spheres I - 6 102 -_— Second order 10 region I J Ill'lllI l 1 1111111 I I JIJJ‘I -1 2 1° 1 10 10 Figure 1.1.-~First normal stress difference and Xe vs. shear rate for 0.2% separan in corn syrup. 0 < k < w. As can be seen in Figure 1, in contrast to Chhabra et al's results, the experimental drag differs from the Newtonian drag well within the second order region for the model fluid (where the shear modulus is not shear rate dependent). Using Chhabra et al.'s logic, the observed change in drag is due exclusively to the elasticity of the fluid. Therefore, one can see if the drag on the two contacting spheres in Jefri's experiments can be properly modeled, then one may be able to predict and determine the effect of elasticity on drag in viscoelastic fluids. Determining how the elasticity affects drag in Jefri's experiments may also help in understanding other phenomena, for instance, the relation between fluid shear modulus and drag reduc- tion noticed by Chhabra, and the drag reduction found in very dilute polymer solutions. Jefri (1983) suggests a model for the drag in his experiments that involves an accurate calculation of the Newtonian rate of deforma- tion tensor over the entire flow field for two contacting spheres translating, in creeping motion along their line of centers, through a Newtonian fluid. To verify the accuracy of these tensor calcula- tions, it is necessary to evaluate velocity accurately over the entire flow field, which, in turn, requires an accurate calculation of the stream function over the entire flow field. These calculations are not trivial and the bulk of the work done here isaidevelopment of accurate numerical results from flow field calculations. Once the flow field details are established, the validity of Jefri's drag model is tested. CHAPTER II BACKGROUND The problem of creeping flow past two spheres translating along their line of centers, with equal velocity in a viscous fluid, was first solved by Stimson and Jeffrey (1926). They give an exact solution for Stokes flow past the two spheres not in contact. Cooley and O'Neill (1969) give exact analytical expressions for the stream function and forces acting on the spheres when they are in contact. Both Stimson and Jeffrey's and Cooley and O'Neill's results were derived for Newtonian fluids. Cooley and O'Neill's results are used extensively in this work. Chhabra' et al. (1980) point out in their work that the dimensionless group that arises as a result of fluid elasticity, in the absence of shear thinning, is the Weissenberg Number. For creep- ing motion around two contacting spheres this number may be defined by: we = wi / (1+ k)R 2.1 where A is relaxation time. Chhabra et al. also mention that any deviation from Newtonian drag in settling experiments should be a unique function of the Neissenberg Number. Therefore, as a first attempt at modeling the results in his experiments, Jefri used a regular perturbation expansion in powers of small weissenberg Number around Cooley and O'Neill's (1969) Newtonian solution for the drag. In the first approximation, only terms of 0(we) are retained, resulting in 1 - F = j - F0 + Ne(i, F1) 2.2 where F is the drag F0 is the Newtonian drag as calculated by Cooley and O'Neill, j is the unit vector in the direction of motion, F1 is the non-Newtonian (elastic) contribution of drag i - F1 was obtained by application of the Reciprocal Theorem as i - F1 = 2uvw jf g . 9:9 dV 2.3 Vf where u is fluid viscosity Voo is velocity of fluid at infinity Vf is the volume of fluid surrounding the particle D is rate of deformation tensor The rate of deformation tensor is given by = 1/2 (vV + VVT) 2.4 NO It is desired to see how well (2.2) will model the results of Jefri's experiments. To find this out, it is necessary to be able to calculate F1 accurately since F0 and We are easily calculated. This translates into being able to evaluate the integral in equation (2.3) accurately which is no trivial task and is what the rest of this work will deal with. The integral in (2.3) is a volume integral where the volume is defined as being the volume of fluid surrounding the particle. One problem in evaluating this integral is deciding what the boundar- ies of this volume are and picking points for numerical integration inside this volume. Another problem (probably the biggest one) is an accurate numerical evaluation of D-ng, more specifically D, at various points in the fluid volume, once the volume is defined. The calcula— tion of 0 involves a stream function expressed as a Hankel transform; also involved are derivatives of the stream function. Each of the above calculations involve some error which can overwhelm the correct answers if not properly handled. It appears that an evaluation of these errors and correct evaluation methods for the flow field expressions involved in calculating D have not been discussed by any other researchers. It is the objective of this work to do so here. Being able to calculate the flow field accurately could have an impact on other fields of study such as analysis of particle deposition on bisphere collectors, etc. Cooley and O'Neill (1969) furnish the stream function for the contacting sphere problem, but they never use it directly in any of their calculations. Instead, by making some simplifying assumptions and by using clever algebra, they derive an integral for the Newtonian fluid drag which is easy to evaluate numerically. This is useful for evaluating F0 in (2.2), but their work does not suggest ways to evaluate the stream function and its derivatives accurately. To cal- culate the rate of deformation tensor, 0, one needs to know over what ranges of space the stream function and its derivatives can be accurately, numerically determined. Davis et al. (1976) did later work on this problem that is of interest here. This was a study of separation of flow around two contacting spheres of equal radii. Separation of flow is identified by a stream surface, where the stream function equals zero, dividing the flow into two distinct flow regions. In their study, Davis et al . show surfaces (at which the stream function equals zero) near the contact area of the two spheres, which divide the flow field into an infinite number of distinct flow regions. In this case, there is an infinite number of nested ring vortices in the contact region separated from the rest of the flow field. The results of Davis et al.‘s work can be used to determine the accuracy of the calculations in the work done here. Also, flow separation for two contacting spheres (of unequal radii) is determined here and compared to the equal radii contacting spheres results, a study that does not appear to have been done prior to this. Reed and Morrison (1974) investigated the flow field around and inside of two touching fluid spheres translating along their lines of center. Their work is different from previous work in that they use biSpherical tangent coordinates directly in evaluating all. relevant flow field quantities such as velocity and strain rates. The calculations are done here in a similar manner. CHAPTER III THE FLOW FIELD-~(ANALYTICAL) It is desired to numerically calculate the various character- istics of the flow field around two rigid spheres in contact, trans- lating, in creeping motion along their line of centers, in a Newtonian fluid. Detailed expressions follow for velocity profile, stress dis- tribution, and vorticity which are developed from similar derivations done by Cooley and O'Neill (1969), Reed and Morrison (1974), and Happel and Brenner (1973). The geometry of this problem suggests use of the tangent sphere coordinate system. The tangent sphere coordinate system (n, g, o) is a rotational orthogonal curvilinear coordinate system related to cylindrical coordinates (z, r, m) by z=_2_§___;r=_ZD.—;¢:¢ 52+le €2+OZ or n = 2r ; E = “-EEL‘—‘S P = P 3'1 r.2_.,22 r.2+22 Figure 3.1 shows how the tangent sphere and cylindrical coordinate systems relate to one another. For this problem, the sphere lying in the +2 direction is assumed to have a radius of one and the sphere lying in the -z direction a radius (M: k, where 0 < k < m. In (n,g) 9 10 mmpmcwucoou pmuwcncwpzu on m>w»m—mc Emumxm mumcwucoou mcmsgm pcmmzmh--.fi.m mcnmwd N om as NH w a o a- m- NH- on- ON- _ F _ _ l b I . - _ — . Hoe. . -w H ,- .., c was. m." . a mo." m . u m H.- mo. 1. NH N 0H. 0.- mma. m mo.- .. on mmo.- H. mo.- I om Hso.- .u mmwo. em 0 o. «Sac. moo.-. 6 mo.- co. m o. c m o.- mm 11 coordinates the +2 sphere corresponds to g = 1 and the -z sphere to g = —1/k. This point where the spheres contact, (r = O, z = 0), corresponds to n = w and the region occupied by the fluid is given by -a < 5 < 1, O < n < m, where a = I/k. Assumptions made about the fluid are it is unbounded, incom- pressible, has constant density and viscosity, and is moving in the -z direction with a speed of V” far from the spheres. The Reynold number is assumed to be approximately zero which allows one to neglect inertial effects. The equations which govern the flow are 92: 32W l M 32" .3: dr uVm [:3r2 + r 3r +'537£ r2 3’2 gp_= 32v _1 av: azvz dz “Von [$521 + r 8r + 322 3'3 and in: in _Vz= ar +r I 32 0 3'4 where vr(r,2), vz(r,z), and p(r,z). The boundary conditions are Vr = 0, v2 = -1 at z = w from the way the fluid has been defined, and Vr = 0, v2 = 0 on the surface of the Spheres due to no slip condition. v¢ is assumed to be zero due to the symmetry in the geometry for this problem. The solution of (3.2), (3.3), and (3.4) is obtained in terms of a Stokes stream function, u, defined by 12 _ 1_84 . z- 1 34 v-rfi.v 7-7 3.5 where u : .32 _l_§_ 32 - - W-[Wirafisz—zJ “”0 3-5 The solution to (3.6) is straight forward when tangent sphere coordinates are used. To convert from cylindrical coordinates to tangent sphere coordinates, the metric or scaling coefficients = =M. .3133: h5 h5 2 . h¢ Zn 3.7 are used. Using these coefficients, velocity in terms of tangent sphere coordinates is given by =,- 2 + 2 2 84 . = 152 + n2)2 ENL vn EFL)— 733’ v5 4n 8n 3'8 Applying the chain rule for partial differentiation to the stream function allows one to transform velocity components rapidly from one system of coordinates to another. The velocity components (Vr’ v2) can be expressed in terms of tangent sphere velocity components (vn, Vg) as 13 Likewise, the Stokes stream function operator can be expressed in tangent sphere coordinates, so (3.6) becomes A“: n€2+n2 152+92j+§_52+92. 8—‘1w=0 ' 2 3n 2n an 35 2n 3€__ 3.10 The boundary conditions for (3.10) are w-fi-‘Pantahml‘ 3.11 4=Tgr%‘%ryzatn=o.g=o. 3.12 A solution for (3.10) satisfying the boundary conditions (3.11) and (3.12) is found from the solution to A2 W = 0 3.13 (3.13) is R-separable in tangent sphere coordinates (see Reed and Morrison, 1974, pp. 577-578). Let - G(€.n) m - (£2 + n31; 3.14 where G(€.n) = N(£) M(n) substituting (3.14) into (3.12) yields the following ordinary differential equations in (3.15) and(3.16) d2N 3.52—— - SZN = O 3.15 14 d2M dM 2 _ 2 2 = The solution of the linear differential equation with constant coefficients (3.15) is N(g) = Clsinhsg + Czcoshsg 3.17 while the solution to (3.16), a form of Bessel's equation, is M(n) = nI'C3J1(Sn) + mush] 3.18 Since 0 is finite as n + 0 from boundary condition (3.12), Y1(sn) must be omitted from the solution because it becomes unbounded logarith- mically as n + 0. Therefore, the solution to (3.13) is w = TO7DT_§773 Jf [c1(s)sinhsg + Cz(s)coshsg] J1(sn)ds 3.19 O Stimson and Jeffrey (1926) have shown that w + 20 is a sufficiently general solution to (3.10) where w is determined from (3.13). Using this fact, the boundary condition (3.12), and some algebraic manipu- lations, yields the following solution to (3.10) 00 2n2 4 = WW2] F(s.€)J1(sn)ds + (£2 + ”277 3.20) 0 where F(s,g) = (A + Cg)sinhsg + (B + Dg) coshsg 3.21 15 A, 8, C, and D are functions of s which must be determined from the boundary conditions. The following relations, obtained using tabulated Hankel transforms, are useful in writing the boundary conditions. 13—27274, =f e-Sm (IEI +%) 41(Sn)ds 3.22 'O (EYE—fifi/Z = [ 59-Slgldl(sn)ds 3.23 ' 0 Using the four boundary conditions in (3.11) and the above Hankel transform relations yields the following four equations Asinhs + Bcoshs + Csinhs + Dcoshs = -2e'5(1 + %) 3.24 Asinhds - Bcoshas dCsinhas + aDcoshas = 2e'sa(a + %) 3.25 Acoshs + Bsinhs + C [51.2115 + coshs] + D [320th + sinhs:] = 2e'S 3.26 Acoshds - Bsinhds C [élgflgé + acoshas] + D E:_0%h_q_s_ + asinhas] = 2e'o‘s 3.27. 16 Solving (3.24) through (3.27) for A, B, C, and D, and then substi- tuting into (3.21) yields the solution of the stream function, u, for the problem. The velocity expressions,(3.8) and (3.9) and the stress tensor expressions to be derived later involve first and second derivatives of the stream function with respect to n and/or 2. This implies evaluating the derivatives of the integral inside w. Let 1, =0] f(x/n,g) J1(x) 9715 3.28 Then the derivatives of the integral are given by 00 _ BI _ 8J1(x) dx 12 - a—n-L — f f(x/n,g) x 3x a, 3.29 o 2 2 13 = —2*§ I =1 f(x/n.€) X‘°'—-‘L-Lr2—a J a x 9’; 3.30 n x n 0 =.§l1 3fLX/n.€) .95 It} 35 of ag J1(X) n 3.31 - J2 - EM) 241951 dx 15 - 3g - J 85 X 3X '57 3.32 - 321 .. °° 32f(x/n.;1 d_x - I6 - E71 - f 3;?— J1(x) Tl 3'33 17 where the following equalities have been used because they ease com- puter implementation of the above integrals x = sn; f = (A + gC)sinhx§_+ (B + gD)coshx§. 3.34 n n M) = _ , 2 321] X x 3x x Jo(X) 01(X). x $154 (2 - x2)J1(x) - x Jo(x) 3.35 Now that the derivatives with respect to n and g of the integral inside the stream function have been established, the various derivatives of the stream function are defined as follows .22- - n _4 1 , , 3n2 3O (E? + n2)3/2 12 '1' [(52 + ”7).”; (EZrT—Isz-T] 11 4n 3 8 n +-(-—z——2-)-€ +11 2' (£2 + n2): 3°36 .34 = n ,, I _ 3n€ at (52 + nYV’" ‘* 142‘ + n‘)°/2 8n2€L I1 (n‘+€‘)’ 3.37 3211) = n l and (at + n[) 5] ' 1 E - 3112 .- 2 I3 + Zfigz + nrydli (£2 + nryal—z 12 _, _ 9n 15Tl3 *[W*Wlh 4 - ~4On2 48ml+ (£2 + n?)2 (E2 + nYIS + (Ez'tn2)“ 3'38 18 82 ____ n - _ 6571 E1 (5‘ + n‘lifz 16 “(Wm-‘1 + 15mg2 _ 3n I (5‘ +n‘)’/‘ WWI—Tl??? 1 4822 82 + "(7—an min-4.21.213 3.39 —qJ—~ Tl. 15+ - 1 - 3112 I Gian (Ezi'niI (5‘ + 411377* ( 2 “ £2 + n2) 5f : 15n22 - _ 3 16 +[(t‘+n7T’/r F. +11 1“” T873553 48n3€ (527+ nilr 3-40 Now the velocity expressions (3.8) and (3.9) can be calcu- lated directly using (3.36) and (3.4). Using the definitions from Happel and Brenner, the gradients of the velocity field, Vv. can be evaluated. The components of this dyadic are 2 2 2 ann = (€28: nilz [E4n:%%+-&§—€FILJ- %%-- (n2 + £2) gag: - 283%] 3.41 = (£2 + r12)2 - 84 (52 + mi) 34 32: W... a. 4n—-——— swam-.4 3.42 l 33 8112; L____I _(€2+n2)2 am 2 2 824 64 VVEn - 8h [E45 5E”'(” + g ) 5E?-+ 2n1——--] 3.43 .44 O) = (£2 + n2)2 am 2 2 32m am— Avat 9n 4g—EM” +€)anaa+2”‘ét' “44> . 4128; #1: [an a - 621.121 aw - 2. 11] 3,45 The other components of this dyadic are zero. The rate of deformation tensor is evaluated in the usual way as D =-% [Vv + WT] 3.46 The components of D are easily stated using (3.41) through (3.45) as on 8n at n .(42+n1[_4nap_.1§2__+_ni1 8;: 3 2 - (n2 + 521-3331” - 24 3%] 3.47 = = (52 "' 02V 2 2 _3 ll} _ 1112 321p DnE DEn 16h (n + E ) (557— n 8n '55?) ,6n%-5gg_g] 3.48 - (£2 '1’ T12)2 3111 2 2 3252 3111 ”t: ' 8n 4F 8n + (n I 5 ) aaan + 2” 55' 3'49 D (52 + n2)2 2 .22 -.1n3;:_§il.§_ - 25 39. 3 50 44 8n n 85 n a: an ° Again, the remaining components of Q are zero. Note that as would be expected tr(D) = 0 for incompressible flow. Finally the original goal of this development was to evaluate the integral vff 13-234 3.5. which in terms of (n,g) coordinates is 211 1 0° {[ .l’ {T Dnna + 0543 I D443 I 2”rm Dug"3 +4055 DnEZ -a _ " fil—zm. 3.3.3.. 21 This concludes the development of the equations that model the flow field. As can be seen there are numerous calculations. The next chapter discusses how these calculations were carried out. CHAPTER IV CALCULATION SCHEME Chapter III gave a derivation of the expressions (velocity profiles, stress distributions, etc.) which model the flow field around two contacting spheres translating in creeping motion through a viscous fluid. The present chapter demonstrates how these expres- sions were used to generate numerical values for velocity profiles, stress distributions, etc. The derived flow field modeling expressions in Chapter III all incorporated the stream function (3.20) and its derivatives (3.36) through (3.40). Given a n and g, calculating the stress function and its derivatives translates to numerically evaluating the integrals (3.28) through (3.33). These integrals contain four variables--A, B, C, and D--which are functions of s the integration variable. A, B, C, and D are determined by simultaneously solving (3.24) through (3.27) and itis the numerical evaluation of these quantities that will be dealt with first. By definition cosh a = (ea + e'a)/2; sinh a = (ea - e'a)/2 4.1 a -a For a > 2.3, ea + e' and ea - e are equal to ea within 1 percent error, therefore, 4.1 can be written as 22 23 cosh a = sinh a = ea/2 for a > 2.3 4.2 Using this fact with min(s,5a) greater than 2.3, (3.24) through (3.27) can be solved algebraically for A, B, C, and 0 resulting in A = 2. [-.-2s (2 + ,1. + 5n + .24-14.2 + g, + 23)] B - ~25 [e'25(2 +‘%2 +-§) + 6-2a5(302 + é%'+ %?11 C = 25 [e'25(2 +-%) + e-ZQS(Za +-%)] o = 25 [e'25(2 + %) -e’2°‘5(2a + $1 4.3 For min(s,sa) less than 2.3, (3.24) through (3.27) can be solved for A, B, C, and D by using a standard equation solving technique incor- porating matrix inversion. Substituting the analystical expressions (4.3) for A, B, C, and D, and the exponential approximation cfi’ cosh and sinh into the integrals (3.28) through (3.33), and then using the following relations for Bessel function integrals taken from math tables (I) j’ -at v _ (2b)vT(V+%) e Jv(bt)t dt - o (a2+b2)"*i/E -at v j. e Jv(bt) t dt = 4.4 o (a2 + 1.2)v + 3/2 I? 24 yields the following expressions for integrals (3.28) through (3.33). Remember these expressions only hold for min(s,sd) greater than 2.3 = 4(E-l)n _ 2[(C2 + n2)é - C] I 2(3- 21m2 + .2)! - c1 _ 4a (a +91. n(¥+ hi w2+¥fifl + _ 21(92 + 712)g ' 9] _ 2(23 +5) [193. n 2); ’ 91 4.5 MeUI <2 3?” + N A m I N v r—i N n N V Nir- + N (A) \ N I NIH |_1 I A N Q + m V r—1 N to N v w + 25 3 I3 = 4(4 - 111 '9” 215" (c2 + n2)5/2 (c2 + n 217/2] - 2[ 'D 1 + ZIICZ +302)é - C] (C2 + n 2)3]2 n(c2 + n2)i n +2(€ _ 2)[ 'C _ 2C - 3C1] +_2_] n(c2 + n2)3/2 n3(c2 + n2)1 (c2 +n 2)5/2 n3 '3 -9h 15 “401(01 '1' E) E ‘5 '1' Tl 62+F)” mz+n 2)7/2 J -n 2 2 I - g] -2[ - 1 + 2119 + n I (92 + nn2)3/2 n(92 + n2)1 n3 J -2(2a+-a) 1 '9 '29 '39 n(92 + n213/2 M3(g + n 211 (92 + n215/2 ,.j% J 4.7 O = -3(€- 1) (E - 2) 1 I4 4”1 (c2 + n2)5/2 I (c2 + n2)3/21 2, -c C(E ~212 I + rII: 2)2 + (C2 + n2)§72] _ 4an [ 3(0 + 51(5 + 20) + 1 (92 + n2)5’2 (92 + n 1 2)3/2 3 _ 2 '9 9 ] 2)§T+ (92 + O 2)3/2 —- 4.8 n (92 + n 26 52 I = -12( - 1 5 ’ 2 V( '2) 5 5 )5 (c2 + n2)5/2 (:2 + r12)772] 2 1 3 + 4[ - ’1 ] (c2 + n2)3/2 (c2 + n2)5/2 + 2(5 ‘ 2)[ 2 2 12; + 2CZ2 2 3/2 "“%"“2 3 2 nu +n) nk +n) « +n)/ 2 2 23c 3 + ] - 4a[ 2” we +¥§n m2+¥fifl (gH+n)3” (g2 + n 2W? (s:2 ‘+ 137/2 2 1 g 1 - 29[2 - + n 2(92 + n2)éf n2(92 + n2)3/2 (92 + n2)3/2 __§g(g + 2a) ] 4 9 (92 + n2)5/2 _ -3(2g - 3) 15 c2(g- 1) 3c I - 4n[ + + J 6 (c2 + n2)5/2 (C2 + n2)7/2 (c2 + )572 2 4 2 1 4c 3c + —-[ ~ + ] n (C2 + n2); (C2 + n2)5/2 (C2 + n2)5/2 -3(2§ + 3g)! 1592§a + §1_ 39 -401} [ + - (92 + n2)5/2 (92 + 2)7/2 (92 + 2)S/2 -Z[ 4 +___4£3_2_:§__—.394 mm n (92 + n2)5 (92 + n )3/2 (92 + n2)5]? where c = 2 - g and g = 2a + g. 27 With the above analytical expressions (4.5) through (4.10), we can numerically evaluate the integrals (3.28) through (3.33) by expressing these integrals as follows X I u a1 00 62 - F . Ia - ‘fl Fa - Fa approximate ds +‘/. a approx1mate ds 0 0 a = 1, 2, 3, . . . , 6 where F is the expression inside the integral Ia in (3.28) through (3.33) a approximate 15 Fa expressed with the exponential approximation of cosh and sinh discussed above Therefore, the integral Ia1 is just the correction for the difference between the actual integrals (3.28) through (3.33) and the approxi- mation of these integrals (4.5) through (4.10) in the range min(s,so) less than 2.3, where the approximation does not hold. Xu is chosen so that min(s,sa) is greater than 2.3, and the integral, Ial’ is numerically evaluated by using Chebyshev quarature. The integral Ia2 is replaced with the appropriate analytical expression in (4.5) through (4.10) which is easy to calculate, thus yielding numerical results for the integral Ia' The expressions for velocity profiles, stress distribution, etc., derived in Chapter III were ultimately all functions of inte- grals (3.28) through (3.33), n, and g; therefore, it is now possible to find numerical results for these quantities. The question now is how reliable are these numerical results? The remaining chapters will deal with this question. CHAPTER V ACCURACY OF CALCULATIONS The previous chapters have dealt with the derivation of the various expressions for modeling flow around two contacting spheres, and how these expressions could be simplified for better implementa- tion on a computer. This chapter will deal with how the numerical results compare with previous work, where there might be errors in the calculations, and what was done to eliminate error as much as possible. The calculations were done in double precision on the Control Data Corporation Cyber 760 computer at Michigan State University. Double precision on the Cyber 760 allows for accuracy to approxi- mately the 29th decimal place according to operating manuals. Two different sets of calculations were done: one for the flow field around two contacting spheres of equal radius (this will be referred to as the a = 1 case in the following); the other for the flow field around two contacting spheres, with the trailing sphere having five times larger radius than the leading spheres (the a = .2 case). It was emphasized in Chapter IV that the flow field modeling expressions were all functions of the integrals (3.28) through (3.33), n, g, and 0:; therefore, it is obvious that any error in numerical 28 29 values calculated from the flow field modeling expressions is a result of error in evaluating the integrals (3.28) through (3.33). Recall from (4.11) (shown here for reference) x 00 Ia1 Ia2 I3 = Jr Fa - Fa approximate ds + Jf Fa approXimate ds 4.11 0 a = 1, 2, . . ., 6 O that each one of the integrals (3.28) through (3.33) could be repre- sented as the sum of an integral that needs to be evaluated analyti- cally, Ial’ and an integral that can be evaluated analytically, Iaz Each one of the Ia1 and 1&2 integrals has an error associated with it, and it is these errors that we shall consider first. The error in the Ia introduced by the IaZS in (4.11) are due to the differences between Fa and F in the region 2.3 = a approximate min(s,5a) 9o mgaoucou--.w.m 69:99; 3. .mmmu N. n a com ANonV L> mo mgaoucou--.m.m mgzmwm 50 51 .mmmu H u a Lo9 N> 9o mcaoucouu-.oH.m mcamwm .mmmo N. u a cow N> 9o mczoucou--.HH.m mgamwd 52 53 09H 02 .mmmu N. 03 _ a Let H om _ N um g cm H .m> >--.NH.m atamwa 09 H ow o.H 54 Rate of Deformation Tensor andTVorticity In this section, the contour plots of the various components of the rate of deformation tensor will be considered together with vorticity contours to facilitate discussion of the identification of stagnant, stretch dominated, rotation dominated, and shear dominated regions of the fluid for the a = 1 and the a = .2 cases. The rate of deformation tensor, 0, (D = 1/2 (v v + VvT)) is considered in r,z coordinates because these coordinates are more easily understood physically than the n,g coordinates. D is expressed more explicitly as 2 = Drr Dr¢ Drz D¢r D¢¢ D¢z Dzr Dzo Dzz __ = dvr/dr 0 $(dvr/dz + dvz/dr 0 vr/r 0 £(dvr/dz + de/dr) 0 dvz/dz _ 6.1 Remember that the r,z components of D can be expressed in terms of n,E coordinates by the methods mentioned earlier, and this is what was done to evaluate these components. Contour plots of Drr’ D zz and Dr2 for the a = 1 case are shown in Figures 6.13, 6.15, and 6.17, respectively, and contour plots of Drr’ D , and Drz for the 22 a = .2 case are shown in Figures 6.14, 6.16, and 6.18, respectively. 55 .mmmo H u a com Amonv LLo pm menopcou--.mH.m 693mm; om om: mHn .mmmu N. u a Low Amoncho we mgaoucoo--.vH.m mgamwu 56 0/5- \O/ o 57 .mmmu H u a com AmonvNNo mo mgaoucou--.mH.o mgamwg om ooH 0H OCH: om: .mmmu N. u a 909 AmonVNNQ 9o mgaoucounu.cH.m mesm9u oml K\\ om: \\\|Ili om .mmmu H a :09 ANonVNLo mo mesoucou--.NH.o mgamwu 59 .mmmu N. n a cow ANonvNLo we menopcou--.mH.m wcamwd 61 Vorticity is defined as .n.=va=(-a-z-'—‘--E%) 6.2 and can be viewed physically as a numerical representation of the amount of spin the fluid is imparting. Contour plots of vorticity for the a = a case, the a = .2 case, and the single sphere case are displayed in Figures 6.19, 6.20, and 6.21, respectively. In the contour plots of vorticity that are given, positive vorticity implies that the fluid would impart spin in the counter clockwise direction and vice versa for negative vorticity. To illustrate this more clearly, if a flat blade whose axis of rotation was perpendicular to the direction of flow were placed in a region of fluid for which vorticity was positive, the blade would spin in the counter-clockwise direction. Since for axisymmetric flow vorticity is independent of coordinate system used to determine it, vorticity was calculated here in n,g coordinates. To check the accuracy of 2, dvz/dz, dvr/dz, dvz/dr, and dvr/dr were calculated at different points in the flow region and compared to the velocity contours to see if they had reasonable values. It was found when checking against values of V2 for r less than .1 that error resulted in the evaluation of dvz/dr, therefore, de/dr was evaluated by using values of V2 at low values of r in dvz/dr = sz/Ar. All other components of D seemed to have reason-' able values. .mmmu H u a com HNonv xuwowugo> yo mgzoucou--.mH.m mgammu 62 .mmmo N. u a com ANonV xumuwugo> 9o menopcou--.oN.o mczmpm 63 oH 64 mNoo.-. .mmmu mcmzam mecvm Low mumuwpgo> 9o mcaoucou--.HN.w mgmmwm 65 It was found that using the above correction of dvz/dr for r less than .1 that vorticity could be evaluated accurately through- out the flow field also. Now, the Drr and D22 contour plots will be considered together because they have many of the same characteristics. Physi- cally, Drr can be viewed as stretching normal to the z axis, and 022 can be viewed as stretching normal to the r axis. Both Drr and D22 reach absolute maximums at the front and back edges of the spheres for the a = 1 case and the a = .2 case, and both Drr and D22 have the characteristic that their absolute magnitude decreases much more rapidly proceeding in the r direction away from the front and back regions of the spheres than proceeding in the z direction away from these regions for both the a = 1 case and a = .2 case. Next, the contour plots of the shear component of the rate of deformation tensor, 0 , will be considered. Note that this is a [‘2 cylindrical coordinate shear rate. For the a = 1 case, 0 , reaches rz equivalent maximums at the equators of the two spheres, while for the a = .2 case only one maximum, at the equator of the trailing sphere, is reached. Drz decreases rapidly when proceeding in the r direction away from the maximum points and when proceeding around the spheres toward the poles, while Drz decreases least rapidly pro- ceeding at approximately a 45°angle away from the maximum points. Much like Drz’ vorticity reaches identical maximums at the, equator of each sphere for the a = 1 case and at the equator of the trailing sphere for the a = .2 case. When comparing the a = .2 case 66 to the single sphere case demonstrates that the vorticity contours for both cases are almost identical with the smaller leading sphere in the a = .2 case having some effect on the contours in the region around it. It was stated earlier that the rate of deformation and vorticity contour plots help identify stagnant, stretch dominated, rotation dominated, and shear dominated regions of the fluid. The criteria for identifying each one of these regions is given in Table 6.1. Table 6.1.-—Criteria for identifying stagnant, stretch dominated, ’ rotation dominated, and shear dominated regions of fluid Magnitude State of Fluid weak weak weak Stagnant weak weak strong Stretch dominated strong weak weak Rotation dominated strong strong weak Shear dominated Using the criteria of Table 6.1 and the rate of deformation and vorticity plots, the stagnant, stretch dominated, rotation dominated, and shear dominated regions of the fluid were identified for the a = 1 case and a = .2 case and are shown in Figure 6.22. It is left to the reader to decide where the boundaries to these regions are. 67 III III II II Figure 6.22.--Stagnant (I), stretch dominated (II), rotation dominated (III), and shear dominated (IV) regions of fluid for a = 1 and a = .2 cases. 68 Surface Shear Rate One of the purposes of this work, as mentioned in Chapter I, was to gain a better understanding of Jefri's experimental results. Following the reasoning in Chapter I, it would appear that Jefri, in contrast to other workers, has found a way to measure the effect of elasticity on drag, for viscoelastic fluids, when the shear modulus is constant, by using settling, contacting spheres experiments. This argument is based on the fact that drag reduction occurs at an aver- age surface shear rate which is low enough that the shear modulus is not dependent on shear rate in these experiments. Jefri arbitrarily has defined an average surface shear rate as 9 = vm/R(1 + k) 6'3 and if this average shear rate is not correct, then the conclusions drawn from Jefri's experimental data are not correct. The calcula- tions done here allow us to determine what the actual average surface shear rate is. The surface shear rates considered so far have been for Newtonian fluids, and although they may not be precisely correct for the viscoelastic fluids considered here, they are probably very close. The dimensional tangential surface shear rate for the con- tacting spheres is Ysphere = -20€n(v“/R) ~ g = 1, - a 6.4 69 The average surface shear rate is determined by integrating the 2 component of the tangetial surface shear rate over the surface of the spheres and dividing by the surface area. The resulting expres— sion for an individual sphere is 2n n i = (I/4n R2) If Jf -20€q((vm/R sine) stinO dado) 6.5 O 0 where O is angle measured from the positively oriented z axis. Using D as calculated here and numerical integration in (6.5) the true average surface shear rate for each sphere in the a = 1 case and the a = .2 case was calculated and is given in Table 6.2. Jefri's average surface shear rates are included for comparison. Table 6.2.--Average surface shear rates Average Surface Shear Average Surface Shear Sphere Rate--Jefri Rate--Actual Leading and Trailing sphere a = 1 case .5vm/R '6766 vw/R Trailing Sphere ' a = .2 case .1667mv/R .1962 Vm/R Leading sphere a~= .2 case .1667 vm/R .1160 vm/R As can be seen in Table 6.2 the actual average surface shear rates are slightly larger than the average shear rate Jefri calculated, but 70 still confirm that Jefri's experiments were done at shear rates where shear modulus is not shear rate dependent. It was mentioned in Chapter I that prior to Jefri's work, Chhabra had done settling experiments with single Spheres in the same types of fluid as Jeffri used, and both Jefri and Chhabra used average surface shear rates to determine if the shear rates were in the region where shear modulus was shear rate dependent or not. It would seem that the maximum surface shear rate should be considered instead of the average surface shear rate. To better understand this, recall that Chhabra's average surface shear rate was defined as 5 = Vm/R 6.6 The tangential surface shear rate for Chhabra is given by v Ysphere = 3/2 -—-Sine 6.7 R while Jefri's tangential shear rate is given in (6.4). In both cases there exists the possibility that the average surface shear is low enough so that the shear modulus is not shear rate dependent, while regions of the fluid are at shear rates where the shear modulus is shear rate dependent. If the maximum surface shear rate were low enough that shear modulus is not shear rate dependent on any part of the sphere's surface, then obviously the effect of changing shear modulus can be ignored. ' The calculations done here allow us to determine if changing shear modulus effects can be ignored in Jefri's experiments. The 71 plots of dimensionless tangential surface shear rates for the a = 1 case, the a = .2 case, and the single sphere case are given in Figures 6.23, 6.24, 6.25, and 6.26, where dimensionless shear rate is defined as 4 /(vm/R) = -2D Since it has been noted that the in' o = .2 case is much like a single sphere, the single sphere tangen- tial surface shear rate shown in the plots is scaled to be on the same dimension as the trailing sphere in the a = .2 case. The plots show that the tangential surface shear rates on a solitary sphere .2 are almost identical to those on the trailing sphere for the a case with the exception being in the polar regions as expected. The plots were used to determine the maximum surface shear rates and the results are given in Figure 6.27. Using the maximum tangential surface shear rates as shown in Figure 6.17 and Jefri's experimental data shows that the maximum shear rates for the a = 1 case are not in the region where shear modulus is shear rate dependent, but the a = .2 case has maximum shear rates that are in the region where shear modulus is shear rate dependent. Attempt to Model Drag In Chapter II, Jeffri's suggested method for calculating the drag for the contacting spheres in his settling experiments is shown. This drag is expressed as i . F1 = l’Fo + we (1.F1) 2.1 .mmmu H u a Low mum; Lamzm wumwgzm Hm9ucmmcmu mmchowmcmewa--.mN.m wcsmwm NH o.N m.H m.H e.H N.H o.H w. m. e. N. o 72 o.H cw cm. «H 73 .mmmu N. u a :9 mgmsam mcwummH cow mums gmwgm mom9g=m meucmmcmu mmch09m=mswa--.¢N.o «gnaw; o.N m.H o.H ¢.H N.H o.H m. w. v. N. _ _ _ _ H _ _ . H _ N - c woN- 74 .mmwu N. n a :9 mcmcam acmemcu 909 mum; cmmsm mummcsm mepcmmcmp mmchowmcmewoii.mN.m 693mm; 0H m m N o m 9 m N H o 4 _ A _ _ _ _ _ _ d 0 am. 75 A.mmmu N. u a :9 mgmcgm mew—9mg“ mm m~9m mean m9 mcmsam 659V .mmmu mcmsam mpmcwm co; wan; gmmgm momwcsm Hmwpcwmcmu mmmpconcmswoui.oN.o mgzmwm v m N H o H. N: mu 9. mu _ _ H _ . 11 _ . _ c NON- 1.3 Figure 6.27.--Values and locations of maximum dimensionless surface shear rates (4 xhog/R) for single sphere case, a = 1 casg§ and a = .2 case. 77 1 ° F1 = '2UVco Jf- 9-9 : D dV 2.3 In order to use (2.1) to model drag, one must be able to evaluate the where integral in (2.3) accurately. This integral can be expressed in n,g coordinates, using (3.52), as the sum of two integrals O “3 Ibl -. = _ 3 . 3 3 2 1 F1 4npvm j, j’ [Dnn‘ + D55 + P¢¢ + 2DnnDnE -a 0 2 + 4 DEEDnE'] (83 2)3 dnda E + n ‘ 1 00 I b2 _ 3 3 3 2 4Wv°° f I who + DEE; + D¢¢ + 2DrinDnE o o + 8n dndg 40550:: 3 (£2 + n2)3 6.8 The n,g coordinates are an inverse coordinate system so the n and g limits of O in (6.8) actually represent infinity in the flow field. The zero limits for n and g in (6.8) have to be replaced because using a zero limit for n and 5 would involve dividing by. zero in numerous places during the calculation of 2. Care has to be taken when using a numerical integration scheme that the values used to replace the zero n. E limits in (6.8) are sufficiently low to 78 incorporate the whole flow field (see Figure 3.1 to understand this more clearly). Trying to choose these low integration limits causes a dilemma because at lower n and g the calculations of D are more susceptible to round off error and the error magnification problem mentioned in Chapter V. The infinity limit for n in (6.8) also has to be replaced with some reasonable number to facilitate carrying out numerical integration. Following is a discussion of how the above limits were established, how integrals Ibl and Ib2 in (6.8) were numerically evaluated, and what the results were for the a = 1 case, the a = .5 case (trailing sphere with radius two times readius of leading sphere), and the a = .2 case. The numerical integration of Ib1 and Ib2 was done by using a two dimensional Simpson's Rule which given x integration limits of a and b with midpoint m and y integration limits of c and d with midpoint n is given by ~ d f f(x.y)dydx = hk/36 [f(a,c) + f(a,d) + f(b,c) + f(b,d) + C {BRET 4(f(a.h) + f(m.c) + f(b.n) + f(m.c)) + 16f(m.n)] where h = b - a and k = d - c. A grid of n, 5 points was developed for each case considered so that using the formula (6.9) successively over the grid the integrals Ib1 and Ib2 could be evaluated. The actual n,g quadrature points used for the a = .2 case are given in Table 6.3. Graph 3.1 79 Table 6.3.--Quadrature points for the a = .2 case n E for IbI g for Ib2 .0001 -.2 .0001 .01 -.175 .0075 .02 -.15 .015 .03 -.125 .0225 .04 -.1 .03 .05 -.09 .04 .06 -.08 .05 .07 -.065 .065 .08 -.05 .08 .09 -.04 .09 .1 -.03 .1 .125 -.0225 .125 .5 -.015 .15 .2 -.0075 .175 .25 -.0001 .2 .375 .3 ,5 .4 1.0 .7 1.5 1.0 80 shows that the quadrature points in Table 6.3 cover the flow field very well for the a = .2 case. Similar quadrature points were used for the a = 1 case and the a = .5 case. The upper limit of 1.1 for n was chosen for each case because the integrands in Ib1 and Ib2 are less than or equal to 5 x 10’5 for n greater than or equal to 1.1. It was noticed from previous calcu- lations that for low g values the calculations of 2 were accurate, but for n values such that r is less than .1 the calculations of D were not accurate at all. Some of the quadrature points shown in~ Table 6.3 fall into this low r region, but the value of the integrand is in line with values of the integand ”at7 points nearby which are not in this region. It appears by observation of points that are accurate near the r axis that the integrand falls to zero as one approaches the r axis, and it turns out that the calculated integrand in the low r valued region is very low. It can be shown any inaccuracy of the integrand in the low r region as calculated here creates less than a .1% error in the value of the integrals Ib1 and Ib2' The error term associated with the Simpson's Rule integra- tion scheme used was evaluated, and the maximum error in the value of the integrals IbI and Ib2’ with the Quadrature points used, was less than .O4%. The results of the integration for the three cases considered are shown in Table 6.4 in terms if j - F1. Also shown in Table 6.4 is the Newtonian drag FO as calculated by Cooley and O'Neill. Recall that the drag is being modeled by 81 Table 6.4.--Fo and F1 for a = 1, .5, and .2 a -Fo/6n v0° -F1/6n v0° 1 1.29 0 .5 2.09 -.0009 .2 10.00 -.0024 l o F = i 0 F0 + we (1 0 F1) 2.1 where We is small. As can be seen from Table 6.4 the proposed drag model (2.1) for the contacting spheres settling in a viscoelastic fluids does predict a reduction in drag, but the magnitude of the drag reduction predicted is not close at all to the drag reduction Jefri noticed in his experiments. In fact, the integrals Ib1 and Ib2 have values on the order of .2, and with the error of at least .2% in the calculations involved in calculating F1, F1 is zero for all intents and purposes. CHAPTER VII CONCLUSIONS AND RECOMMENDATIONS In the first chapter, it was stated that the objective of this work was to obtain an accurate numerical representation of the flow field around two contacting spheres settling, in creeping motion, in a Newtonian fluid. We have demonstrated a method here in which the stream function, velocities, vorticity, and rate of deforma- tion tensor can be calculated accurately throughout the whole flow field. The average surface shear rates and maximum surface shear rates were established for the settling contacting spheres here. The maximum surface shear rates show that, contrary to what was originally believed, constant shear modulus cannot be assumed for Jefri's (1984) settling contacting spheres in viscoelastic fluid experiments. However, the assumption that Jefri has shown drag reduction that is due solely to elasticity effects still remains valid. The second objective of this work was to test how well a first order perturbation expansion of the Newtonian drag in terms of small weissenberg Number fit Jefri's experimental results. It was found that this method of modeling drag shows virtually no drag reduction and did not come close to matching the magnitude of drag reduction found in Jefri's experiments. It is obvious that another 82 83 method of modeling the drag in Jefri's experiments must be hypothe- sized. One such method is hypothesized below. It has been shown by other workers (e.g., Morrison and Harper, 1965) that changing the no slip boundary condition to a slip boundary condition results in drag reduction for non-Newtonian fluids. There- fore, it is recommended that the next attempt at modeling the drag reduction in Jefri's experiments should be to change the no slip boundary condition to a slip boundary condition and see if the behavior of the calculated drag is comparable to the experimental results. The results of the work here would be valuable in estab- lishing the slip boundary condition because the slip would most probably be proportional to surface shear stress. BIBLIOGRAPHY 84 BIBLIOGRAPHY Cooley, M. D. H., and O'Neill, M. E. J. Camb. Phil. Soc. 66 (1969): 407-415. Chhabara, R. P., Uhlherr, P. H. T., and Roger, D. V. J. Non-Newtonian Fluid Mech. 6 (1989): 187-199. Davis, A. M. J.; O'Neill, M. E.; Dorrepal, J. M.; and Ranger, K. B. J. Fluid Mech. 77 (1976): 625-544. Happel, J., and Brenner, H. Low Reynolds Number Hydrodynamics. Englewood Cliffs, N.J.: Prentice-Hall, 1976. Jefri, Mohammad Amin. “Translation of Two Contacting Spheres in a Viscoelastic Fluid." Ph.D. dissertation, Michigan State University, 1983. Morrison, S. R., and Harper, J. C. I + EC Fund 4 (1965): 176-181. Reed, L. P., and Morrison, F. A. Int. J. Multiphase Flow 1 (1974): 573-584. Stimson, M., and Jeffrey, G. B. Proc. R. Soc., A111 (1926): 116-119. Illl) 168 9122 H I“ Y. H "I II n H I | H || | 1293 O3 3 TIHIHHIIH