Willhll'iillillllIIMW‘l u‘ 109 879 THS um mm m lflfllflllllllmwllllzlllwlfll ”'5" 3 12 This is to certify that the thesis entitled Fluid Flow Through a Partially Filled Cylinder presented by Raymond N. Greenwell has been accepted towards fulfillment of the requirements for Ph. D. Applied Mathematics degree in W707 Major professor Date June 8, 1979 0-7639 OVERDUE FINES ARE 25¢ PER DA PER ITEM ‘ Return to book drop to remove ‘ this checkout from your record. FLUID FLOW THROUGH A PARTIALLY FILLED CYLINDER BY Raymond N. Greenwell A DISSERTATION submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics 1979 ABSTRACT FLUID FLOW THROUGH A PARTIALLY FILLED CYLINDER By Raymond N. Greenwell The laminar flow of viscous incompressible fluid through a partially filled circular cylinder is studied.’ For this problem the Navier-Stokes equations can be transformed into Laplace's equation, and the region of solution is the interior of two intersecting circles. The region is mapped to the unit circle by a conformal transformation, so that the problem may be solved by the Poisson integral formula. The integration is performed by a combination of residues and numerical integration. It is found that by not filling the cylinder all the way, one can increase the flow rate by as much as 27.3%. Two special cases are also studied, using different methods. The first is the case where the region of solu- tion is the interior of two circles intersecting at right angles. In this case the computations simplify and the integration can be performed analytically. The second case is the nearly half filled cylinder. A perturbation method is used on the well known solution of the half filled circular cylinder. The results obtained from these special cases agree with the previously obtained results. DEDICATION This thesis is dedicated to my dearest wife, Karla, without whose love, support, and encouragement this project would have been completed a lot sooner. ii ACKNOWLEDGMENTS I wish to express my gratitude to Professor C.Y. Wang, who suggested this problem to me and Who has provided invaluable assistance and guidance. I also wish to thank Professor Larry Segerlind for verifying my results with his finite elements program, and Dr. James Lyness for his suggestions on numerical integration. Finally, I am indebted to Mrs. Mary Reynolds for typing this manuscript, and to Professors NOrman Hills Chi Y. Lo, Glen Anderson, and David Yen for acting on my thesis committee. iii TABLE OF CONTENTS Chapter Page I INTRODUCTION 1 II FORMULATION OF THE PROBLEM 4 III SOLUTIONS OF SPECIAL CASES AND RELATED PROBLEMS 12 IV SOLVING THE BOUNDARY VALUE PROBLEM BY CONFORMAL TRANSFORMS 15 V CALCULATING THE FLOW 28 VI INTEGRATION TECHNIQUES 31 VII RESULTS 41 VIII THE PARTIALLY FILLED CYLINDER WITH a = ——1— 47 fl IX PERTURBATION SOLUTION OF THE HALF FILLED CYLINDER 64 x DISCUSSION 77 APPENDIX: ANGLE OF INCLINATION FOR LAMINAR FLOW 79 BIBLIOGRAPHY 81 iv Figure 1.1 LIST OF FIGURES Flow through a partially filled circular cylinder The effect of surface tension Body forces on the fluid The region of solution Flow through two parallel plates Cross section of the partially filled cylinder 2 +./{-a2 2 The transformation wl = - -a z The transformation w2 = w? The transformation w3 = cw w3+-1 The transformation w = w3-l The transformation given by equation 4.6 The curve for integrating equations 6.17 and 6.18 FIOW’(Q) vs. fill level (a) Velocity profile for a = .707 Velocity profile for a full cylinder . 2 l+-It The transformation 0 = —————§ l-it Page 10 13 15 17 18 19 20 36 43 44 46 SO Figure Page 9.1 Flow through the almost full circular cylinder 65 9.2 Cross section for the almost half filled cylinder 66 9.3 Flow through the almost half filled unit cylinder 74 9.4 Plotting the flow using perturbation methods 76 vi CHAPTER I INTRODUCTION Let us consider the steady flow of a viscous, incompressible, Newtonian fluid through a partially filled inclined circular cylinder (see Figure 1.1). we will assume that the flow is laminar, that the cylinder is infinitely long, and that surface tension can be ignored; these assumptions will be examined more closely in the next chapter. When the pipe is full, this type of flow is called pipe Poiseuille flow or proper Poiseuille flow [1]. Figure 1.1. Flow through a partially filled circular cylinder We wish to consider the following question: for a given inclination, to what level should the pipe be filled in order to maximize the flow rate? At first glance this question may seem trivial; one would think that filling the cylinder all the way would maximize the flow. However, this is not true. The velocity of a viscous fluid is zero on the walls of the pipe. Even though a partially filled cylinder has a smaller cross section of fluid, the fluid may be flowing faster due to a shorter solid boundary where the velocity is restricted to zero. In other words, the fluid on the free surface is not adhered to any walls, and it may flow faster than it would in a full cylinder. The increased velocity may be more than sufficient to compensate for the loss in cross sectional area, and thus the flow will be greater than for a full cylinder. Civil engineers have known experimentally for some time that turbulent flow through a circular cylinder is maximized when the cylinder is partially filled. In 1946 Camp [2] graphed the flow as a function of the fill level. These graphs were reproduced by Chow [3], who calculated that the depth for maximum flow was .97 times the diameter, where the flow is about 3% more than for the full cylinder. An increase of 3% is not very significant, but we might expect a greater increase for laminar flow, since turbulent flow attains near maximum velocity mudh closer to the walls than laminar flow. Therefore, the effect of a free surface should be felt more strongly in laminar flow. Chemical engineers have analyzed this problem by methods different from those used in this thesis. Chaudhury [4] and Buffham [5] used bipolar coordinates and Fourier integrals to study flow in Open circular channels, but they did not search for the maximum.flow rate. Charles and Redberger [6], Gemmell and Epstein [7], and Yu and Sparrow [8] have worked on the problem of maximizing oil flow through a pipe by partially filling the pipe with water. The solution for a full or half filled cylinder has long been known; it is rederived in Chapter III. The solution for a partially filled but not half filled cylinder is much more difficult, since axial symmetry is no longer present. We will solve the problem by means of conformal transformations and techniques of numerical analysis. Two verifications will also be performed; the first is the solution of a particular case in which the computations can be done by a different method, and the second is a per- turbation solution of the half filled cylinder problem. CHAPTER II FORMULATION OF THE PROBLEM The critical parameter in determining whether flow is laminar or turbulent is the Reynolds number LV.’ \J 2.1 R = where E is some length which characterizes the problem, V is some characteristic velocity, and v is the kinematic viscosity of the fluid. For the flow through a pipe, V is taken to be the average velocity and L is taken to be four times the area of the cross section divided by the perimeter (for a full circular cylinder, this is just the diameter ) [9]. The flow through a full cylinder is found experimentally to always be laminar when R < 2320 [10], though the flow may remain laminar for values of R as high as 24,000 if the pipe is very smooth and care is taken to avoid disturbances [11]. We can see from equation 2.1 that the flow will be laminar if the viscosity is high enough, the pipe is small enough, and the average velocity is low enough. For instance, if crude oil, with a typical viscosity v = .5 cmz/Sec, were to flow through a pipe with a diameter of 50 cm at an average velocity of 20 cm/Sec, then R = 2000 and the flow would be laminar. We assume that the cylinder is infinitely long. In practice, this means that the cylinder is long enough that end effects can be ignored. For a full pipe, this assumption is valid when the fluid has travelled a distance of approximately .06R times the radius of the pipe [12] . We also wish to ignore the effect of surface tension, Which tends to make the free surface curved rather than flat (see Figure 2.1). The height 11 that the fluid rises - 1:33 h Figure 2.1. The effect of surface tension when it comes into contact with the surface is [13] 2v (l-sin 6) (2.2) h = / F5 P9 YSA = Yrs + YAF °°s 9 where YSA = surface tension between the solid and the atmosphere YFS = surface tension between the fluid and the solid YAF = surface tension between the atmosphere and the fluid 9 = angle of contact between the solid and the fluid 9 = acceleration of gravity = 981 cm/sec2 p = density of fluid. The values of y and p are measurable physical constants. For example, at 20°C the surface tension between water and air is = 72.8 dyn/cm. The density of water is VA]? .998 gm/cm3. If we assume the worst and let 9 = 180° (actually a is usually much smaller for water), then equation 2.2 gives h = .5 cm If the diameter of the cylinder is much larger than h, then the effect of sur- face tension can be ignored. Let us set up a Cartesian coordinate system with the z axis parallel to the axis of the cylinder and the y axis perpendicular to the free surface (see Figure 1.1). The notion of the fluid is described by the Navier-Stokes equations for viscous, incompressible, time independent flow [14]. In Cartesian coordinates they are: 2 3 a v. 3 3v. 1 hp I 2.3 0 = v 2) -—-—$—-- - - Z) v ———-+ F. k=l axka p axi k=1 k axk I I = 1,2,3 where x1 = x, x2 = y, x3 = 2 vi = ith component of the velocity p = pressure exerted on the fluid i = ith component of the exterior body force p = density of the fluid V = kinematic viscosity of the fluid When the flow is due to gravity, there exist body forces F2 = -g cos a in the y direction and F3 = g sin a in the z direction, Where g is the acceleration of gravity and a is the angle of inclination of the pipe (see Figure 2.2). Since the pipe is infinitely long, it is reasonable to assume that the velocities are independent of z. F3:9'S[n0( ,‘ F‘s-Sucosm Figure 2.2. Body forces on the fluid NOW consider the case where a = 0, that is, the pipe is on a level surface. In this case, the fluid would not flow at all, that is, v1 = v = v = 0. If we 2 3 now change a to a positive value, so that F3 is posi— tive, this does not affect equation 2.3 for i = l or 2. Thus, it is reasonable to assume that even in this case v1 = v2 = 0. Also, the pressure is independent of 2, since the atmosphere above it is at the same pressure at every point in the pipe, so %§’= 0. Since there is only one nonzero component of velocity, and since this component does not depend on 2, we will denote v3 = v(x,y). When we make all the simplifications in equation 2.3, the result is 2 2 2.4 o= v(i-‘21+§—‘2’-) +F3. ax By Define a non-dimensional velocity and non-dimensional lengths ’=—‘—‘-’l— ’=’.‘ =2 V 2 x L Y L' F3L Where L is the radius of the pipe. Then equation 2.4 becomes 2 , 2 , 2.5 L1...2.+§._V__2.= 1, 6(X’) 6(y’) This is a special case of Poisson's equation. we shall henceforth drOp the primes. NOtice that in our new coordinate system the pipe has unit radius. Also, the velocity is negative, since v is in the opposite direction of the original v. Part of the boundary conditions for viscous fluid is that v = 0 on the boundary between the fluid and the wall. If we make the change of variable 2 2 2.6 u = E—IIX’ - v then equation 2.5 becomes 1=_a_21,I2_v= aflxzwz- Him- 2 2 2 4 u 2‘ 4 u) ax 6y ax By 2 2 _ l §_s.+ 1 2.2. ’ 2 ‘ 2 2 ' 2 6X By 2 2 2.7 g—3'-+ §_B.= O 2 2 ax 6y with boundary condition 2 2 2.8 u = £—i%x- on the fluid-wall boundary. For the partially filled cylinder, there is also a boundary condition on the free surface, namely, that the stress tensor avi av. must be continuous at the boundary of the liquid and the atmosphere [15]. H is the viscosity of the fluid (not to u/P). The be confused with the kinematic viscosity v tangential stress tensor on the free surface y = constant has components 6v __1) 6x2 5v _2., T21 = "9521 + “(5x1 = 0 lO sz 5v3 6V T23 = -P623 + u(-a§-3-+ 5;? = us"; Since there is negligible tangential stress between the free surface and the atmosphere above it, we must have 21:0 by This condition can be achieved by solving equation 2.5 with v = 0 on the boundary of the region formed by reflecting a cross section of the partially filled pipe across the free surface (see Figure 2.3). Due to symmetry, the lines of constant velocity will be perpendicular to the free surface. line of constant velocity Figure 2.3. The region of solution 11 We thus solve equations 2.7 and 2.8 over the region in Figure 2.3. The solution over either the tOp half or the bottom.half of the figure can be taken as the solution to our problem. It should be noted that equation 2.7 is the Laplace equation. With the boundary values prescribed, the problem is known as the Dirichlet problem. The boundary condition given by equation 2.8 occurs in other circumstances. The solution to equations 2.7 and 2.8 also is the solution to the torsion problem [16]. Furthermore, it gives the deflection of a membrane, such as a soap film, stretched over the region [17]. CHAPTER III SOLUTIONS OF SPECIAL CASES AND RELATED PROBLEMS Several special cases of the problem posed in Chapter II have already been solved. Sokolnikoff and Sokolnikoff [l8] solved the torsion problem for a lune formed by two circular arcs of equal radius and intersecting in the interior at a right angle. Another stated problem is solved by‘Wigglesworth and Stevenson [19], who studied the problem of circles of unequal radii intersecting in the interior at a right angle. Lal [20] studied the related problem of flow through a tube whose cross section is formed by the intersection of two circles, where the area under consideration is outside one circle but inside the other. The solution to the full cylinder problem is well known. Let the z axis run through the center of the cylinder. Then equations 2.7 and 2.8 become 31 93.2,az_u=o ' 2 2 BX 6y u = % on the unit circle x2+y2 = l. 12 13 By InspectIon, the solution 15 u = %, or 3 2 v = x24-22 _ u = x24-y2 _ l _ r2-1 ° 4 4 4"4° This is the familiar parabolic velocity profile of pipe Poiseville flow. The total amount of flow per unit time is -w/8, found by integrating the velocity over the cross sectional area of the fluid, which is the unit circle in this case. In order to compare this value with the flow rate for the partially filled cylinder, denote 3.3 Q = -(flow rate)/(w/8). Thus Q = 1 for the full cylinder. The half filled cylinder is easily handled by using the arguments in Chapter II to conclude that the solution is just the bottom half of the solution for the full cylinder problem. Thus Q = % for the half filled cylinder. As an indication of the significance of the partially filled cylinder problem, let us examine a similar but simpler problem. Consider two-dimensional planar flow between two infinite plates, commonly called plane Poiseuille flow. Let the two plates be at y= 0 and y= -1 (see Figure 3.1). m :1? Figure 3.1. Flow through two parallel plates 14 The only non-zero component of the velocity is v1 = v(y). USing arguments similar to those in Chapter II, equation 2.3 reduces to 2 3.4 1321 = 1 dy with boundary conditions The solution is easily found to be Integrating this solution between -1 and 0 gives the flow per unit time per unit width OJLZY £230 -1 Q"Ja._12+2dy=6+4 -1='i’2' ° NOW suppose we raise the upper boundary an infinitesimal amount so that the fluid no longer touches it. We still have the boundary condition v(-l) = 0, but y = 0 is a free surface and needs a different boundary condition. USing the same arguments as in Chapter II, we consider this problem to be the bottom half of the flow between two plates at y = -l and y = 1. Thus our new boundary condition is v(l) = 0, and the solution is found to be 2 =L-l V 2 2 0 2 3 0 _ z__} _L_ =1 0‘1 2 2dY’6 2-1 3' Thus the flow is quadrupled when the fluid has a free surface at the t0p. This shows that the presence of a free surface can cause a tremendous increase in the rate of flow. CHAPTER IV SOLVING THE BOUNDARY VALUE PROBLEM BY CONFORMAL TRANSFORMS We now turn to the problem of the partially filled cylinder. Let a = the distance between the center of the cylinder and the free surface of the fluid, Where 0 < a < 1 corresponds to the more than half filled cylinder, and -l < a < 0 corresponds to the less than half filled cylinder (see Figure 4.1). we will solve the problem by a complex transformation from the region Figure 4.1 to the unit circle. Since a harmonic function I ’l I, \\ a h ’ y a; / llu.‘. u I - I (J ' ‘ \ l \ I I a>0 \ , a<0 \~ , ‘0. r.’ Figure 4.1. Cross section of the partially filled cylinder. a1 = sin-1a c2 = cos-la 15 16 (that is, a solution of the two-dimensional Laplace's equation) is still harmonic When the variables are transformed under a conformal mapping [21], and since we know the solution to the Dirichlet problem on the circle, we can use that solution to find the solution to the problem at hand. We compose a series of transformations (see Figures 402-405). _ z + “I -a2 40]. W1 "' 2 z - -a 4.2 w2 = w? where p = W -1 2(w-—cos a) with branch 0 g.e < 2W 4.3 w3 = cw2 Where c = cos e-ki sin e = cis e and 9 = v sin-1a w-cos-la -+1 4.4 w = :3__1 . 3 The bilinear transform given by equation 4.1 changes the two intersecting circles into intersecting lines (see Figure 4.2). The transform given by equation 4.2 changes the angle between the two lines to n by an appropriate choice of p. Since raising wl to the p power multiplies its argument by p, we set 1 l p(2w-—cos- a)-p cos- a = w CW] kn ) a 5‘“ cos'k \\‘\; =00 -fi-a‘i moom> a O cpl-VI dfli -| I N N | +- I l I a: DJ M N 18 a V‘ajr..f fl‘uu‘vlgmu.lm \‘ \cos 0‘ 1Tcos"a \ : 3(Tr-cos"a) Bza-wh-a—‘ai B=c;5(n(an- cof°'q)) =cisla1‘r-cos"¢) 3(‘T'CO5 a) 0‘0 c=o D: 0+ "aai D:cis( TTCOquJ ) =cis(cos"oc) 3L(17-405 0:) Ez-l ___ . 11‘ E Cls(3(1T-cos"a)) Figure 4.3. The transformation w2 = v? 1r 2(1r -cos-la) 19 SK RE D \V‘ I: = . v(afr—cos'b) ___- B c.5(airr-cos"a)) i=2) C = 0 D: I = - Trcos"a I =- D C's(alrr-cos"a)) E ' E z“ 5(3(1T1I:0$"d)) Figure 4.4. The transformation w3 = cw2 where c = cis e . -1 1r sIn a 9 = -1 TT-COS a \\\\\\-\\\\\\\ mUflCD> .L_.o_.8 mooCDZP 0.4.2..-._ 21 4.5 p = W -l . 2(w-cos a) We take the branch 0 g_e < Zn. The transform given by equation 4.3 rotates the half plane to the region left of the y axis by multiplying by c, which has a magnitude 1. Its argument 9 must be the difference between E and -l the original angle W COS a 2(w-—cos'la) e = 71’ _ 1f cos-la _ 1_r (w—Zcos-la-Zsin‘la+ gsin-la) 2 (1T - COS-la) 2 1r - cos-1 a = 71 223.1123. = 7r sin-1a 2 1 ' 7r - cos-la 1r - cos- a The transform given by equation 4.4 transforms the y axis to the unit circle. Equations 4.1-4.4 can be summarized as 2 P c z + -a + l 2 4.6 w— z" "a l..- . -l c = cis( v sIn-1a) “Ir-cos a p = F (branch 0 g,e < 2?) 2(r-—cos-la) A diagram of this transformation is shown in Figure 4.6. 23 We will also need the inverse of the transformation given by equation 4.6. l _ 2 w+l 13) ‘A' a (1+(c(w-l)) 1 _v_v_:r_l_§_ (C(W-l)) 1 Let z = x+iy . i w=u+1v=re 9. We need to write x and y in terms of r and e. 1 w+l f5) 2 (1+(c(w-1)) l 4.8 x=Re -a w+l i5_ (C(W-1)) l where 49 w+l=j£cose+ll+rsinai ' w-l (rcose-l)+rsinei r -;-i2r sin 8 r2-2r cos 9+1 Notice that equation 4.9 is a complex expression of the form 13 A+ Bi with A g 0, so its argument is 7r+tan- (A) where Ago, or 1r+tan-l(w—a) when rill. Let l-r2 l 1 w+1 ._ _£n (————) = .F_+_1__ P = p C(w-l) 4.10 A Re(c(w-l)) Re e 1 2 2 2 . 2 '- Re exp[l' {£n[(r2'1) +4r Slnz Q12 p (r -2rcos 9+1) + i[7r+tan.l(2r Si; 9)1 -£.n C}] l-r 24 where we have used equation 4.9 in the last step. We will rewrite the expression in equation 4.10 4 ll (r2--l)2+4r2 sin e = r24-2r cos 94;; (r2-2r cos 9+1)2 r2-2r cos 64-1 Recall that [cl = 1, so . . -1 4.12 Ln c = i arg c = 1” Sin 1a . TT-COS- a Putting equations 4.11 and 4.12 into equation 4.10 yields 2 4.13 A = Re exp[l {% 2n (r2+21r~9°s 9+1) p r -2r cos 9+1 . . -1 + i[7T+tan-l(2r 511129)] -i 11’ SIn -2; }] . l-r w-—cos a When r = l, we cannot use equation 4.13 since wi-l) 2r sin e arg (w—- 1 ) = tan-l( 2 isn't defined. However, it l-r can easily be defined in a consistent manner. When r = 1, equation 4.9 yields w4-1 _ -i2 sin e _ -i sin e w-l 2-2cose-l-cose' Thus, when r = l, él‘ O < 9 < n w4-1 2 4.14 arg (m) = g n < e < 2? When W’= l, we cannot use equation 4.14, but it isn't necessary to use it, since we already know that w = 1 corresponds to x = -a2, y = 0. Equation 4.13 can be rewritten 25 l w4-1 5 )) 4.15 A = RetETaij = R cos w 2 .1. r 4-gr cos 94-1 2p R = ( ) r2-2r cos 9+1 __l_)w+ Tr sin-la] co= glarmw—T— 1 - -1 TT-COS a r . “tan-MW) r 7, 1 l-r + arg(::_ i)= ) 3—27T- r=1.0iifip+ 1] l-cos e 1-cos 9 2 l 2 [(i+cos g)p+2(l_2a2)(l+cos e)p+l] =_1_:.§_ "COS l-cos e 4 _]_. i 2 l+cos 9p_ 1+cosfiZp [(1-cos e) 2a(lucos e) + 1:] Define _l_ _ 1+cosfi 2p 4'23 B - (l-cos e) ' Then we can define 4 24 m, = 22.123 = l-a2 m4+2<1 -2a2)(52+41 4 4 (Bz-ZaB+l)2 l-a2 ————4 (1 + 2 435 ) B - ZaB + 1 Finally, we can use Poisson's formula [22] to write our solution to the Dirichlet problem on the circle as 1 Zn 2 4.25 udcp = f(e) f2” 2 2n 0 o 1-re1(9'“p) Io In the first integrand on the right hand side of equation 6.4 we make the substitution 2 = eyI, dz = ieflpdp = izdp, yielding 6.5 15(9)] 2‘22 — Eli-1f ——§E-.— c (l-rele/2)iz 1 c z--re16 where c is the unit circle traversed counterclockwise. The integrand has a simple pole at z = rele, so the value of the integrand in equation 6.5 is 6.6 22%22 x 2vi = 4Wf(6) 33 Putting this into equation 6.4 gives 6.7 I22 G(®)d¢ = 4Wf(e)-2nf(9) = 2rf(9) 0 We substitute equation 6.7 into equation 6.2. We then take the real part of equation 6.2, yielding 2 6.8 I277 f(cp) z 1'13 dcp 0 1+r -2r cos(9-cp) 27" l-r2 =I [f(cp)-f(e)] 2 J dcp+21rf(9) 0 1+r -2r cos(9-cp) In addition, the integration in equation 4.25 is difficult to perform because f(m), as given by equation 4.24, has a sharp peak when 8 = a m l. The denominator of the integrand is zero when 82-2aB+-1 0. Thus the denominator is close to zero when .1. _ 1n+cos @ 2p _ ~ We have used the definition of 5 given by equation 4.23 l+cos cp 2p 6.10 l-cos T — a cos m - a2p__l a2p+ 1 When a is slightly less than 1, p, as given by equation 4.5, is slightly greater than %. Thus one can see from equation 6.10 that cos 6 is a small negative number. Equation 6.10 has two solutions: We attempt to use the scheme described earlier of adding and subtracting a function which is identical to the original function at the point of difficulty. The function we will choose may seem arbitrary, but as the calculation continues it will become evident why this particular function is suitable. 2 . 6.12 (”2” [1+ 2432 J 2 1" 0 B -2aB+l 1+r -2r cos(9-o) 7r 4293 1--r2 =1 1* 2 2 0 B -2aB+l 1+r -2r cos(9-m) 4a2(B2P+ J_.) l - r2 1+ 2 2 2 dcp (B -2aB+l)(a p+1) 1+r -2rcos(9-cp1) + [27 [1+ 4aB J l-r2 2 2 F B -2aB+l 1+r -2r cos(9-m) 4a2(BZP+ 1) l -r2 1* 2 2 2 d“) .(B -2aB+l)(aP+l) 1+r -2r cos(9-m2) l - r2 Tr . 4a2 (B2p+ l) + 2 I 1+ 2 2p de 1+r -2r cos(9-cpl) 0 (B -2aB+l)(a +1) 2 2 2P + 2 l-r I21r 1+ 24a (B +jl)2 Jdcp 1+r -2r cos(9-€02) 1T (B -2aB+l)(a p+1) We wish to change the variable of integration from co to B. Rewrite equation 6.9 as 2 6.13 coscp= 1--——— B2p+l 35 Taking the differential of equation 6.13 yields 2p-l 6.14 -sin cpde = 129—7 dB. (B p+1) We again use equation 6.13 to eliminate sin m in equation 6.14. . 2 6.15 Sln @ = :_ -cos m ==j;~/{- CL--7i%—-)2 B +1 -.,_2.£_ —BZP+1 The sign in equation 6.15 is plus for 0 g_¢ g_w and minus for w g_m g_27. Substituting equation 6.15 into equation 6.14 yields 13-1 - 29 6.16 -2pB dB - :15 4-l)@p. We now see the reason for the factor (B2p4-1)/(a2p4-l) in equation 6.12. As equation 6.16 shows, we need a factor of (6213+ 1) when we change variable, and (62P+ 1) /(a2P+ 1) is equal to l at the troublesome point B = a. thice from equation 6.9 that as m 4 0 or 2w, B 4 a, while as m 4 w, y‘4 0. Thus we can write the second to last integral in equation 6.12 as 2 2p 6.17 (‘2 [1+ 24a “3 +92 Jan 0 (6 -2aB+l)(aP-l) _ 8a2 0 pp-l aZP+1 a BZ-ZaB+l dB 1 “”2932. " 4‘" dB . aZP+1 O B2-2aB+1 36 Similarly, 2W 2 2p 6.18 I [1+ 2 4a “3 +1% Jan 1r (B -2aB+1)(a p+ 1) pp‘l -1r+——E—ea2 . _ 2p 2 ‘22 a +1 0B -2aB+l We integrate equation 6.17 (and equation 6.18, which is identical) by residues. Due to the branch 0 g.9 < 2v of the function zp-l, we integrate over the curve shown in Figure 6.1 (see [24] for a similar calculation). Figure 6.1. The curve for integrating equations 6.17 and 6.18 37 Let _ ap'l _ pp'l B ’2aB+l [B-(a+i./{-a2)][B-(a-i./{-a2)] If CO, C, L1' and L2 are the paths shown in Figure 6.1, and if we take rO small enough, R big enough, and 6 small enough that the poles at a4-i./l-a2 and a«-i --a2 are on the inside of the contour, we get 6.20 IL 9(B)dB+j gdB+fL g(B)dB+fCO Shams l C 2 (3+ i ‘A-a2)p-l + (a-i ./{-a2)p-l 2i ./I -a2 -2i \/1 --a2 = Zwi On L1, B = rele, flp-l = rP-lel€(p'l), dB = eIé dr On L2, B = re1(27-6) = r -i€' Bp-l = rp-lei(p_1)(27_€)' e Thus rp-l ei(P-l)€ eie . . x r2e21€-2arele+-l R 6.21 (L1 g(B)dB+J‘L g(B)dB = fr dr 2 O r rp-l eI(p-l)(27r-€) -i6 0 + ~¥ o c e dr IR rZe-ZIE IE -2 are- +1 The integrands are continuous functions of r and E in the closed region rO g.r g_R, 0 g_€ g_g, so 38 6.22 lim [ g(B)dBi' g(B)dB] 6+0 J«1‘1 IL2 p-l 2 r dr r -2ar+l . R _ 21(p-l)w - (l-e ) fro The integral of 9 over C, where B = Re19 can be I written Rp-l ei (P-l) 9 2w—E . ie 6.23 g(B)dB = . - Rle d9 ‘fc Io R2e219 - 2ARe19 + 1 1' j (B)dB 'rP f2" eipe de 1m 9 .= I . . 6+0 0 0 nzezle - 2aRele + 1 p Ilim‘f g(B)dB] g 2 24” 6+0 0 R - 2aR - 1 The expression on the right hand side of equation 6.23 approaches 0 as R approaches m, since 0 < p < 1. Similarly, . . 2W eipe 6.24 11m g(B)dB = -Irp . . d9 6+0 J2C0 0 J0 rge219-2aroele+-l P rO 2r 1-2ar —r2 Ilim J‘C g(B)dB) g 6+0 0 O O The right hand side of equation 6.24 approaches 0 as rO approaches 0. Thus we conclude from equation 6.20 and 6.22 that o flp-l 0 B2-2aB+1 = ____1r___ [(a+ 3'. MI -a2)P-l - (a-i ./l-a2)p-1] . 6.25 (1-e2i(P’1)"T) J“ as 39 Taking the real part of equation 6.25, O Bp-l 0 B2 - 2aB + 1 ——-TT—— Re[(a+ i ./1-a2)P—l - (a — i MI -32)p-l] -az 6.26 (l—cos 21r(p-1))j‘ dB To simplify equation 6.26, 6.27 Re(a+ i .4 ~12)!"1 (p-l)£n(a+i -a2) Ree (p-l)i cos.1 a =Ree l = cos[(p- l)cos_ a]. Similarly, 2 P_l (p-1)i(2w-cos‘1a) 6.28 Re(a-i -a) = Re e 1 = cos(p - l) [21T ‘cos- a] = cos 27r(p-l)cos[(p-l)cos'l a] + sin 27T(p- l)sin[(p- 1) cos-1 a]. Then equation 6.26 becomes 6.29 I. 2 BP-l dB = -—L— [cos[(P-l)cos-l a]} 0 6 -2aB+ 1 \A_a2 - cos 27r(p-l)cos[(p-l)cos-1 a] - sin 27r(p—1)sin[(p-l)cos-1 a] /(l-cos 21T(p-l)) = —72—- {c08[(p-l)cos'l a1} .a2 - cot 1T(p-l) x sin[(p-l)cos_l a] . 40 With the help of equations 6.17, 6.18, and 6.29, equation 6.12 can now be numerically integrated on a computer. CHAPTER VII RESUDTS The numerical integration was carried out by the program SIMP [25]. This program is an adaptive procedure based on Simpson's rule. It was altered to do the double integration indicated by equation 5.1. An adaptive pro- cedure uses a finer mesh in regions where the function is changing rapidly, taking smaller and smaller subintervals until sufficient accuracy is attained, which was, in our case, .05%. There are internal limits in the program to prevent the integration from continuing forever if the function is extremely ill-behaved; the program will not take more than 2000 function evaluations or a subinterval smaller than 2"30 times the original interval length. When the integration of equation 4.25 was originally performed, the program reached the limit of 2000 function evaluations. This happened even when the limit was increased to 10,000 evaluations. Only When the techniques discussed in the beginning of Chapter VI were introduced did the function become smooth enough for SIMP to integrate. The techniques discussed in the latter part of Chapter VI were not needed for the present study. Even with a = .98, 41 42 it took more computation time with these techniques than ‘without them (130 seconds vs. 78 seconds). If it had become necessary to investigate values of a closer to l, the techniques would eventually need to be utilized. The results of the numerical integration of equation 5.1, are plotted in Figure 7.1, where the flow rate Q is plotted against the fill level a. As expected, the flow rate increases for small values of a. However, 0 becomes larger than 1 (that is, the flow is greater than that of a full cylinder) for .38 < a < 1. Q reaches a maximum of 1.273 when a = .724. This means that the flow can be increased by as much as 27.3% when the cylinder is partially full! For a < -.48, we do not have data for the flow rate because of difficulties in numerically integrating equation 4.25; the program continually reached the upper limit of 2000 function evaluations. These difficulties have not yet been resolved, but we do have some ideas about their cause. One can see from equation 4.12 that as a approaches -1, the argument of C approaches n. Thus care must be taken in choosing the pr0per branch of the function zl/P. Fortu- nately, this difficulty does not affect our results, since it is clear from Figure 7.1 what is happening when a < -.48. We have included in Figure 7.1 the graph of the approximate value of Q for values of a close to -l [5]. The velocity profile for the near optimum case of a = .707 is shown in Figure 7.2 (this case is of particular interest and is examined in the next chapter). For the full 43 .320 op‘l’imum value Lam)- 0= .73” Q= [.373 |.l., LO" Jfiu 05‘) .7" ‘6... .5+ .49 OPPrOximdl'e Q {or shallow new [5] .3v .3" J2. J ; : :2 : . i ‘ t' I?’<, -l.o '18 -.6 :9 -.a o .3 H .6 .5 1.0 Figure 7.1. Flow (Q) vs. fill level (a) 44 Figure 7.2. Velocity profile for a = .707 45 cylinder case, one can see from equation 3.2 that the maximum velocity, achieved at the center, is v = .25. In Figure 7.2 we observe a very large region where the velocity is greater than .25. This helps explain why such a significant improvement is achieved by not filling the cylinder. For purposes of comparison, the velocity profile for the full cylinder is shown in Figure 7.3. 46 Figure 7.3. Velocity profile for a full cylinder CHAPTER VIII 1 THE PARTIALLY FILLED CYLINDER WITH a = —- In order to verify our numerical results, we examined a related problem. Suppose 02 = 45° in Figure 4.1. Then the circles in Figure 4.1 intersect at right angles on the exterior, and a = cos 02 =-JL== .707107. As we mentioned in Chapter III, Sokolnikoff and Sokolnikoff [18] solved the related torsional pro- blem in which the interior angle is a right angle (or 02 = 135° in Figure 4.1). We will solve our problem using their method, though the computations are much more difficult in our problem. Suppose we have a mapping 2 = h(w) which takes the unit circle in the w plane to a region in the 2 plane, where z = x4-iy. According to [18], a solution to the problem vzw = 0, with w = %(x?-+y2) on the boundary of the region in the 2 plane, is given by 8.1 w(x.y) = Re {3%} I 1-1--(-%-2:-2559-2-do-I-constant] V where y is the unit circle traversed counterclockwise, 0 is a point on the unit circle, and w is the point in 47 48 the w plane corresponding to 2. Once we find t, we obtain the corresponding solution to the fluid flow problem as [%(X24'y2)-W] = 0 on the boundary <: u haw 8.2 l. 4 < II The constant in equation 8.1 may be determined by setting v = 0 at a point on the boundary. We can use the transformation already developed in Chapter IV. When a = l/Vr, equation 4.6 yields simple values for c and p: cis - cis 60° II o P. m (ma 8.3 c 6:) u l ODIN 2(1r-g) Substituting these values of a and p into equation 4.7 yields £(1+(___ w+l )3/2) 8 4 z = h(w)= C(w 1) ' (w+1 )WQ-l C(W’- 1) We then construct 1 0+1 3/2][ 8+1 3/2] l+(---1--) 1+(:-_—-—) 5 [ C(U'-1) C(O- 1) 0+1 3/2__ 6+1 3/2_ °‘°'1)) 1] EEG-1)) 1] Equation 8.1 is integrated by making the change of variable +.2 8.6 C=_l._.%. l-it 8.5 h(o)fi(o) = [ ( 49 When t traverses the imaginary axis from +mi to 0 and then traverses the real axis from 0 to +m, t2 traverses the real axis from -m to m, and 0 traverses the unit circle counterclockwise (see Figure 8.1). Thus 2 2 0+l=l+it +1-it = 1 °‘1 1+it2-(1-it2) cis 90° t2 . O 8.7 {-c-‘(lail-l—fi)3/2=( 1 2)3/2 = (C18 310 )3/2 cis 60° xcis 90° t t = cis 315° = -cis 135° t3 t3 When t goes from 0 to m (i.e., t is real), the conjugate of equation 8.7 is 6+1 )3/2 = cis 45° 8.8 ( . 6(6-1) t3 When we integrate from mi to 0, let t = is = cis 90°s, t2 = cis 180°s = -sz. Then . 2 8.9 0’ =l’_-Ls.2_ 14-18 and .1L1__3/2 _ 1 3/2 cis 60° x cis 90° x cis 180°s (cis 30°)3/2 = cis 45° 52 33 The conjugate of equation 8.10 6+1 )3/2 = -cis 135° 5(5-l) s3 8.11 ( I: ==+co 1+it2 I-it2 The transformation 0 = Figure 8.1. 51 Also, from equations 8.6 and 8.9, 8.12 do = 42t2 2 dt (l-it) = -413 ds . (l+isz)2 Substituting equations 8.7, 8.8, 8.10, 8.11, and 8.12 into equation 8.5 and substituting the result into equation 8.1 gives 8.13 .1.J~ 11121182101, 27ri 0 - V w ' O ' O ]_2_.(l+c15345 )(1 _CIs 1335) I0 s s (-4i§)ds _ 27ri cis 45° —cis 135° _ . 2 °°( 3 '1” s3 ‘1) (l+is2)2(-l——l—s-2--w) S l+-is 1(l_ci5135°)(1+cis45°) + 1 I“ 2 t3 t3 (4it)dt 21ri -ci5135° cis 45°_ . 2 l-it = 1 I“ (33+c_i§45°)(s3-ci3135°) 8 ds v 0 (s3-cis45°)(s3+cis 135°) l+is2 [l-isz—w(l+isz)] 1f” (t3+cis 45°)(t3-cis 135°) t dt _ + 1.1) 7’ 0 (t3-cis45°)(t3+cis135°) l-it2 [1+it2-w(l-it2)1 We will integrate by parts, so we need to factor the integrand. SBA-cis 45° = (s-+cis 15°)(s-cis 75°)(s4-cis 135°) s3-cis 45° = (s-cis 15°)(s+-cis 75°)(s-—cis 135°) $3-+cis 135° 53 - cis 135° (s4—cis 45°)(s-—cis 105°)(s+-cis 165°) (s-cis 45°)(si-cis lOS°)(s-—cis 165°) 2 (l+isz)(l-is -w(1+is2)) = (s -cis 45°)(s-+cis 45°)(sz(l+ww)+-cis 90°(l-w)) 52 (1— it2)(l+ 11:2 -w(l-it2)) = (t-cis l35°)(t+cis 135°)(t2(l+w)-cis 90°(l-W)) We can then rewrite equation 8.13 as 8.14 LI MM (10' 21m. o-w =21. (s+cis 15°)(s-cis 75°)(s+cis 135°)(s-cis 45°) 0 (s+cis 105°)(s-cis 165°)s /[(3-cis 15°)(s+cis 75°)(s-cis l35°)(s+cis 45°) (s-cis 105°)(s+cis 165°)(s-cis 45°)(s+cis 45°) (s2(1+w) +cis 90°(l-w))]ds + %J' (t+cis 15°)(t-cis 75°)(t+cis 135°)(t-cis 45°) 0 (t+cis 105°)(t-cis 165°)t /[(t-cis 15°)(t+cis 75°)(t-cis l35°)(t+cis 45°) (t-cis 105°)(t-cis 165°)(t-cis l35°)(t+cis 135°) (t2(l+w) -cis 90°(1-w))]dt = 12512 (s+cis 15°)(s-cis 75°)(s+cis 105°)(s+cis 135°) 0 (s-cis 165°)s /[(s-—cis 15°)(s+cis 45°)2(s+cis 75°)(s-cis 105°) (s-cis 135°)(s+cis 165°)(32(l+w)+cis 90°(l-w))]ds +3151. (t+cis 15°)(t-cis 45°)(t-cis 75°)(t+cis 105°) 0 (t-cis 165°)t /[(t-cis 15°)(t+cis 45°)(t+cis 75°)(t-cis 105°) (t-cis 135°)2(t+cis 165°)(t2(l+w) -cis 90°(l—w))]dt 53 To evaluate the first integral in equation 8.14, set 8.15 (Si-cis 15°)(s-cis 75°)(s-+cis 105°)(si-cis 135°) (s-cis 165°)s , /(s-cis 15°)(s+-cis 45°)2(s4-cis 75°)(s-cis 105°) (s-cis 135°)(s+cis 165°)(sz(l+w)+cis 90°(l—w)) A B C D = + + + . (s+cis 45°)2 s+CIs 75° s-cis 15° s-+cis 45° E + F + G s-cis 105° s-cis 135° s-tcis 165° Hs+-I 1' 2 s (1+w)+cis 90°(l-w) As we shall see later, certain values of w will require special treatment. Let us denote f(s) = (sd-cis 15°)(s-cis 75°)(s-+cis 105°)(s+-cis 135°) (s-cis 165°)s fl(s) = s-cis 15° f2(s) = (si-cis 45°)2 f3(s) = s+-cis 75° f4(s) s -cis 105° f5(s) = s-cis 135° f6(s) = s4-cis 165° sz(l+w) +cis 90°(l-w) f7(s) We can evaluate most of the constants by multiplying equation 8.15 by the denominator of the term with the desired constant and evaluating at the root of that denominator. Thus 54 8 16 A = f2(S)f3(S)f:E:;f5(s)f6(S)f7(S)!s = Cis 15° C = f1(s)f3(s)f:(:)£5(s)£6(s)f7(s) s = ‘Cis 45° D = fl(s)f2(s)ff%:;f5(s)f6(S)f7(s) s = -Cis 750 E = f10, y710 4x -2+4y §21 if 4x2-2+4y2=0, y740 L In the first case of equation 8.25, the argument is in quadrant III but tan-l( E4‘VQ y 2) is in quadrant I, 4x -2+4y which is why we add W. In the second case, tan-l( -4‘ng ) 2 2 Mc-2+4y is in quadrant IV, but we need to add 2w because we have chosen the branch 0‘g_9 < Zn. In the second and third cases we have assumed y # 0. In the second case, if y = 0, the argument should be zero rather than Zn. In the third case, if y = 0, then x = ipvfi/Q, in which case ‘w = :1. Denote 2 2 . 8.26 6 = arg(4x -2+4y -4\/2 yI) Then equation 8.24 becomes 1 8.27 c+-IB = l giznl(4X2-2+4y2)2+32y2]2+ié) e (4x2-4\/2 x+2+4y2)2/3 - ((4x2--2+4yz)2+32y2)1/3 . ) ‘ 2 2 2/3 C13‘3 5 (4x -4.¢§ x4-24-4y ) 60 Then equation 8.25 becomes -511 8.28 w - §-1 where g = (cis §)(a+-iB). There is one other difficulty to take care of. The partial fraction expansions in equations 8.15 and 8.20 are not valid for values of ‘w Which cause one of the roots of 32(1+w)+cis 90°(l-w) or t2(l+w)+cis 90°(w-1) to be equal to the root of any other factor. The roots of these two expressions are _ = - o 2L2}; 8.29 5+, 5 i915 45 w4-l _ . l-w If we compare these solutions with the roots of the other factors, we find that the only value of W’ in the unit circle which causes any difficulty is w = 0, that is, at the center of the unit circle. In that case 8.30 32(1+w)+cis 90°(l-w) = sz+cis 90° = sZ-cis 270° = (s-cis 135°)(si-cis 135°) and 8.31 t2(l+w)+cis 90°(l-w) = tz-cis 90° = (t-cis 45°)(t4-cis 45°) Putting equations 8.30 and 8.31 into equation 8.14 and simplifying yields 8.32 Here gral 61 % I; (Si-cis 15°)(s-cis 75°)(s-+cis 105°)(s-cis 165°)s /(s-cis 15°)(s-+cis 45°)2(s4-cis 75°)(s-cis 105°) (s-cis 135°)2(s+cis 165°) ds + % I; (t4-cis 15°)(t-—cis 75°)(t-cis 105°)(t-cis 165°)t, /(t-cis 15°)(t4-45°)2(t+-cis 75°)(t-—cis 105°) (t-cis 135°)2(t+cis 165°) dt - 31" A + B + C + D - w 0 s-cis 15° s4-cis 45° (Si-cis 45°)2 s-+cis 75° E + F s-cis 105° s-cis 135° G H + (s-cis 135°)2 + s4-cis 165° d8 =12; [A zn(s-cis 15°)+B £n(s+cis 45°) -c /(s+-cis 45°)4-D Ln(s4-cis 75°)4—E Ln(s-—cis 105°) + F £n(s-cis 135°) -G/(s-cis 135°) + H £n(s4-cis 165°)] . . O we have used the fact that, in this case, the s inte- and t integral are identical. To evaluate the constants, define d(s) = (s4-cis lS°)(s-—cis 75°)(s4-cis 105°)(s-cis 165°)s s - cis 15° d1(8) d2(s) = (s-+cis 45°)2 d3(s) = s4-cis 75° d4(s) = s-cis 105° d5(s) (s -cis 135°)2 d6(s) = s4-cis 165° 62 We conclude as before d(s) _ |_- o 8'33 A - d2(s)d3(s)d4(s)d5(s)d6(s)Is - Cls 15 _ d(S) I _ _ ' o C ‘ dl(s)d3(s)d4(s)d5(s)d6(s)|s ‘ C15 45 _ d(S) _ ' o D ‘ d 14(s)dz(s)d (s)d5 (s)d6 (s)ls ‘ '°13 75 _ d(S) I _ . o E ‘ d1(s)d2(s)d3(s)d5(s)d6(s)ls ‘ C13 105 G d(s) = cis 135° - I _ dl(s)d2(s)d3(s)d4(s)d6(s)ls _ d(Sng - ° 0 H ‘ d 13(s)dz(s)d (s)d4 (s)d5 (s)|s ‘ 'Cls 165 To find B and F, we define _ d(s) 8.34 X(S) - d1(s)d2(s)d3(s)d4(s)d5(s)d6(S) B + F s-+cis 45° s-cis 135° The last equality in equation 8.34 follows from the first equality in equation 8.32 B F 8'35 x(0) = cis 45° - cis 135° B F X(2 Cls 135°) = 2 cis 135°4-cis 45° + cis 135° Adding the above equations yields ' o 8.36 B = IX(O) +X(2 C18 %35 ) cis 45° + 2 cis l35°4-cis 45° 63 Equation 8.35 then gives us B . 8.37 F = (EYE—45? - X(O))CIs 135° At last we have everything needed to calculate v(x,y) for any point in the upper half of the region in Figure 4.1. The velocity was integrated numerically over this region using SIMP, as described in Chapter VII. The result was that Q = 1.273, which is exactly what the methods of Chapters IV through VI predict. The velocity distribution for this case was plotted in Figure 7.2. CHAPTER IX PERTURBATION SOLUTION OF THE ALMOST HALF FILLED CYLINDER In this chapter we shall use perturbation techniques to solve the problem of flow through a cylinder filled just under or just over halfway. We were unable to use such a technique for the almost full cylinder. This is because the full cylinder is quantitatively different from the partially filled cylinder, since the latter has a free surface and the former doesn't. Consider a circle of radius ./1+-€2, where 6 is small compared with l, and center at (0,6), as shown in Figure 9.1 (though 6 > 0 in Figure 9.1, all of our analysis is valid for E < 0, that is, for a more than half filled cylinder). The equation of this circle is 9.1 x2+y2-2y€-1=0. Writing equation 9.1 in polar coordinates, we get 9.2 r2-26r sin 9-l = 0 64 Figure 9.1. Flow through the almost full circular cylinder . 2 . 2 . 9.3 r = 6 sIn 9 +./{4-6 sIn 9 (SInce rug 0) 2 = 14-6 sin 9 +'%r sin2 9 + 0(64) USing the same argument as in Chapter IV, we can consider this problem to be half of the flow between two intersecting circles, as shown in Figure 9.2. Equation 9.3 gives the equation for the bottom half of Figure 9.2, i.e., When sin 9 g_0. The equation for the t0p half, i.e., when sin 9,; 0, is 9.4 x2+y2+2y6-l=0. 66 l . '- I rel-€51n6+-§Srn6 r. 6 , A a . a Aresmew‘gsm 9 Figure 9.2. Cross section for the almost half filled cylinder In polar coordinates, r24-26r sin 9-1 = 0 2 9.5 r = l-6 sin 9 +‘%r sin 2 9+0(64) Thus we can state our problem.(with accuracy to C(63)) as 9.6 V2v = 1 2 v = 0 When r = l-6|sin 9) +‘%? sin 9 ExPanding v(r.9) in a Taylor series about r = l, 2 9.7 v(l-El sin 8] + 32— sin2 e) = v(1,e) 2 - .€_ - 2 927. + (-E[51n.el +-‘2 Sin 9)ar (1.9) 2 2 + 2(‘elSin 9! +6— sin2 e)2 a v(l,9) +... 2 arZ 67 For general r, we can also expand _ 2 3 9.8 v(r.e) - v0(r.e)+6vl(r.e)+€ v2(r.e)+0(€) Putting equation 9.8 into equation 9.7 yields 62 2 9.9 v(l-—€lsin 9| +‘7f sin e) _ 2 — vo(l,e)+-Evl(l,e)4-€ v2(1,e)+... E2 av0(1.e) av1(l.e) + ) . __ . 2 + (-€Isa.n 91+ 2 81!) 6)( 3r + 6 Br 2 2 a v (1.9) + %(-€I sin 91 +-e-2- sin2 e)2( O 2 +...) +0(€3) ar . av0(1.e) = vo(1.e) + E vl(l.e) -lsa.n 9! 3r . 2 3V (lye) 3V (1.9) 2 O . + e [v2(1.e)+="-“2‘—a 3r ‘lsm el—l'ar_' 2 + Slgfi O 2 ].+ 0(63) ar We now compare like powers of E in equations 9.6 and 9.9. First we have the terms of order 60. 2 _ 9.10 v vO - 1 In polar coordinates, 2 2 9.11 v2=L§+%_a§;+_lz._§_2_ 5r r as The system in equation 9.10 has the solution 2 _ r --1 Next, we compare terms of order 61. 2 _ 9.13 v vl - 0 av . o _ VIN-.9) “ISln 9| ar(loe) - O 68 which implies 9.14 v1(1.e) = Js_1r21_iL The solution to equation 9.13 is a harmonic function; we seek a solution of the form a 9.15 vl(r.e) = -Q'+ Z) akrk cos k6 2 k=1 We have not included the sin ke terms in the expansion since our solution must be an even function of e. We will solve for the constants ak by expanding the boundary condition in equation 9.14 as a Fourier series. a0 an 9.16 [sin e] =7+ 23 a cos k9 k k=1 w w ak = i I cos kxlsin xldx = Z I cos kx sin xdx w -w w 0 k = 0.1.2,... If k = l, 2 V . a1 = % I cos x Sin xdx = O 0 Otherwise, 1 MT . . 9.17 ak = E J [Sin(k4-1)x-Sin(k-l)x]dx 0 _ 1 {-cost4-ll§.+ cos§k-—lzx] W — v ki-l k-l 0 0 if k is odd #lk-IZ-l - TIE—l if k is even, say, k = 21. 69 Thus [sin 31 _ 1 _2_ 2 °° 1 9°18 2 ‘ 2 [1r + 1r 1.21 (n+1 21- 11) °°s “91 _ 1 _ 2 ° 1 " 77 3; 1:51 (21+1)(2z-1) C°s “9' Thus co 2!, -1-2 r The terms of order 62 are 9.20 vzvz = o 2 avo (1.6) av (1.6) sin 6 - __l___—— sinze 62 V0 (109) + = O :2 2 Br Substituting equations 9.12 and 9.19 into equation 9.20 gives 9.21 v2(1.e) = :sigii x % ~ sigifi % ‘ i #2; T2k4-l§?ZR-l) °°s ZRBISin 9' = _(l-czs 28) ‘ $121 [(2k+1)](<2k-1)(a—2q+ “:31 am °°S m9” =;4.L+C_°§_£Q-14.T§3[k l(:-2()+Z)1amcosme)] k=1 (2k)2 m=1 where we have expanded cos 2kelsin 9] as a Fourier series. 9.22 70 _ 1 7T . am - fi f—w [cos me cos 2keISin 9]]de 2 Tr . = - I [cos me cos Zke s1n e]de 7’ o 1 1T . . = 1? I [cos(m-2k)6 Sin e+cos(m+2k)e Sin e]de 0 _ .1. 7T - 2k 1 ' 2k 1 - 2” f0 [s:.n(m- + )e-s:|.n(m- - )9 + sin(m+ 2k+ 1):: - sin(m+ 2k- 1)9]de If m-2k is odd, then m-2k+l, m-Zk-l. m+2k+l, and m- 2k- 1 are all even and am = 0. Thus we have am 7! 0 only when m-Zk is even, i.e., when m is even. Letm 9.23 Thus 9.24 = 2n. 1 1T . . a2n -— 2? IO [s:1n(2n-2k+1)e-Sin(2n—2k-l)e + sin(2n+2k+l)e-sin(2n+2k-l)6]de _ ; {-cos(2n—2k+1)a __ cos(2n+;k+l)e - 21r 2n-2k+1 2n+2k+1 + cosfln-gk-lze + cosL2n+2k-IZQ]ITT 2n-2k-l 2n+2k—l O _ 1 I 1 + 1 _ 1 _ 1 1 -11" 2n-2k+1 2n+2k+l 2n-2k-1 2n+2k-1 =12 (2n)2+L2k)2-1 7’ [(2n+1)2-(2k)2][(2n-1)2-(2k)2] _ _ 1 cos 29. "$7 ‘ +[ 122 + 23 k=l (2k) -1 TTHZk) -1 n=1 _ g [(ZnL2 + (2k)2 - 1]cos 2ne J 1T [(2n+ 1)2 - (21021 [(Zn- 1)2 - (21021 1 cos 29 8 °° k = _ _ + + E 4 4 n2 k=l ((2):)2-1)2 + $22 2 [Z k((2n)2+(2k)2-1) w n=l k=l /< (2k)2 - 1)((2n+ 1)2 - (2102) ((2n- nz - (2k)2)]cos 2ne where we have applied Fubini's theorem to justify switching the order of summation. Since we are looking for a solution of the form of equation 9.15, we conclude from equation 9.24 that 1 cos 2g 2 b0 m Zn 9.25 v (r,n) = - - + r + ——-+ Z) b cos 2ne r 2 4 4 2 n=1 2n b =& z k((2n)2+(2k)2-l) 2n #2 =1 /((2k)2 - 1)((2k)2 - (2n+ 1)2)(<2k)2 - (2n- 1>2) To get the amount of flow, we integrate the velocity over the upper half of the region in Figure 9.2. . 2 sinze 1r 1-6s1n 6+6 -—2—- I I v(r,e)rdrde O O 9.26 Q 2 3 Q0 + 601 + 6 Q2 + 0(6 ). If F(r,e) is an antiderivative of v(r.e)) 4 We have divided by -w/8 as in equation 3.3. Hence- forth we will drop the primes. We can now construct a new figure corresponding to Figure 9.1 but with radius 1 (see Figure 9.3). Figure 9.3. Flow through the almost half filled unit cylinder 75 In order to compare this result with the results of Chapters IV through VII, let us write equation 9.34 in terms of the parameter a. -6 9.35 E = Substituting equation 9.35 into 9.34 yields 2 a 5 .3. ...; n! =1”) l-a Equation 9.36 is plotted in Figure 9.4 for -.2 g.a g .2. 0n the same graph we have plotted Q as found by the methods of Chapters IV through VII. we see that the agree- ment is excellent for small values of a. 76 Q as given by Ike meH-ods of chap'l'ers I! 'I’ItrougI! II Q as given by equ‘I‘ion cI336 A L J 1 Ag; A A A A A A A A L A L L L V r v t v v r r r v two '3‘ 412' {sf-.3» '-.30 'J‘ ma not '54 I d ...| Figure 9.4. Plotting the flow using perturbation methods CHAPTER X DISCUSSION Chapters IV through VII give the solution to the problem of flow through a partially filled cylinder when a 2 -.48. As we mentioned in Chapter VII, there are unresolved difficulties of integration when a < —.48. We h0pe to resolve these difficulties in the near future. Relaxing the assumptions set forth in Chapter II opens up many new and interesting areas of research. Since the underlying theory of turbulent flow is not well understood, it would be very difficult to relax our assumption that the flow is laminar (though this problem has been solved empirically, as is mentioned in Chapter I). We could, however, study what might happen if the free surface is not flat. Often waves develop in the surface of fluid flowing in an open channel. These waves may affect the flow rate, and they would also complicate the problem considerably (see [26]). The shape of the pipe could also be changed. For example, if the pipe had an elliptical rather than a circular cross section, we could solve the problem using 77 78 elliptical coordinates. We could also consider the case where the pipe had a slight curvature. Dean [27] studied the problem of flow through a fully filled toroidal cylinder. The computations are quite complicated; the level of com- plexity would increase even more if the cylinder was not completely filled. A very natural follow-up study would be to examine experimentally the results reported in this thesis. It should be possible to set up an experiment and determine the flow really does increase when the cylinder is partially filled. It is our intent to perform this experiment some time in the near future. We have a great confidence that the experimental results will confirm What we have demonstrated mathematically. APPENDIX ANGLE OF INCLINATION FOR LAMINAR FLOW We wish to determine the maximum angle of inclination that will cause flow through a circular cylinder to be laminar. For simplicity we will consider flow through a full cylinder. From Chapter II, F3 = 9 sin a, where d is the angle of inclination and g = 981 cmz/sec. USing the normalization from Chapter II and the solution to the full cylinder from equation 3.2, 2 A.l V’=;v1§=r4-1. FZL Thus A 2 v = [l-—r2)L2g sin 0 ° 4v Since we need the average velocity for the Reynolds number, A.3 v £ I I (l r IE q Sln c r dr d6 7’ o o avg 4 _ ngsin a - 8v ° Using equation 2.1 with L = 2L, . 2L v 3 . R=_m=Lg:1na< 2320 V 4v 2 21.4 a < sin'1(93-8—§V—) . L 9 79 80 For example, if crude oil, with a typical viscosity of v = .5 cmz/Sec, were to flow through a pipe with radius 25 cm, equation A.4 tells us that a must be less than .009° to guarantee laminar flow. If this same oil ‘were to flow through a pipe of radius 1.33 cm, it would be laminar for any angle of inclination. BIBLIOGRAPHY 10. BIBLIOGRAPHY Von Mises, R. and Friedrichs, K.O., Fluid Dynamics, Springer Verlag, (1971), pp.194-205. Camp, T.R., "Design of Sewers to Facilitate Flow," Sewage works Journal, V61. 18, (1946), pp.1-16. Chow, Ven Te, Qpen Channel Hydraulics, McGraw-Hill, (1959). pp.134-136. Chaudhury, T.K., "On the Steady Flow of Viscous Liquid under Constant Pressure Gradient in a Cylinder of Lenticular Section," Revue Roumaine des Sciences Techniques Série de Mecanique Appliquee, Vol. 9. (1964). pp. 759-766. Buffham, B.A., "Laminar Flow in Open Circular Channels and Symmetrical Lenticular Tubes," Transactions of the Institution o§;Chemical Engineers, VOl. 26. (1968). PP.T152-T157. Charles, M.E., and Redberger, P.J., "The Reduction of Pressure Gradients in Oil Pipelines by the Addition of Water: Numerical Analysis of Strati- fied Flow," The Canadian JOurnal of Chemical Engineering, Vol. 40, (1962). PP.70-75. Gemmell, Alan R., and Epstein, NOrman, "Numerical Analysis of Stratified Laminar Flow of Two Immisible Newtonian Liquids in a Circular Pipe," The Canadian JOurnal of Chemical Engineering, VOl. 40, (1962). pp.215-224. Yu, H.S., and Sparrow, E.M., "Stratified Laminar Flow in Ducts of Arbitrary Shape," A.I.Ch.E. JOurnal, V010 13' (1967), pp. 10‘160 Goldstein, 8., Modern Developments in Fluid Mechanics, Vol. 1, (1965). p.297. Synge, J. L., "Hydrodynamic Stability," Semicentennial Publications of the American Mathematical Society, Vol.2 (Addresses), (1938), pp. 227- 269. 81 ll. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 82 Pai, Shih-J., Viscous Flow Theory, Vol. 2 - Turbulent Flow, (1957), p.3. Goldstein, 8., 0p. cit., p.304. Batchelor, G.K.,.An Introduction to Flpid Dynamics, Cambridge university Press, (1970), pp.66-67. 596-597. Shinbrot, Marvin, Lectures on Fluid Mechanics, Gordon and Breach, 1973, p.108. Aris, Rutherford, Vectors, Tensorsyiand the Basic Equations of Fluid Mechanics, Prentice-Hall, (1962), pp.109-110. 231-232. Love, A.E.H., A Treatise on the Mathematical Theory of Elasticity, Dover, (1944), p.313. Timoshenko, S., and Goodier, J.N., Theoryiof Elasticity, McGraw-Hill, (1970), pp.303-306. Sokolnikoff, I.S., and Sokolnikoff, E.S., "Torsion of Regions Bounded by Circular Arcs," Bulletin of American MathematicallSociety, vol. 44, (1938), pp.384-387. Wigglesworth, L.A., and Stevenson, A.C., "Flexure and Torsion of Cylinders with Cross-Sections Bounded by Orthogonal Circular Arcs," Proceedlggs ‘2; the Royal Soclety of London, A-l70, (1939), pp.391-4l4. Lal, Krishna, "On the Steady Laminar Flow of a Viscous Incompressible Fluid Through a Tube whose Cross Section is Formed by the Inter-section of Two Circles," JOurnal of Technology, vol. 11, (1966), pp.21-25. Navanlinna, Rolf, and Poatero, V., Introduction to Complex Analysis, Addison-Wesley, (1969), pp.191-196. John, 1., Partial Differential_§quatlons, Springer- Verlag, (1975), p.130. James Lyness, private communication. Churchill, R.V., Complex Variables and Applications, McGraw-Hill, (1960). PP.168-l71. Shampine, Lawrence, 1., and Allen, Ridhard C., Numerical Computing: An Introduction, W.B. Saunders Company, (1973). pp.238-242. 26. 27. 83 Duncan, W.J., Thom, A.S., and Young, A.D., Medhanics of Fluids, American Elsevier Publishing Company, (1970). PP.424-458. Dean, W.R., "The Stream-line Motion of Fluid in a Curved Pipe," PhilosophicalfMagazine, Vol. 5, (1928),