THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING THE FLOW DUE TO A ROTATING DISC WITH A CENTER SINK Charles Go Richards A dissertation submitted in partial fulfillment of the requirements for'the degree of Doctor of Philosophy in the University of Michigan Department of Engineering Mechanics 1964 December, 1964 IP 688

Doctoral Committee: Associate Professor William P. Grae'bel, Chairman Professor Samuel K. Clark Associate Professor Walter R. Debler Professor Arthur G. Hansen Assistant Professor James 0. Wilkes

ACKNOWLEDGEMENT The author wishes to express his gratitude for the assistance and advice of the members of his Doctoral Thesis Committee, A special expression of gratitude is due to Associate Professor William P. Graebel, Chairman of the Committee. The author is also grateful for the technical assistance given him by Mr. William Go. Huizenga and Mro, Milo W. Kaufman. ii

TABLE OF CONTENTS Page ACKNOWLTEGEMENT......... *4*4..4 ii LIST OF TABLES.... a a........e0... a 4 a a O V LIST OF FIGURES.......,... 4 o o. *... a'a vi CHAPTER 1 INTRODUCTION..4*........a.o........Oa *.*. 0*.. 1 2 THE DIFFERENTIAL EQUATIONS GOVERNING THE FLOW DUE TO A ROTATING DISC WITH A CENTER SINK,.0......4 o 5 2,1 Statement of the Problem........................ 5 2,2 Dimensionless Form,..... 4................... 12 3 THE FINITE DIFFERENCE FORM OF THE EQUATIONS GOVERNING THE PROBLEM OF THE ROTATING DISC WITH A CENTER SINK.. a 0 a a. * a o 4 a o * a a * 44 a 4 *. a a4 17 3o1 Introduction,...................... 4 * 4. 4 4 17 3,2 System of Grid Points,......................... 19 3,3 Finite Difference Approximations to the Equations of Motion..,...,. 4....... * * * 4 4 * 4 20 3.4 Finite Difference Approximations to the Boundary Conditions........................................ 26 355 Stability of the Method.....,.......... 29 4 NUMERICAL RESULTS,..ao.....oo o....... 531 4,1 The Dividing Streamline,..,.............e a 5...o 31 4,.2 Preservation of Circulation..4..,,... * *.....* 33 4,3 Influence of the Sink at Large Zoo....a.o.oao, 34 4,4 Influence of the Parameter,..0.,...... 38 4,5 Technique of Solution..... 4 a........... o..... 46 4,6 The Rotating Disc with a Center Source,,,,,...... 50 5 EXPERIMENTAL PROGRAM,....,.,.......... * ~ * *... 51 5 1 Apparatus........... 4 o 4 * 0. 0 0...., 51 5,2 Procedure.....,.,......a. o...........O * 52 5-3 Results.* #a* a,.. -..... a a a * 4 -a a * 4 a 54 iii

TABLE OF CONTENTS (CONT D) CHAPTER Page 6 CONCLUSIONS........................................... 60 7 SUGGESTIONS FOR FUTURE WORK,..,,..,,.,.,.. 61 APPENDIX A - GRAPHS OF STREAMLINES AND VELOCITY PROFILES,,,: 4 63 APPENDIX B - COMPUTER PROGRAMSO.O.4.O,..o......... 79 BIBLIOGRAPHY.. 3 v4...........**e...............9.................... 98 iv

LIST OF TABLES Table Page 4o1 Comparison of the Computed Values of U with the Quasi-Potential Values of U at the Point R = 6.0, - 6.............................................. 37 4.2 The Location of the Radial Distance Ro from the Origin of the Point of Intersection of the Dividing Streamline with the Disc............................... 43 4.3 Computed Values for a and C for Use in Equation (4o31) 44 4~4 Computed Values for the Coefficients in the Power Series Given in Equation (4o41)...................... 46

LIST OF FIGURES Figure Page 2.1 Plot of w vso r Across the Axis of Symmetry........ 9 4l1 Pressure Gradient on the Disc as a Function of Radial Position........................................ 42 4.2 Simplified Flow Diagram for the Computer Program...... 47 5.1 Experimental Apparatus................................ 53 5~2 Superposition of the Streamlines (Experimental)....0.0. 56 5.3 Photograph of the Streamlines for o = 2.5 rad/sec and Qo = 5o0 scfm........................................ 59 A-1 Streamlines for = 1.0............................... 00 64 A-2 Streamlines for ( = 2.0,...0.... 0000..............00 65 A-3 Streamlines for 3 = 350.......................... 66 A-4 Streamlines for ( = 4,0.0..........,...ooo....oo.ooo 67 A-5 Streamlines for 3 = 50 o................................ 68 A-6 Streamlines for 63 = 60o............................. 69 A-7 Streamlines for = 700 eoo.... J o o...........0.0... 70 A-8 Streamlines for ( = 12o0............................ 71 A-9 Tangential Velocities (V) for ( = 3.0................ 72 A-10 Vertical Velocities (W) for 3 = 300...........00.00000 73 A-11 Radial Velocities (U) for B = 3,0..................... 74 A-12 Axial Velocity (W) Profiles for Various Values of ( at Z = 702......................0a....aa......... 75 A-13 Axial Velocity (W) Profiles for Various Values of ( at Z = 30 ooooooooooo.oooo......oo........ oooo.oooo........o 76 A-14 Radial Velocity (U) Profiles for Various Values of f at R = 06.............O.......................... 77 A-15 Streamlines for ( = - 2,0..,,.,...,.*....,..,....... 78 vi

CHAPTER 1 INTRODUCTION A number of studies dealing with rotational motion have been conducted in the field of fluid mechanics. In many of these studies, the flow outside of the boundary layer had circulation. Taylor(13) considered the problem of a liquid tangentially entering a chamber a distance R1 from the center of rotation and emerging through an orifice of radius R2 (R1 > > R2) after passing down a converging cone with vertex angle 2a o Rosenblat(9) investigated the viscous flow between two parallel infinite discs which were torsionally oscillatingo Stewartson studied the effect on a viscous fluid between infinite parallel discs rotating steadily with the same speeds in opposite directions~ Stewartson(l2) also studied the flow between two parallel coaxial discs rotating almost as if solid as well as the flow inside a circular cylinder which is almost rotating like a solid body. Long(7) worked on the problem of a vortex in an infinite viscous fluid and obtained both numerical and experimental results0 Earlier, Long(8) had obtained the solution for the flow of a rotating, frictionless, incompressible fluid due to a strong source or sink at the axis of rotation. The experimental and analy(4) tical results which he obtained were for a tall cylinder, Boewadt handled the problem for the flow of a viscous fluid in solid body rotation over a stationary flat disc. Barua(l) studied the case of a source in a fluid which is rotating with solid body rotation~ -1

-2The second type of problem dealing with rotational motion, for which there is no circulation outside of the boundary layer, is illustrated by the problem of a rotating disc in an infinite fluid which was solved by von Karman(l14) and Cochran,(5) In von Karman s work, which appears to be the most closely related to the present work, it was assumed that the radial and tangential velocities near the disc were linearly proportional to the radial distance from the axis of symmetry and that the axial velocity and the pressure were functions of the vertical distance above the disc. Making this similarity transformation, von Karman was able to reduce the set of partial differential equations governing the flow to a set of four coupled ordinary differential equations with the dimensionless vertical coordinate as the independent variable. The Navier-Stokes equations for axially symmetric flow, which arise for the present case of the rotating disc with a center sink, rarely can be solved analytically in closed form. For this reason, the partial differential equations governing this flow in the present work have been solved using finite-difference methods0 A very impressive demonstration of the power of finite-difference techniques has been given by Fromm in various publications (Physics of Fluids, LASL Reports) and at various meetings. At the July, 1964, International Union of Theoretical and Applied Mechanics Symposium on Concentrated Vortex Motions in Fluids held in Ann Arbor, Fromm showed motion pictures of the development of the Karman vortex street as predicted by his numerical results and as photographed experimentally.

-3The agreement was generally excellent, except far downstream where the computer results began to break down after long periods of time. This was due to the fact that he made the computation assuming a series of cylinders which were spaced far apart upstream and downstream from the cylinder under consideration and waves from the cylinder downstream finally arrived at the downstream boundary of the region under consideration, At first, one was able to observe the symmetrical pair of eddies form behind the cylinder. Then after the application of a small disturbance (to accelerate the formation of the vortex street and thus save machine time) to the vortex field, the wake lost its symmetry, and the vortex street was formed~ Grid space size naturally precluded any fine structure details, but large scale agreement for the time duration was remarkably good. The present work, in which the flow is due to a disc rotating in an infinite fluid with a sink at the center of the disc, differs from previous efforts in that in previous work if a sink was present, the fluid far away was in solid body rotation, or if the fluid was not rotating far away from the solid surface considered, there was not a sink.. For the present works the computation was carried out on the IBM 7090 at The University of Michigan0 In addition to the steady state solutions, the transient solutions were also obtained for various values of a parameter Q (= q 2 73) relating sink, 2T rotational, and viscous fources. An experimental program was carried out in conjunction with the computer study in the hope of obtaining experimental evidence to support the numerical results. Unfortunately,

-4the experimental flow regime did not coincide with the flow regime studied numerically. This was due to the fact that the experimental value of f had to be sufficiently large so that the effects of the sink could be observed, while the value of f for the numerical study had to be sufficiently small so as to not require excess computer time. However, the extrapolated experimental results appear to support the computational results, as is discussed in Section 5.3.

CHAPTER 2 THE DIFFERENTIAL EQUATIONS GOVERNING THE FLOW DUE TO A ROTATING DISC WITH A CENTER SINK 2.1 Statement of the Problem The formulation of the problem is as follows -An infinite disc is assumed to be rotating at a constant speed in an infinite fluid. The conditions are taken to be the solution as given by von Karman for the rotating disc without the sink, The sink is suddenly "turned on", and the resulting flow studied0 It is natural to assume axial symmetry about the axis of rotation of the disc, so one chooses a cylindrical coordinate system (r, @, z) in the usual manner with the origin at the center of the disc and the z-axis coincident with the axis of rotation. The Navier-Stokes equations in cylindrical coordinates for unsteady axi-symmetric flow are au au auu v2 _ 1a _ 2u _ 1 u u u2.~ "t W -ccl 3~ ---- - --- ra (+ la) 6r at r z pr ap Lr z2a rr r2 _v a v av uv [2v 2v 1 av v ( — + ~ —+ - + Vra + —+ -- (2olb) -t ar =z r z2 rLr r2 aw aw a+ w r2w a2w 1 aw -+ u +- w — 2+ + - 2 - (2 lc) at ar z p az -ar az r ar while the continuity equation is $- -O t- o(2 2) r + — + az = 0-5

-6where u, v, and w are the radial, tangential, and axial components of the velocity, respectively, p is the fluid density, v the kinematic viscosity, and p the pressure of the fluid. Differentiating the first of Equations (2.1) with respect to z and the -,third with respect to r, subtracting the results and using au aw (23) rl ~r(253) one obtains the tangential vorticity equation an,,4,4_ 2:v 6.+ +2 (2 4) + u + w v 2 (2,4) 6t 6 r z r r z _ r 6z r r r The continuity Equation (2,2) permits an axisymmetric stream function, 4, to be defined such that u 1' and w -.. (2~5) r az r ar Substituting (2o.5) into (2.3) yields 1 a2~ 1 a2~ 1 ab r az 2 r2 6r (2,6) r z r r r 6r The system of equations which were used in obtaining the solution is then given by (2.lb), (2,4), (2,5), and (2.6) and repeated here for convenience is.

-7av av av uv 2v 2v 1 v vlb) + + u +w - =2 at ar az r [r az r ar r 6a + u 6a + w an _ u+ 2 av 5S L ll (2.4) at ar az r r ar ar2 az2 r ar r2 1 a 1 a24 1 a = r az2 + r ar2 r2 ar (2~6) and 1 ab. 1 ab u = a; W=- = a j(2o5) The boundary conditions for the problem are the following - On the disc, the no-slip condition must be satisfied, i'e, the fluid must move with the disc, Hence, on the disc v = rn and u = o (237) Also, fluid must not penetrate the solid boundary; thus the disc must be a streamline, i.e, = constant for z = 0 (2o8) This ensures that at z = 0 w = 0 (2.9)

-8Along the axis of symmetry, one requires that the axis be a streamline since the radial velocity must vanish there. Hence at r 0: 0 - 0 (2,10) Also, at r = O (2.11) Along the axis of symmetry, a minor difficulty is encountered with the velocity definitions in terms of the stream functions. For instance, the second of Equations (2,5) must be taken as w(0,z) = lim (2012) r 0 r- r (2o12) Using l'Hospital's rule, (2.12) becomes a2* at r=O: w a2 (2.13) This expression must be used then in computing w on the line r = 0 o The vorticity, Rj, on the axis of symmetry due to (207) is given by at r =lim awl (2.014) hr y[~ o rJ Now if - were not zero at r =0, then the velocity profile would have a cusp at the axis of symmetry. If one were then to transfer

-9the system to rectangular Cartesian coordinates, would not be ax Figure 21 Plot of w vs r Across the Figure 2.1. Plot of w vs. r Across the Axis of Symmetry. continuous at x = 0, and hence the shear stress, which is proportional to the first derivative of the velocity, would also be discontinuous. This cannot happen in a physical system, so an additional requirement is that = at r = 0 (2.15) ar with the direct result that at r = O = O (2.16) The boundary conditions far from the disc are not as obvious. It would seem that the effects of the sink should decay at least

-10algebraically with increasing z; thus it would appear natural to require that the radial and tangential velocities be zero as in the von Karman problem, i. e, u = v = 0 at z= oo (2o17) From these and also the continuity equation, (2,2), it follows that 6w o at z oo (2.18) az Therefore W w(r) at z = oo (2 l9) The vorticity Equation (2..4) for steady state then gives XLa+ -I= 0 (2 20) ar2 r r r 2 or B aw Thus w -1 Ar2 - B n r + C (2.22) If it is required that w be bounded for large as well as small r A = B = 0 (2.23)

-11Then from (2.21) 0 = O at z = o (2.24) Finally, for large values of r, it again seems natural to require that the flow approach that given for the von Kgrman problem of the rotating disc without a sink because the sink is very far away and its influence is negligible. Hence, the requirement is that the downward velocity component, w, is independent of its radial position and that the tangential velocity component, v, is a linear function of its radial position, i.e., for large r: w - (2.25) 6r and r 0 (2.26) These are the radiation conditions used by von Kirmkn. Their use is suggested here since the influence of the sink should not be felt near the disc for large r The difference between the value of the stream function on the axis of symmetry and its value on the disc is proportional to the discharge, %, of the sink. For cylindrical coordinates, in this case, Qo 2 H [*disc Vaxis of symmetry (2.27)

-12Using (2,10) vdi = - (2,28) disc 2t 2..2 Dimensionless Form Solving the equations can be simplified significantly by a suitable dimensionless form of the system of partial differential equations and associated boundary conditions. A natural length in the present problem is related to the thickness of the boundary layer as used by von Karman. This yields Z = z r; R =r; i Ur~ = t(; - _ v - U - V; W =- (2029) where c is the rotational speed of the disc0 U V, W, and r are functions of R, Z, and T, and T is a function of R and Z o Substitution of (2.29) into the system of Equations (2olb), (2.4), (2,5), and (2.6) produces the dimensionless system.ar. a ar 6r Ur 2va V [2 ] r (2[30) T + U R + W - R R Z = L - R) av+ Uav av V (2,31) aT aR az R RJ

-13 - au aw 1 a2y 1 a~ 1 a2 r = = 2 2 + (2o32) 6Z 6R R 6R2 R2 6R R6Z2( and U= R Z W (2,33) where 2 a6 2 2 (2034) In a similar manner the boundary conditions are transformed into the following dimensionless forms: on the disc at Z = 0 V = R (2o355) Tdisc = (2q36) U = 0 (2q37) on the axis of symmetry at R = 0 T = o (2~38) U = V = 0 (2.39) r = 0 (2,40) and W = -.2 (2o41);R2 is used as means of calculating W here so that these values of W will be on hand,

-14Far above the disc at Z oo U v= 0 (2.42) r =0 (2043) For large radii at R = oo bR = ~ (2.44) av Ry O (2.45) The dimensionless system of partial differential equations with its associated boundary conditions is the form used in obtaining a numerical solution. Note that no Reynoldsi. number or other dimensionless parameter is present in the system of differential equations. The only dimensionless parameter is P which enters through the boundary conditions, One also notes that the values of 5 need not be restricted to positive values, but may be given negative values as in the case of a point source at the origin0 A physical interpretation of P is that it is the ratio of sink and rotation forces to viscous forces~ There are three simultaneous partial differential equations, (2o30) - (2.32), in the three unknowns, r, V, and' o Inspection of the boundary conditions indicates that on the boundaries Z = 00 and R = 00 there does not appear to be sufficient information0 This will be discussed in Section 3~4, The foregoing dimensionaless form is suitable if the only length scale present is engthe boundary layer thickness0 If, however, there is a characteristic length, L, present,

-15such as in vortex chamber problems, then a more suitable transformation is R= z (2.46);L LZ and U V W =; U =L; -v=; W (2.47) along with the remaining relations of (2.29) to transform the system of partial differential equations into the form ar ar ar ur 2V av [2r a2r 1 ar r -W- " - (2.48) aT aR Rz R R az [R2 az2 R R2 av av av w [a2v v2v v av v - + U- -+ W — + = Re 2v + (2249) 6,T 6R z R 6Z R aZ R R R and r = R/[R + aR32R2 (2150) where Re is the inverse of the Reynolds number and is given by Re L2 (2.51) cbL2 Likewise (2.33) becomes R 3/2 a3 Re3/2 a U W= —--; w aR (2Q52) R (Z R bR and (2..41) becomes

-16at R =0 W - R 3/2 a2T (2.53) e aR2 The boundary conditions either remain unchanged or must be modified in order to consider other problems with walls, etc.

CHAPTER 3 THE FINITE DIFFERENCE FORM OF THE EQUATIONS GOVERNING THE PROBLEM OF THE ROTATING DISC WITH A CENTER SINK 3.1 Introduction The partial differential equations that govern the problem of the rotating disc with a center sink were given in the previous chapter as ar Uar + War Ur 2V aV V2 _ 1 r (31) aT aR az R R aZ R2 6V 6V ~V UV [v 1 a7 U r +W + Ba L[ 1 V (3.2) and =1 1a - a- 1 - ar 7- - + D (353) R aR2 R2 aR R aZ2 along with the boundary conditions and the two equations yielding the radial and axial velocities from the stream functions. The problem considered by Wilkes(l5) (two dimensional convective flow in a cavity), while physically different, bears a high degree of mathematical similarity to the present problem. His vorticity equation and energy equation compare with (3.1) and (3.2), V being analogous to temperature in his case and the Coriollis and centripetal forces to internal energy and bouyancy force. This suggests a similar mathematical approach. The form of the equations in cylindrical coordinates presents difficulties for values of R approaching the axis of -17

-18symmetry (R = O) which are not present in the cavity well problem, while the infinite domain presents difficulties for large values of R, as well as Z, which are also not present in Wilkest problem. Nevertheless, the methods developed by Wilkes(15) and Fromm(6) prove quite useful in the numerical solution of the present problem. Equations (3 l) and (3.2) give the relationship of the first order time derivatives of the circumferential vorticity and circumferential velocity to the space derivatives of the same quantitieso These two equations are of the parabolic type. On the other hand, Equation (3a3) gives the vorticity as a function of the stream function and does not contain a time derivative. This equation is of the elliptic type~ Different methods of solution are indicated by this fundamental mathematical difference. The equations are next put into finite difference form0 The finite difference forms for (3.1) and (352) are used to predict the new values of vorticity and tangential velocity after an increment of time, An T Using this new value of the vorticity, the finite difference form of (3.3) is used to compute the new stream function. This process, along with the application of the boundary conditions, is repeated until steady state conditions are reached, This method was chosen in preference to other methods - say a relaxation method or a non-transient alternating direction method - because it was felt to be more economical with regard to computer time than these other methods

-19The methods used for obtaining the finite difference approximation to the solutions of the partial. differential equations are essentially the same as given by Wilkes(l5) and Fromm,(6) but are included since there are some minor differences. 3~2 System of Grid Points For the von Karman problem, the solution given by Cochran(4) indicates that the boundary conditions at Z = oo are very closely approached by the values of the solution at Z = 4.4~ With this in mind, it was felt that for use on the computer a value of Z = 7~2 was sufficiently large to use as the location of the infinite boundary in the vertical direction. For small values of P v io,e.. sink strength, this was felt to be sufficient, However, one can see from the results shown in Figure A-12 that the distance should be increased for the larger values of P. The location of the radial infinite boundary had to be at least as far away, if not farther. For a square region with grid spacings of 0o6, the total number of grid points is thus 169, With the present form of the program the time required to reach steady-state conditions for a grid of this size is 14 minutes on the IBM 7090~ An increase of R at infinity to a value of 10.8 yields a square grid of 247 points and an increase in computer time to at least 20,5 minutes in order to reach steady-state conditions. In using finite difference methods, one usually truncates the infinite series at the point for which the truncation error is

of the order of the square of the change in the independent variable, which in this case is the grid spacing. Hence the finer the grid spacing, the better the accuracy. However, there is a limit as to how fine the spacing can be made in practical computation, since eventually round-off error becomes larger than the improvement due to a finer grid, and the computer time required is at least linearly proportional to the number of grid points. It is shown in Section 4.4 that the radial position of the dividing streamline is of the order of P1/3 for small P. Thus a crude estimate of computer time would be of the order of 15(0.5 f)1/3 minutes for a grid spacing of LAR = zZ = 0.6 using the present program, a minimum number of grid points, and small P. As P becomes large, it is reasonable to suppose that the time will depend on an even higher power of Thus economics and computer capabilities limit the range of B values. For this study, it was felt that B = 12 was economically the maximum feasible value, since these values exhibit all of the basic physical aspects of the problem. 3.3 Finite Difference Approximations to the Equations of Motion The finite difference forms of the equations of motion, (3.1) - (3.3) and of the equations giving U and W, i.e., (2533),

-21are obtained in the usual manner using a Taylor's series expansion. For the vorticity and tangential velocity equations, i.e., (3,1) and (3.2), the equations are put into a form suitable for the implicit alternating direction (I:AD.) method (for instance, see Birkhoff(3)). In the method, one considers two successive half-time steps, each of length AT/2. Over the first half-step, all derivatives in the Rdirection are approximated implicitly (ioe, at the half-step) while all derivatives in the Z-direction are approximated explicitly (ioeo, at the present value)-, During the second half-step the procedure is reversed. In the following finite difference equations an asterisk(*) superscript denotes the value of the function computed at the end of the half-step in time while a prime(') superscript denotes the value at the end of the full time step.. Variables with no superscript are evaluated at the beginning of the time step. The subscripts i and j denote the grid coordinates of the point in the R- and Zdirections, respectively., Thus for given values of i and j, say I and J, the space -coordinates are R = IER Z = Jnzs (3o4) where AR and A~Z are the sizes of the grid spacing in the R- and Z-directions, respectively. In practice, (352) was first evaluated, and then, using these results, (3.1) was evaluated, Thus, for the first half-step

in time, (332) is written Vij V*i,1j i+ - Vi,j+l Vi, jl Ui,jVij AT/2 UJ ma +J U. (+) i -2Vj+l + V v 1 - 2Vi + V (,z)2 (+2 + 2__ (355) (AZ) (inR)2; AR 2AR (i/R) 2 and for the second half-step in time, (3~2) becomes 1vi v: w Vi, i+l Vi j-1 UijVij + U. +.. AT/2 iJ 2ZR iJj 2iAZ iaR V. - 2V.j + V V. 2V. + V........... ii+l,j -, i Vi ( z)2 (aR)2 AR 2AR (iAR)2 In (356) some V's appear without superscripts due to the fact that they do not appear as derivatives in (352), and hence are evaluated at the beginning of the time step, as are U and W e For the first half-step in time, (3o1) is written._ _,,rii - rii ri+.i ri-,, +l - ri, Ui,,iri,', AT/2 + Uij. + Wi j Z1T/22YkR 2z - iLR 2V v. V.r- + V j V j+l - i, j-l -i+lj - j i,j i-lj JiAR AR (Az)2 (: R)2 (307)

-23For the second half-step in time, (3.1) becomes * "r' **. r*- r*: - F' U _ AT/2 i, AR 2 2Z iAR _i,__ j__ __ __ __ _ - i, j i-i+'j 2Vr rVij+l - 21Vi, + r iAR PAR (a)2 (;AR) Since no time derivative is explicitly stated in (3.3), the finite difference approximation is given by 1 i - 2'i + d 1 - iAR (R)i 2 AR R-l, +. 1'+,j+i: ii., (3.9) iAR ((z))2 (3,9) could be solved for j and a relaxation procedure used to solve it, but the successive over-relaxation method (S.O.R.) converges more rapidly, thus saving computer time. Using the S.O.R. method, one modifies and rearranges (3.9) to the form (s+l) (s) A F Rr fint Tdif + 2An) [_i i( sR)3 gi+e b - ) - 2(1+ ) () (3.10) +r~ 2' -2- J+ 1~~~~~~'~~~'J;

where the relaxation parameter, A, lies between 1.0 and 2.0 and is best evaluated using a trial and error method on the computer, and where (~)2 D- ) (3511) The superscript (s) is the number of the iterationo Thus in (3.10) one sees that the new value of Ti.j is computed using the present values of the T's and of ij It was found by trial and error for f =5 that the optimum value of A was 1.720 and that 50 iterations were sufficient to reduce the magnitude of the change in the stream function to the order of 10 7 of the value of the stream function at that pointo All of the finite difference forms given above are used only at the interior points of the grid, i.e,, i and j always lie between their minimum and maximum values - never at these values. The values at the boundary are computed using formulas given later, The forms for (2.33) are slightly more involved, due to the differing information available on the various boundaries. Thus if the maximum values of i and j are M and N, respectively, one has for i 0 N -22 j2 U1 i,j-2 8-1 i9,j 8i +l ij+2 (312) j9 i jAR 12 AZ

-25 - for i f 0, j = 1 (adjacent to the disc), -2i0 - 3 + 6. - i.3 U i I'll i,2_ i, 3(3.13) for i # 0 j =N- 1 UiNl iAR iN 3 6i N-1 6 iN-2 +iN3 (35 14) > > R 6 z for M 2 i - 2 W1 i-?2,j - 8wi-lsJ + 8Yi+,j - Yi+2j ( wi'W AR 12 ~LR (3o15) Wij = iAR 12 AIR for i 1= -. 1 _ 2 j - 3Y1 + 62 - 3 3-1 0j lj 2 (5.16)! J AR 6 iR for i= M - 1 w-l 32w-Mj 5wY lj - 6_M.2j + wMSJ (siy W -1 (3 17) -1, j (M-l)AR 6 AR No expressions are given for U and W on the boundary j = N or for W on the boundary i = 0 o The values of U and W are needed in their explicit forms only at interior grid points since they are used in Equations (3,5) - (3,10) which involve values of U and W at interior grid points only, as mentioned above, The forms

-26for U and W on these boundaries were constructed and used in the program for completeness, but will not be given here, since they did not interact with any other part of the solution. 3.4 Finite Difference Approximations to the Boundary Conditions Having given the finite difference forms of the equations of motion, it is now necessary to prescribe the boundary conditions in a similar manner. Thus on the disc, at Z = O, one has for j = 0 V. = i.0 (3518) 130= (3019) 1~o U. = 0 (3.20) 1,0 Here, an explicit condition for Fi o is missing. It is supplied by computing the value using the condition which arises from (2,9), ioe.o at Z =0 o 0 (3,21) aR so r = T - - (3.22) z =o az z=o R z L z=o or in finite difference form ro 8 i. i j 2 (3.2.)

-27Next, on the axis of symmetry9 R =0, one has for 1i 0 uo - vo (3 24) uO Vj (3.25) rOv: o (3~26) and for completeness w.= - (3227) Far above the dis5 at large Z, one has at j = N Ui N ~ (3.28) vN 0 (3,29) N = 0 (3030) Here one notes that no condition is given for the stream function. Again) one makes use of (2~33) and constructs the finite difference form i N (3-31) 3~

-28Also) it was desired to see if the vorticity would go to zero on this boundary if its value were not fixed at zero. Here, the vorticity was computed using (3032), T 2T +! T. Y - +. i'+l,N iN N i -l,N i N il N iN. (AR)3 2i2(AR)3 + 9Ni''iN2 (3.32) 2i(ABR)(AZ)2 and the result used to check that the vorticity does indeed tend toward zero far above the disc. The check proved satisfactory, However, Equation (3.30) was used in the final computations. Finally, the boundary at large R must be considered, Equation (2,44) is replaced using (2,33) in finite difference form, Thus (2.44) becomes M-2,-... M-.l, j M2,j j 2 4M where (3~33) was obtained from the form 2 l(3534) aR R aR Equation (2*45) is first written as av v _R - )

-29and then put into the finite difference form V M Mj + 9VMlj - 4i 5VM-2 + VMj (336) The value for the vorticity at this boundary was simply that computed from 1 1-2jM + - M,j ~i R (,j M,j)2 where use has been made of (2044)) 3.5 Stability of the Method All of the finite difference approximations neglect terms at least of order (AR)2, (Z), or AT and higher. It was found that only when the ratio of the time increment to the square of the grid spacing was less than or equal to one-half was the process numerically stable, io,e, AT < 1 Age -I 2 C 1 (3,38) (grid spacing) 2 for stability. Otherwise large meaningless numbers were generated with the result that the computer "dumped" the program. The values used in this study were AT = Ol-5, AR = 15 = o 6 (3~39) giving a ratio of 5/12.

-30No criterion for the stability of the implicit alternating direction method has yet been developed. An attempt was made by Wilkes,(l5) but his analysis which is based on Fourier analysis, is valid only in the small-J* according to Birkhoff,(2) and does not involve the boundary conditions, The conclusion one may draw from Wilkest analysis is that instability was not proved, but neither was stability. Also, the original system of differential equations is non-linear, and his stability analysis is performed on a system that, in a sense. has been linearized0 *ie., near some fixed value of time; although this value can be any arbitrary value of time, it is not allowed to change by more than an amount A/T o

CHAPTER 4 NIUMERICAL RESULTS 4,1 The Dividing Streamline One notes that in all the plots of the streamlines (see Figures A-i1 to A-8 in the appendix) there is a streamline which intersects the disc and upon which the stream function has the same value as on the disc, The streamlines between it and the axis of symmetry all enter the point sink at the origin while the remaining streamlines all tend to infinity in the direction of large R 0 This streamline which intersects the disc divides the flow into two parts and hereafter will be referred to as the dividing streamline0 The vorticity on the disc is given by r (4 l) since W = 0 everywhere on the plate. Just above the disc, between the origin (sink) and the dividing streamline, the radial velocity component is negative while on the disc it is zero. Hence in this region the vorticity is negative, i,e,,o 0 < R < Ro e: r < (4,2) where Ro is the location of the dividing streamline on the disc and e is some positive number, On the other side of the dividing streamline, just above the disc, the radi.al velocity component is -31

-32positive, Thus the vorticity here is positive, ioem Ro + <; < - > 0 (4.3) Therefore, the vorticity becomes zero somewhere along the disc, If the vorticity is not zero at the point of intersection, then it must be zero elsewhere. Now if r = 0 elsewhere, then from (401), = 0 at this other point. But if 0 at, this other point, then another streamline must intersect the disc at this other point since the velocity on the disc is zero and the velocity profile starts out normal to the disc. But since two streamlines which intersect have the same value, this second intersecting streamline can only be the dividing streamline itself, i.e., they are one and the same, Hence, the conclusion is that the dividing streamline must be normal to the plate at the point of intersection. A consequence of this is that the shear stress component, Tzr, is also zero at this point since it is proportional to the vorticity along the disc. Inspection of the plots of the streamlines shows that the dividing streamline for large values of Z slopes toward the origin, and becomes more vertical as it approaches the disc. The latter is due to the centrifugal force being greater than the sink force in this region, and the former is due to the reverse of this. This influence of the sink on the flow far away is discussed in more detail in Section 4.3.

-33 - 4,2 Preservation of Circulation In viewing the plots of the tangential velocity component, V, (see Figure A-9 in the appendix) one notes that there is a bulge in the curves near the sink, i,e,, the velocity increases with decreasing Z to a maximum value and then decreases to the value of the disc4 If one looks at the streamline plots, one sees that near the axis of symmetry, a particle descends fairly vertically until it nears the disc (and the sink) at which time it begins to curve inward towards the sink and in so doing decreases Its radial distance from the axis of symmetry. While descending, it had picked up a tangential velocity component due to the viscous effect of the rotating disc, and hence had a circulation about the axis of symmetry. Looking at the second Navier-Stokes equation, (2olb), one sees that it can be rewritten as a(rv) +u a(rv) + )W a (r) 1 (r) + (rv) (44) at ar az l ar2 r 6r az2 Far above the disc u and v are both small, so for small values of r one may approximate (4,.4) by t.rv) 1 6(rv) + (rvr) (4z5) at a r2 r ar az2 In (4.5) one sees that as a particle between the axis of symmetry and the dividing streamline descends, its circulation, rv, is slowly increasing due to the viscous effects, mentioned above, and that this

-34rate of increase is of the order of the viscous terms which are very small far away from the disc. Now as this particle nears the sink, it moves very rapidly inward, thus rapidly decreasing its radial distance, r, from the axis of rotation. This very rapid change in radial position is accompanied by a rapid increase in the particle's tangential velocity component, The viscous forces initially give the circulation, Then, as the particle nears the sink, the inertia forces dominate the motion0 This then is the explanation for the bulge in the curves near the axis of symmetry, This effect was also observed experimentally when soap bubbles were used to trace the particle paths, 4.3 Influence of the Sink at Large Z As mentioned previously and as can be seen from the plots of the streamlines, the streamlines far away from the disc, ieo,, for large values of Z, curve slightly away from the vertical direction toward the sink. This can also be seen from pictures taken in the experimental phase of the program, although the effect is exaggerated by the three-dimensionality of the streamlines which were photographed. At first it seems as if the flow at large distances above the disc and at large radial distances should be entirely downward as in the problem solved by von Karmgn. However, a little reflection on potential theory recalls that, for a threedimensional sink flow, all the flow is directly toward the sink - even very far away, Hence in cylindrical coordinates, there is a

355 radial component, For a three-dimensional sink in cylindrical coordinates, if the disc is not present, the radial velocity at any point is given by U (4.6) 4t(r2 + z2)3/2 where Q is the sink discharge. In the rotating disc with a center sink, since only flow above the disc is considered 2jr(r + z 33/ Non-dimensionalizing (4.7) using (2.29) yields U = - R (4.8) Usink = (R2 + z2)3/2 While these arguments are useful for a qualitative picture of the flow, they are limited by the fact that the centrifugal forces near the disc tend to confine the source, the amount of confinement being inversely dependent on P, However, general trends are indicated. For example, for a value of P = 5.0, R = 6o0, and Z = 6.0 one has from (4.8) Usink = -0o0494 (4,9)

-36And from the von Kgrman problem as run on the computer Uvon 0.0215 (4. 10o) vonK Then U + UvonK 00279 (4.11) As given by the present computer program, the corresponding value was found to be U = - 0.0243 (4.12) Recall that in von Karman's problem, he used the transformation u = raiF(Z) (4.13) while in the present work, the transformation is u = aI U(R,Z) (4.14) Thus for conversion from one set of values to the other, one uses U = RF (4.15) Comparisons are shown for other values of P in. Table 4.1. The signs of the second and third columns agree, and the magnitudes of the velocities in the two columns are of the same order. Thus

-37TABLE 4.1 COMPARISON OF THE COMPUTED VALUES OF U WITH THE QUASI-POTENTIAL VALUES OF U AT THE POINT R = 6.o0 z = 6,o Usink + UvonK computed % difference 1 + 0o01164 + 0,00915 27 2 + 0,00185 + 0o00008 96 3 - 0oo0081 - 000827 2 4 - 0,0178 - 0.o0166 7 5 0.0279 - 0,0243 13 6 0.0375 - 0.02912 22 7 0- 0475 - 0o0391 18

-38this model can be used to approximate the location and shape of the dividing streamline for f of such values that the confinement is not too severe (i.e., for sufficiently large B ). 4.4 Influence of the Parameter j Looking at the results for W, the axial component of the dimensionless velocity, one sees that the values of W near the axis of symmetry increase with increasing P, while the values of W away from this axis change very little. For a given value of B, the dimensionless dividing streamline intersects the disc at some fixed R * However R and f are both dimensionless. As before i=Qo2 (4.16) and R = r (4.17) Experimentally, Qo and cD can be varied independently. Thus it is possible to obtain the same value of f with numerous sets of values for Qo and X. Suppose one decreased Qo and increased o in such a manner as to give a constant value of B. Then, since R

-39here is constant, from (4.17) one can see that the physical radial position, ro, of the dividing streamline moves inward toward the axis of symmetry, while from w - VM W (4, 18) one can see that the downward physical velocity component, w increases in magnitude, If, instead, n were to be decreased and Q increased, the equations would predict that the reverse occurs, Now near the axis of symmetry, the sink dominates the flow, Since the dimensionless dividing streamline moves radially outward for increasing values of P, holding the rotational speed u constant while increasing Qo means that more fluid flows into the sink. Thus near the sink the physical velocity components u and w should both increase. But if X is constant, then from (4o18) and u 47 U (4.19) U and W must increase in order to give an increase to the physical velocities, u and w o The computer solutions show that this is indeed what happens, However, next suppose' is increased by increasing w with Qo fixed, As before for increasing B, the computer shows that both U and W increase in this case4 Then from (4j,18) and (4).19) both u and w must increase due to the larger U and W in addition to the increase in c o

-40Now from the computer data, the dividing streamline position, Ro, is some function of P,i oe., Ro = f(p) (4~20) and from (4.17) the physical position, ro, of the dividing streamline is given by rO ff f(f) (4.21) or r f Y2 - (4.22) From (4.22), one can see that if the function f(P) is of any degree in P less than unity, then ro will decrease with increasing c o If this is true then the physical components, u and w, must increase with c (Qo fixed) as predicted above. Now from the Navier-Stokes equation, (2.1a), when viscosity is neglected, one has for steady state - = P a au v2+w (4,23) ar ar +z r On the disc this reduces to =P [ au v2 (424) ar ar r

-41If one now considers the potential flow of a fluid which is in solid body rotation with a sink on the axis of symmetry at the center of the disc, one finds that (4.24) may be written as t ~P [2:2r4r- ](4.25) If one plots the separate contributions of the sink and solid body rotations to the pressure gradient, as well as the total pressure gradientversus the radial position on the disc, as is done in Figure 4.1, one sees that indeed the pressure gradient has a minimum which is non-zero. The minimum pressure gradient is found by differentiating (4.24). Now for the minimum pressure gradient, differentiating (4.24) with respect to r and setting the result equal to zero yields a [ u -v2] =0 (4.26) Using (2.29) this becomes 2 2U 2V V V2 u)2 + 2u 2V av v = (4.27) - 6R 6R2 R 6R R Putting (4.8) and (2.35) into (4.27) yields 1 i- i (4.28)

-42p ar 2 Vy- total ( 2r 2 ) /l J\ \ i <sold body rotation'(o sink ( -r r Figure 4.1. Pressure Gradient on the Disc as a Function of Radial Position for Inviscid Flow.

-43 - solving for Ro Ro = 646O 1'3 (4.29) or Ro = 1,4681 (4o30) which approximates the location of the minimum pressure gradient on the disc. Looking at the data for the position of the dividing streamline and the corresponding values for P in Table 42-2, suppose one makes the simple assumption that R = Cs (4o31) then ln Ro = alnp + in C (4.32) TABLE 4~,2 THE LOCATION OF THE RADIAL DISTANCE Ro FROM THE ORIGIN OF THE POINT OF INTERSECTION OF THE DIVIDING STREAMLINE WITH THE DISC [B ~Ro 1 1,49992 2 1,83384 3 2. 10654 4 230103 5 2.45902 6 2,62331

-44Then c is found from in (R)2 in (R ) a = oo (4.33) in (P)2 - In (p)1 The range of a for the data obtained from the computer is given by Table 403 and is 0..28998 = = 0354737 (4o34) So one may try = 1/3 (4,35) Thus it appears that ~Ro ~~ P~~ 113 / ~(4,36) is a reasonable approximation for small P, and hence the above conclusions based on the assumption that a < 1 are substantiated, TABLE 4,3 COMPUTED VALUES FOR ca AND C FOR USE IN EQUATION (4.31) C 1 0,289980 1.49992 2 0..341917 1,44689 3 0.306958 1.50354 4 00297597 1.52318 5 0~354737 1*38935 6 0.o324893 1.46566

The values of Ro and C were found using the computer program EXPON,001o The location of Ro was based on the data at Z = 0o,6, i.eo,, it was assumed that the streamline intersected the disc in a vertical line which was straight for a distance of at least o06 above the disco Similarly, C was found to be C _ 1047 (4.37) so that one has Ro _ 147 1/3 (4~38) In view of the above results, it is interesting to note that for an infinite, inviscid, incompressible fluid rotating with solid body rotation about an axis on which a sink was located, Long(8) found that all the flow going into the sink came from within a region bounded by a cylinder of radius = o I/ (4-39) In the same paper, Long quotes a personal communication from Sir Geoffrey Taylor in which Taylor gives r = Qo 1 (440) for a jet due to a source at the axis of a rotating fluid.

-46If now, instead, one tries to express 5 as a power series in Ro ioeo = Z an Rn (4.41) n one finds that the coefficients an are as shown in Table 4.4. Thus the first seven terms in (4Q41) are given by = 1362,74 Ro - 3962 5 R2 + 4753.88 R3 - 3012o45 R4 + 1063,96 Ro - 198..628 Ro + 153175 42) TABLE 4.4 COMPUTED VALUES FOR TEE COEFFICIENTS IN THE POWER SERIES GIVEN IN EQUATION (4o41) n a n 1 13624 74 2 - 3962.5 3 47535488 4 - 3012.45 5 1063596 6 - 198.. 628 7 1543175 4,5 Technique of Solution As mentioned previously, the solution was carried out on the IBM 7090 computer at The University of Michigan. The program was broken down into the main program and five subroutines. A simplified program flow diagram is shown in Figure 4.2.

-47Read and Print Data Initialise Stream Functions and Velocities (Begin.) Set Time Equal to Zero Turn on Sink Increment Time Compute EI Using Successive Over-Relaxation (AOPT.) Compute U & W (SUBUWG.) Compute V (SUBV.) I Compute' (SUBG. ) -::a~~~ ITrue Check Time T < T tFalse Print Final Values U. V, W, A, FY, T. End of Program Figure 4,.2. Simplified Flow Diagram for the Computer Program.

-48In more detail, the main program first, in the subroutine BEGIN., initializes the values of the stream function, i.e4. the von Karman stream functions are calculated, The main program then transfers to the subroutine SUBUWGo where all values of U and W are calculated as well as the vorticity on the disc, The program next "turns on" the sink and then moves to the subroutine AOPT, where the new stream functions are computed from Equation (3,10) using the successive over-relaxation technique mentioned earlier. Using these new values of the stream function, the program returns to SUBUWG, After computing the new values for U and W and for r on the disc, in SUBUWG~., the program computes the new tangential velocities, V, using (3~5) and (3o6) in the subroutine SUBVQ Finally, using (357) and (3.8) the new values of the vorticity, r are computed in the subroutine SUBG. The program then computes the values for the variables on the boundaries which are allowed to change and returns to AOPT, to compute the new stream functions, and the cycle starts again, The process is repeated until a pre-set time is reached, In applying the boundary conditions where the value was allowed to change with time, ieo, the stream function at Z = o or at R = o } it was found that if the values were recomputed at the end of each iteration in the successive over-relaxation process (AOPT,) the process became unstable and large meaningless numbers were generated. If, however, the values were recomputed after each time step, then the process appeared to be stable with meaningful results0

-49Also, in applying the conditions at Z = X, three methods were tried, First, U was set equal to zero, and the corresponding stream function was calculated using (3o31). The values of r on this boundary were then computed on the basis of the stream functions0 The second method was to set r equal to zero on this boundary and calculate the stream function here on this basis, The value of U here corresponding to the stream function derivatives was then computedo It was found that the first method was slightly superior to the second in that the values of both U and r were smaller adjacent to this boundary. However, a third method was the one used. This was simply to set U = 0 on this boundary and calculate the stream function using (3531)o In addition, F was set equal to zero here, This is the more rigorous approach, but the other two approaches were tried in order to see if less rigid boundary conditions would possibly result in a good approximation in a shorter time. The third approach gave the smallest values of U and F near this boundary for a given amount of computation time. Due to the fact that the grid spacing was large and the number of grid points was relatively small, the results obtained on the computer are only indicative of the flow behavior for the problem of the rotating disc with a center sink. Better quantitative results would be obtained with a much finer grid spacing, many more grid points and longer runs on the computer,

-504.6 The Rotating Disc with a Center Source The case of a source in a fluid rotating with solid body rotation has been considered by Barua.(l) In his problem, the fluid is rotating about the x-axis and the source is replaced by a discharge through a sphere x2 + y2 + z2 = a2 (4~43) located at the origin~ However, the flow from the source is assumed to remain irrotational, thus giving rise to a discontinuity at the interface of the irrotational and rotational fluids, In the present case of the rotating disc with a center source, no assumption of irrotationality is made, As mentioned previously, the computer program was set up such that it could be used with a source instead of a sink, This is done simply by using a negative value for P o The result of using Pi= - 2,0 is shown in Figure A-15 in the appendix.

CHAPTER 5 EXPERIMENTAL PROGRAM. 5 1 Apparatus The experimental program was carried out in the Fluids Laboratory Building at the North Campus. A room 12 feet long, 10 feet wide and 11-1/2 feet high was constructed using 2-inch x 4-inch beams covered with masonite on the sides and a thin polyethylene sheet 04008 inch thick on the ceiling. All joints, cracks, and holes were then sealed using masking tape. The experimental work was done inside this room in order to eliminate any drafts. Inside this room, a 1/4 inch thick piece of tempered masonite approximately 4 feet in diameter with a 3/8 inch hole at its center was mounted in a horizontal plane on a vertical hollow shaft which could be driven by a motor. The hollow shaft was connect-, ed to a vacuum pump located outside of the roomo The connecting line between the disc hole and the vacuum pump contained two Cox flowmeters connected in parallel used to measure the discharge, The flow was controlled with metering values mounted in series with the flow meters. The drive unit was a Servo-Tek thyratron regulated drive unit. Thus the disc speed could be varied as well as the sink discharge. The streamlines were visualized using a smoke generator which fed smoke into the region above the disc through fifteen 1/8 inch diameter glass tubes spaced two inches apart. The smoke was generated by blowing air through a 1 inch diameter pyrex tube containing a burning oil-of-wintergreen-impregnated cigar4 -51

A second method of visualizing the flow consisted of dipping eight 19 gage hypodermic needles, mounted in parallel at the end of the air line, into a glycerin-soap solution. raising the needles far above the disc, and blowing small bubbles into this region. The path lines of the bubbles were then observed. Tufts of wool were also tried, but even very fine threads were too heavy and too stiff to be of any value, A cutaway drawing of the experimental set-up is shown in Figure 50ol 502 Procedure It was found that the experiments could be carried out only in the evening with the lights turned off above the room~ If there was a strong' breeze blowing outside, there was a noticeable effect inside the room. This was attributed to a cooling effect of the subsequent drafts inside the laboratory passing over the ceiling, which was only a thin plastic sheet, and to possible entry of these drafts into the room~ The lights above the room caused a very noticeable convection current in the room when lit. In order to obtain photographs, a spotlight was aimed from outside the room through. a glass window0 The spotlight beam was at right angles to the direction from which photographs were taken0 To minimize the heating of the air inside the room from the spotlight, it was switched on only during the period when the camera was in use0 The procedure for making a run was to turn off the lights and blowers in the portion of the building near the room0 Ten to

-53MIRROR SMOKE MAN I FOLD PYREX TUBEHYRATRONDRIVE UNIT / CONTAI N ICNGT\L v SPEED CONTROL CIGAR AIC SPOTLIGHT FLOVACUUM PUETERS CAMERA PRESSURE REGU LATOR PYREX TUBE CONTAINING CIGAR AIR SUPPLY LINE VACUUM PUMP FOR SINK Figure 5.1. Experimental Apparatus.

-54fifteen minutes were then allowed for the air inside the room to become quiescent. The vacuum pump was then turned on along with the smoke generator. The smoke was introduced at a velocity comparable to that of the surrounding air so as not to disturb the flowo A short time later. the disc was set in motion and pictures were taken. After about a minute of runningf the apparatus was turned completely off, and a period of ten to fifteen minutes allowed to pass before making another run, This was done to minimize turbulent effects due to the nearness of the walls of the room to the disc. This wall effect also limited the speeds at which the disc could be run for even short periods of time. 5,3 Results Pictures taken of a vertical plane of the flow show some change in the streamlines with changes in the variables Qo and X. These changes are indicative of what is happening but do not yield good quantitative results. This is due to the fact that the actual streamlines are three-dimensional while the photographs are plane. Photographs taken looking down on the disc gave little detail because the streamlines were fairly vertical until they neared the disc where they were quickly dissipated due to the increased velocity of the air in this region0 Actual observation, just as in the attempts made using soap bubbles, appears to be the best method at the present time0 With the disc speed, c 9 fixed at its minimum value one can see the dividing streamline move outward as the sink discharge, Q0, is

-55increased. It is more difficult to see the dividing streamline move inward when Qo is fixed at its maximum value (5 scfm) and X increased~ It is almost impossible to see the dividing streamline move when Qo is fixed at a value of 2~0 scfm or less and w is varied. Actually, one does not "see" the dividing streamline at all; one only observes the streamlines in the vicinity change their slope near the disc as Qo or X is varied, The effect is noticeable in the region 0 = r < 6 inches. The streamlines farther away move so little that it is difficult to tell whether they are being affected by the changes in Qo and c or by the little disturbances that are undoubtedly present in the roomo The fact that the disc was not perfectly plane could be seen at high rotational speeds (1.5-20 rad/sec) as pulses in the smoke streams that caused a wavey form, When a large amount of smoke was present in the room, it was possible to see the air above the disc being pulled downward as the disc was started up from rest. The readings taken directly from the flowmeters were felt to be sufficiently accurate for this work since the maximum pressure drop before entering the flowmeters was 0O55-inch of mercury. This coupled with a correction for changes in atmospheric pressure amounted to less than a 4 percent error in the sink discharge values. The superposition of the streamlines for w = 245 rad/sec, Qo = 5,0 scfm on the streamlines for X = 4.0 rad/sec, Qo = 2.5 scfm is shown in Figure 52,4 Here one sees that due to the three

3 = 5470 -. = 8650 I I.I SINK / DISC Figure 5.2. Superposition of the Streamlines (Experimental) for P = 8650 (w = 2.5 rad/sec, % = 5 scfm) and P = 5470 (a = 4.0 rad/sec, % = 2.5 scfm). The Streamlines are 2 inches Apart at the Top.

-57dimensionality of the system, the dividing streamlines are difficult to detect, but the qualitative effect is immediately apparent. Now in the experiment 1,5 0 scfm (501) 2.5 < <a= 5.5 rad/sec (5.2) and v Z 1,8 x 10-4 ft2/sec (5.3) Using (4~.30) along with (2.29) one finds ro = 1o.468 (o) /t = 176 (5.4) (One notes here that ro is independent of the fluid used, ioe.., the results should be the same for any fluid.) Thus from (501), (542), and (5,4) 1.99 r 3 08 in (545) is the predicted range for the experiment. However, the experiment shows that (see Figure 5..2) 6 r< 10 in,(56) indicating that indeed (5.~.4) is valid only for small ~ o

-58Similarly, from the von Karman results computed by Cochran one sees that the boundary layer thickness (based on the fact that the radial and tangential velocities are very small) is approximately of the dimensionless height z 6 60o (5.7) or 0413 = z = 0,61 in (5,8) is the predicted range of the boundary layer thickness, However, as can be seen from Figure 5r3, the effect of the rotational motion of the disc affects the flow at considerably larger distances of the order of two or more inches,

C) 0) tVf' c4 J~J 0 U~ 1111 I I 11 g 10 _ E~~~~(3 I Z | | | |, | | | | | | | | -~~~~~~~~~~~~~~~~~~~~~~~~~~~~,f.. o ( _1) S~ -4 -4 I I sll I I I_ t.Ns q~.... tN 4f........ sf —{~~~~~~~~~~~~~.. _ Z R R S S _ S hI~~C~ ~ __ Z R R | S S \M~~~~~~~~~t~ f _ | | | | | E | ~ ~ ~ ~ t~J c 1:~ ~ ~ i I11I1 _ 40. _ | | | | | | _ tf;4~~~~~~~~~~~~~~~~~~~~~~~~r~~ P~~~~~~~~~~~~~~_ 0 ||| ~~~~ii~~~~~~~i~~~~~ _SS8 R~R -..,:::_ _:- I GR

CHAPTER 6 CONCLUSIONS Summarizing the results, it is found that (1) the numerical study indicates that for small values of f the non-dimensional radial distance of the dividing streamline from the axis of rotation is proportional to the cube root of the parameter f, which is a measure of the ratio of rotational and sink attractive forces to the viscous force; (2) the flow fairly far away from the sink at rather large heights above the disc is still influenced somewhat by the.sink as evidenced by the streamlines curving slightly toward the sink in this region; (3) the dividing streamline intersects the disc at right angles; and (4) the location of the minimum pressure gradient along the disc appears to be at the intersection of the dividing streamline and the disc, A phenomenon best described as a tendency to preserve circulation is demonstrated both experimentally and numerically, although the effect is diminished near the disc by viscosity. The circumferential velocity of a fluid particle located between the axis of symmetry and the dividing streamline first increases due to the viscous effects as it approaches the disc. This gives rise to a circulation, As the particle is drawn towards the sink, its radial distance from the axis of rotation decreases. Then, due to the predominance of the convection effects over the viscous effects, the circumferential velocity increases. -60

CHAPTER 7 SUGGESTIONS FOR FUTURE WORK The results presented in the present work are indicative of the general behavior of the flow but the fine structure has not been obtained. A finer grid spacing would improve the accuracy of the solution and would serve as an indication as to whether the method converges to the correct solution (i.e., if decreasing the grid spacing by an order of magnitude gives the same significant figures, one would be optimistic about the correctness). Also, the problem of the rotating disc with a center source can be very easily investigated with no additional experimental equipment required and using the present computer program. As mentioned previously, the computer program is also set up to include vortex chamber problems. Specifically, it is possible to consider a cylindrical chamber in which the wall is stationary or rotating, the top is stationary or rotating, and the bottom is stationary or rotating and has a source or sink at its center. Each of these conditions may or may not be present independently of the others. It appears from the results of the experimental program that the use of a much larger room to better ensure a quiescent atmosphere would be very desirable. A stronger pump would enable increasing the range of 5. Also, a better method of measuring the position of the dividing streamline should be devised. One might try finer streams of smoke more closely spaced and use a telephoto lens on the camera. The main difficulty seems to stem from the

fact that the air velocities are very small and hence any disturbances, even though small, affect the flow considerably. A suggestion prompted by Equation (5.4) is to use a different fluid, say water, because the location of the dividing streamline is independent of the fluid used. A small diameter disc, say about 12 inches, could be immersed in a very large diameter tank, say 10-12 feet, such as is available in the Fluids Laboratory on the North Campus. This would allow use of better flow visualization techniques, but would still require long waits between runs, and would still be dependent on room temperature variations~

APPENDIX A GRAPHS OF STREAMLINES AND VELOCITY PROFILES Only the graphs of the velocity components for P = 350 are shown here. For other values of B, the graphs have the same general shape, The graphs of the streamlines for all values of f are also shown. -63

-640.5 1.0 3.0 5.0 10.0 15.0 20.0 12 I0 9 N 6 115 1 2 3 4 5 6 7 8 9 10 11 12 SINK RADIAL - R Figure A-1. Streamlines for 3 = 1,0.

-65 - 0.07 1.0 2.0 5.0 10.0 15.0 20.0 12 I I II 10 9 8 7 NI ~6'~ 5 4 3 I 2 3 4 5 6 7 8 9 10 11 12 SINK RADIAL —- R Figure A-2. Streamlines for, = 2.0.

-660.07 1.0 2.0 3.0 5.0 15.0 20.0 II 10 9 8 7 =. 3 SINK- I 2 3 4 5 6 7 8 9 10 11 12 RADIAL- " R Figure A-3. Streamlines for B 3.0.

0.0 03 1.0 3.0 4.0 5.0 10.0 15.0 12 II 10 98 7 6 — _I~~~~~~~~~.543 2 SINK 1 2 3 4 5 6 7 8 9 10 II 12 RADIAL R Figure A-4. Streamlines for 8 = 4.0.

-680.0 1.0 13. 5.0 10.0 15.0 10 9 -" -, 6 17 6 SINK- I 2 3 4 5 6 7 8 9 10 11 12 RADIAL -- R Figure A-5. Streamlines for B = 5. 0.

-690.0 2.0 4.0 6.0 10.0 15.0 12 II SINK 3 4 5 6 7 8 9 -0 1 2 10SINK 2 3 4 5 6 7 8 9 10 12 RADIAL - R Figure A-6. Streamlines for 3 = 6.0.

-700.0 a2.0 4.0 7.0 10.0 15.0 12 0-0 10 9 6 N 5. 4 SINK I 2 3 4 5 6 7 8 9 10 11 12 RADIAL — R Figure A-7. Streamlines for P = 7.0.

-710.07 3.0 5.0 10.0 120 15.0 12 I1 IO 9 8 7 6 5 4 4 2 I 2 3 4 5 6 7 8 9 10 II 12 SINK RADIAL —-- R Figure A-8. Streamlines for f = 12.0.

1 2 I~ 2 - -i —--- I —— r — r I, I I i I! a II I0 9 8 7 FI I 4 4 4 3 2 - 0 0 I 2 3 4 5' 6 7 8 9 10 II 12 13 14 15 16 17 18 19 20 21 22 23 24 25 RADIAL —- R Figure A-9, Tangential Velocities (V) for B - 3,0.

-73 - 12 Il 10 9 8 -IJ I23456 7 8| 0 I1 4 tl I ) 5 4 3 2 1 2 3 4 5 6 7 8 9 10 II1 2 RADIAL - R Figure A-10. Vertical Velocities (W) for, = 3.0.

I0 10 9 8 5 4 3 2 2 3 4 5 6 7 8 9 10 - 12 RADIALFigure A-11. Radial Velocities (U) for P = 3.0.

-750.0 - 0.8 I. oI I I I I I I I I I I 1 2 3 4 5 6 7 8 9 10 11 12 R Figure A-12. Axial Velocity (W) Profiles for Varioug Values of B at Z = 7,2.

-760.0 -2.0 I 2 3 4 5 6 7 8 9 10 11 12 R Figure A-13, Axial Velocity (W) Profiles for Various Values of P at Z = 3.0.

121 lOI 5 -3 7.0 3.0/ / 1~~~~~~~~~~~~~~~.0/ 5.0 1 -7.0 -6.0 56.0 4.0 -3.0 -2.0 -1.0 0.0 Figure A-14. Radial Velocity (U) Rrofiles V~osVle fPa.6 Varius alue orB.tR t0. 4

-780.0 1.0 5.0 10 9 - 8 t 645' 4 -220 SOURCEI 2 3 4 5 6 7 8 9 10 11 12 RADIAL - R Figure A-15. Streamlines for f = -2.0.

APPENDIX B COMPUTER PROGRAMS The computer programs are in the MAD language and are given here. The data which must be submitted along with the main program (MAINoOOl) and its subroutines is as follows: Boolean Symbols BB Sets the value for BB1 and BB3 after the first time step, BB1 Prints the current values of GAMMAS and GAMMAPo BB2 Transfers to START and transfers to NINTH. BB3 Prints the current values of VSTAR and VPRIMEo BB4 Prints the current values of A, PSI, and DELTA.,. for each iteration. BB5 Transfers to FIRST. BB6 Increases the value of ITMAX and calculates PSI for a more accurate estimate~ BB7 Prints current values of U and Wo BB8 Transfers to NINTH. Reads and prints previously computed values of BETA, PSI, U, V, W, and GAMMA from a punched deck obtained from a previous run, BB9 Prints U, V, W, GAMMA, PSI every DIT time step0 BBFINE Changes to a finer grid. BBFIX Prints I, J, PSI, U, V, W, and GAMMA for fixed increments of I and J at intermediate time steps. BBMAIN Prints and punches I, J, PSI, U, V, W, and GAMMA for steady state. BBPSI Prints PSI and DELTA for each time step. BBSLCT Prints I, J, PSI, U, V, W, and GAMMA for selected values of I and Jo -79

-80BBSKIP Transfers to NINTH and then to SECOND. BPUNCH Punches steady state values. BVTOP Sets the values for V(I,JMAX). BZEROU Sets U(IMAX,J) = O. BZEROV Sets V(IMAX,J) = 0O BZEROW Sets W(IMAX,J) = 00 GAMTOP Computes GAMMA(IJMAX) VSTILL Sets V(I,O) = O. VXI Sets the values for V(IMAX,J). Data A Optimization parameter for the successive over-relaxation technique. BETA Value of the stream function on the dividing streamline. BETMAX Maximum value of BETA, DBETA Increment of BETA, DIM(O) = 2 JMAX + 3, JMAX + L, DIT Increment for ITER. DITMAX Increases ITMAX and gives the maximum value for ITER, DTAU Increment for time, DXI Radial grid spacing. DZETA Vertical grid spacing, EPSI = 0,002 G Values of G from the von Karman problem. H Values of H from the von Kamrman problem.

-81IEND Last value of I for which intermediate values are printed. IFOUR Selected value of I for which intermediate values are printed, IMAX Number of grid spacings in the radial direction. IONE Selected value of I for which intermediate values are printed.. ISPACE Increment of I for which intermediate values are printed, ISTART First value of I for which intermediate values are printed. ITHREE Selected value of I for which intermediate values are printed. ITMAX Maximum number of iterations for the successive overrelaxation technique. ITONE First value for ITER used in selecting the time step for which the intermediate values are to be printed. ITWO Selected value of I for which the intermediate values are printed. JEND Same as IEND, only for J. JFOUR Same as IFOUR, only for J, JMAX Same as IMAX, only for J. JONE Same as IONE, only for J. JSPACE Same as ISPACE, only for J, JSTART Same as ISTART, only for JO JTWO Same as ITWO, only for J, RE Inverse of the Reynolds number for the problem. TAUMAX Time for which steady state is reached.

-82$ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT MAIN.001 RFLOW DUE TO A ROTATING DISC WITH A SINK PROGRAM COMMON U,V,VSTARVPRIME,W,GAMMA,VSTILL, RE, PSI,AAA 1,CEEDEE CADEE PBCAB,CPC,DXI21,DXISQ,IDXDZIDXICU,VER2,MINUSI,PL 2US2I,OIDXI OIDXIS, H,G,A, IDXI,BB,DTAUD2,DXISQ,DXISQ2,DZE 3TSQ,DZETA2,DZETS2,DZSQIDZETS4,PL2ARE,AOARE,DZETA6,DZET12,DXI 412,DXISQ1,DXI2,DXISQH,ARE,EPSI,BETA,DXI,DZETA,I,J,ITMAX,ITER 5,,B1'BB2,BB3,BBFINEISTARTJSTART, IENDJEND 6,BB4,BB5,BB6,BB7,BB8,BB9 -7,DIM,ANSW1'ANSW2,ANSW3,BETMAX,DBETA,DTAU,OLDATA,TAU,TAUMAX,DI 8TMAX,K,BEEIMAX,JMAXL.IMAXLJM AX,ENDBET 9,BBP.SI,ENDTAU,SPACE DIMENSION U(1900,DIM),V(1900,DIM),VSTAR(1900DIM),.VPRIME(I9UU 1,DIM),W(1900,DIM),GAMMA(1900, DIM),GAMMAS(1900,DIM),GAMMAP(190 20'.,DIM),PSI(1900,DIM), DELTA(1900, DIM), AAA (250),CEEE(250),DEE(25 3'0 ),CAPB(250.),CAP-C ( 25-0),DXI2I(250),.IDXISQ(250),IDXDZ(250),'IDXI 4CU(250),OVER2I(250),MINUSI(2'50),PLUS2I(250),OIDXI(250),OIDXIS 5(250),H(50),G(50),A(l11) IDXI(250) 6,DIM(3) EQUIVALENCE.(.VSTAR(0),GAMMAS(0)),(VPRIME(0),GAMMAP()), 1(GAMMAP(0),DELTA(0)) INTEGER IMAX,JMAX,DITMAX,LIMAXLJMAX BOOLEAN BPUNCH BOOLEAN BZEROW,VXI,BZEROV, BZEROUGAMTOP,BVTOP BOOLEAN BBFINE BOOLEAN BBSKIP INTEGER IE-ND,JEND INTEGER IONE,ITWOITHREE,IFOUR INTEGER JONE,JTWO,JTHREE,JFOUR INTEGER DIT,ITONE INTEG.ER ISTART,JSTART,ISPACE,JSPACE INTEGER I,J,KITMAX,ITER VECTOR VALUES SPACE=$5HDXI= E13.6,7HDZETA= E13*6*$ VECTOR VALUES ENDTAU=$5HTAU= E13.6,6HDTAU= E13.6*$ VECTOR VALUES OLDATA=$(2I3,5E13.6)*$ VECTOR VALUES ENDBET=$6HBETA= E13.6*$ VECTOR VALUES-ANSW1=$(S2,I3,S1, I3,5E16*6)*$ -VECTOR VALUES ANSW2$(2 I 3,5E13.6)*$ VECTOR VALUES ANSW3=$1H69S3,1HI S3,1HJSO1,3HPSI,S13,1HUS14, 11HV,S15, lHW,514, 5HGAMMA//*$ BOOLEAN BB,BB1,BB2 B B 4 B B 5, BB6 BB7BB8'BB 9 BOOLEAN VSTILL BOOLEAN BBMA I.N,BBF I X, BBSLC.T BOOLEAN BBPS FIRST READ AND PRINT DATA WHENEVER BB8 READ FORMAT ENDBET,. BETA ~READ FORMAT ENDTAU,TAU,DTAU READ FORMAT SPACE,DXI,DZETA PRINT RESULTS BETA PRINT RESULTS TAU,DTAU,DX.I,DZETA THROUGH OLDA,FOR I=01, I.G IMAX THROUGH O.LDA, FOR J=1O,1,J.G.JMAX OLDA READ FORMAT OLDATA,I,J,PSI ( I,J),U( I,J),V( I,J),W( I,J),GAMMA( I, 1J) PR INT FORMAT ANSW3 THROUGH OLD, FOR I=0,1,I*G.IMXAX THROUGH OLD, FOR J=O,1J,J,GJMAX OLD PRINT FORMAT ANSW1,I,.i PSI'( I,J),U( I:,J),V( I,J),W( I,J),GAMMA(

-83 - 1J) END OF CONDITIONAL RCOMPUTE CONSTANTS WHENEVER BBFINE DTAU=DTAU/4,0 DXI=DXI/2.0 DZETA=DZETA/2*,0 END OF CONDITIONAL DTAUD2=2. /DTAU DX I SQ=DX I*DX I DZETSQ=DZETA*DZETA RCOMPUTE VALUES FOR A FINER GRID -WHENEVER BBFINE EXECUTE BEGIN. PRINT FORMAT ANSW3 THROUGH FINE,FOR I=O1,IIGIMAX THROUGH FINEFOR J=Oi1,J*G.JMAX FINE PRINT FORMAT ANSW1 I,JPSI (I,J),U(I.,J),V(I,),W( I J) 1GAMMA(I,J) TRANSFER TO SKIP END OF COND.ITIONAL IDXI (O)=O.O DXI2I(O)=O * IDXISQ( 0 )=O.O IDXDZ(OC)=O~. IDXICU( 0 )=O.O OVER2 I ( O)=O O MINUSI (O)=O.O PLUS2I(O)=O.O OIDXI (O)=O.0 OIDXIS(O)=O.O THROUGH REPEAT, FOR I=11I.G.IGIMAX IDXI (I)=I*DXI DX12I(I)=:1.O/(2.0*IDXI(I)-*DXI) IDXISQ(I )=IDXI ( I)*IDXI (I IDXDZ(I)=IDXI ( I )*DZETA IDXICU( I )=IDXI(I )*DXISQ OVER2 I ( I )=.O/(2 O I ) MINUSI(I )=10.-OVER2I(I.) PLUS2I (I)= 1.O+OVER2I (I) OIDXI( I)=1.O/IDXI (I) REPEAT OIDXIS(I )=OIDXI(I)*OIDXI(I) SKI P DX ISQ2=2.0/DX'I SQ DZETA2=2.0*DZETA DZETS2=2.0C*DZETSQ DZSQI=1.O/DZETSQ DZET54=4.0/DZETS2 ARE=DXI SQ/DZETSQ PL2ARE=2.0* ( 10+ARE) DZETA6=6. 0*DZETA DZET12=12. 0*DZETA DX I 12=12 O*DX I DX ISQ1=DXI-SQ2/2 ~ 0 DXI 2=2.O*DX IXI' DXI SQH=DX I SQ1/2. 0 DXI I=1.O/DXI LIMAX= I MAX-1 LJMAX=JMAX-1 WHENEVER BBFINE

-84TAU=O 0 TRANSFER TO NINTH END OF CONDITIONAL WHENEVER BBSK IP TRANSFER TO NINTH R COMPUTE PSI, V, OR VALUES FOR A FINER GRID EXECUTE BEGIN. WHENEVER BVTOP THROUGH TOPV,FORI=1 1 I.E. IMAX TOPV V( I,JMAX)=IDXI( I ) END OF CONDITIONAL WHENEVER BZEROW THROUGH.WWALL,FOR J=,1, J.G.JMAX WWALL W(IMAX,J)=O.O END OF CONDITIONAL WHENEVER VXI THROUGH BVXI,FOR J=1,'1,J.G,*JMAX BVXI V( IMAXJ)=IDXI (IMAX) OR WHENEVER BZEROV THROUGH VEZERO,FOR J=1,1,J.GJMAX VEZERO V(IMAX,J)=O.0O OTHERWISE CONT INUE END OF CONDITIONAL WHENEVER BZEROU THROUGH UEZERO, FOR J = 1, 1, J. G -JMAX UEZERO U(IMAX,J)=0.O END OF CONDITIONAL FOURTH TAU=O.O WHENEVER BB2, TRANSFER TO START NINTH BB2=OB WHENEVER BBSKIP,TRANSFER TO SECOND THROUGH EIGHTH,FOR I=O,1,I.G.IMAX EIGHTH PSI( IO i 0=BETA WHENEVER BZEROW.AND.BZEROU THROUGH ZEROW,FOR J-=ltl,J.G.JMAX ZEROW PSI ( IMAX,J ) =BETA OTHERWISE JMAXSQ=1.0/ (JMAX*JMAX) THROUGH ENDPSI.FOR J 1,1,-J.G,.JMAX DELTA(IMAX,J)=PSI (IMAX,J)+ (1.O-J*J*JMAXSQ ) *BETA WHENEVER PSI (IMAX, J+1-).G.DELTA(IMAX,J) PSI(IMAX,J)=DELTA(IMAX,J) OTHERWISE TRANSFER TO- SECOND END OF CONDITIONAL ENDPSI CONTINUE END OF CONDITT'ONAL BBSKIP=IB SECOND TAU=TAU+DTAU PRINT RESULTS TAU RCOMPUTE OPTIMUM VALUE OF A EXECUTE.AOPT. WHENEVER B85,TRANSFER TO FIRST WHENEVER BBFI.NE,TRANSFER TO START WHENEVER BZEROU,TRANSFER TO ENDU THROUGH UEND,FOR J=1,19J.E*JMAX UEND U(IMAX,J )OIDXI( IMAX)*(PSI(IMAX,J+1 )-PSI(IMAX,J-1) )/DZET.A2 1*RE.P. 1.5 ENDU CONT I NUE

-85THROUGH PSITOP,FOR I =1,1,IG I MAX PSITOP PSI (I,JMAX) =(4.0*PSI (I,LJMAX)-PSI ( I,JMAX-2) )/30 WHENEVER BZEROW.AND. bZEROU, TRANSFER TO SOLID THROUGH PSIEND,FOR VALUES OF I=IMAX THROUGH PSIEND,FOR J=1,1,J.E.JMAX PSIEND PSI ( I,J) =O.5*(5.O*PSI ( I-1,J)-40*PSI ( I-2,J ) 1+PSI(I-3,J)+(3.0*PSI(I,J)-4.0*PSI (I-1,J)+PSI (I-2,J))*0.5/I) SOL ID CONTINUE THROUGH GWTOP,FOR VALUES OF J=JMAX THROUGH GWTOP,FOR I=1,1,I.E.IMAX WHENEVER IE.1 W( 1,J ) =-OIDXI ( I ) *(-2.O*PSI (O,J )-3.0*PSI ( 1,J)+6.0*PSI ( 2,J)-PSI 1 (3,J ) ) / (6. O*DX I ) 2*RE.P.1*5 OR' WHENEVER I.Gol. AND.I. L.LIMAX W ( I,'J) =-OI DX I ( I) *-( P S I.( I -2, J) -8. O * P S I ( I - l, J ) + 8. O * PS I ( I +, J ) - P S 1I(I+2,J))/DXI12 2*RE.P.1.5 OTHERWISE W(' I,J)=-OIDXI('I)*(2.0*PSI('I+1,J)+3.0O*PSI(I, J)-6.0*PSI(I-1,J)+ 1PSI(I-2,J))/(6.0*DXI) 2*RE.P..5 END OF CONDITIONAL WHENEVER GAMTOP GAMMA (I, J)=( (PSI (I+1,J)-2. O*PSI (I,J)+PSI(I-1,J))/IDXICU(I)1(PSI (I+1,J)-PSI (I-1,J))*DXI2I(I)/IDXI(I)+O..5*DZSQI*( 28.0*PSI(I,J-1)-7.0*PSI(I,J)-PSI(I,J-2)) 3*OIDXI (I ) 4)*RE.'P.1.5 END OF CONDITIONAL GWTOP CONTIN UE THROUGH CORNER,FOR VALUES OF I=IMAX THROUGH CORNER,FOR VALUES OF J=JMAX WHENEVER GAMTOP GAMMA(IJ)=((2,*PSI(IJ)-50*PSI (tJ-1)+4.0*PSI(I.,J-2)-PSI(I, 1J-3) )*OIDXI (I)*DZSQI+ 2 0.5*(PSI (I,J)-PSI (I-1,J)-PSI(I I-2,J)+PSI -3,J) )/ 3IDXICU(I)-DXI2I(I)*(2 0*PSI(IJ)-PSI(I-,'J)-2.0*PSI(I-2,J)+PS 4I (I-3,J) )*OIDXI (I) 5 )*RE.P.1.5 END OF CONDITIONAL CORNER CONTINUE START CONTINUE R COMPUTE NEW U AND W VALUES AND NEW GAMMA(I,O) (I.E. DISC) VA RLUES EXECUTE SUBUWG. WHENEVER BBFINETRANSFER TO EXVEE WHENEVER BZEROW,TRANSFER TO WENDO THROUGH WEND,FORJ=1,1,J.G.JMAX WEND W(IMAX,J)=W(LIMAXJ) WENDO CONT INUE THROUGH NGAM,FOR VALUES OF I=IMAX THROUGH NGAM,FOR J= 1,1, J. E.JMAX NGAM GAMMA( I,J)=OIDXI (I )*(RE.P.1.5)*(PSI ( I J+1 )-2*0*PSI (I J)+ PSI I IJ-1) )*DZSQI WHENEVER BB2, TRANSFER TO NINTH R COMPUTE NEW V VALUES ~WHENEVER BB3, PRINT COMMENT$1 VSTAR(O,O)...VSTAR(IMAX,JMAX) 1 FOLLOWED BY VPRIME(O,O)...VPRIME(IMAX,JMAX)$

-86EXVEE EXECUTE SUBV. WHENEVER BBFINETRANSFER TO EXGA'M WHENEVER VXI,TRANSFER TO EXGAM WHENEVER BZEROV,TRANSFER TO EXGAM THROUGH VEND, FOR VALUES OF I=IMAX THROUGH VEND,'FOR J=1,I,J.E.JMAX VEhD V( I,J)=(3.0*DXI*OIDXI ( I )*V( I,J)+9.0 *V( I-l,)-45*V( I-2,J)+ 1V(I-3,J) )/5.5 R COMPUTE NEW GAMMA VALUES WHENEVER BB1, PRINT COMMENT$1 GAMMAS(O,O)...GAMMAS (IMAX,JMAX) 1 FOLLOWED BY GAMMAP(,00)...oGAMMAP(IMAX,JMAX)$ EXGAM EXECUTE SUBG. WHENEVER BB9 THROUGH SIXTH, FOR ITER=ITONE,DIT, ITER.G.DITMAX WHENEVER TAU.G. ITER*DTAU-EPSI.AND. TAU.L. ITER*DTAU+EPSI PRINT RESULTS TAU,DTAUW,'DXI,DZETA PRINT FORMAT ANSW3 WHENEVER BBSLCT THROUGH TENTHFOR VALUES OF I=IONE,ITWO, ITHREE,IFOUR THROUGH TENTHFOR VALUES OF J=JONE,JTWO,JTHREE,JFOUR TENTH PRINT FORMAT ANSW1,,I,J,PSI(I,J),U(I',J),'V(I, J,),W(I,J).,GAMMA(I, 1J) END OF CONDITIONAL WHENEVER BBFIX THROUGH THIRD, FOR I=ISTART,ISPACE,I.G.IEND THROUGH THIRD, FOR J=JSTAR-T,JSPACE,J.G.JJEND THIRD PRINT FORMAT ANSW1,I,J, PSI(I,J), U(IJ), V(I,J), W(I,J), GAM 1MA( I,J) END OF CONDITIONAL TRANSFER TO SEVEN END OF'CONDITIONAL SIXTH CONT'INUE END OF CONDITIONAL SEVEN WHENEVER TAU.L.TAUMAX, TRANSFER TO SECOND WHENEVER BBMAIN PRI.NT COMMENT$1 STEADY STATE VALUES $ PRINT RESULTS TAU,DTAU,DXIDZETA PRINT RESULTS BETA WHENEVER BPUNCH PUNCH FORMAT ENDBET, BETA PUNCH FORMAT ENDTAU-,TAU-,DTAU PUNCH FORMAT SPACE,DXI,DZETA THROUGH PUNCH,FOR I O, I.-G.IMAX THROUGH PUNCH-,FOR J=O,1,J.G.JMAX PUNCH PUNCH FORMAT ANSW2,Iq,J,PSI(I9,J), U(I,J), V(I,J), W(IJ)g GAMM 1A( I,J) END OF CONDITION'AL PRINT- FORMAT ANSW3 THROUGH FINAL9FORI=O l I.G.IMAX THROUGH FINAL, FOR J=0,1,J.G.JMAX FINAL PRINT FORMAT ANSWi,I,J, PSI(I,J), U(I,J), V(IJ), W(I,J), GAM. 1MA(I,J) END OF CONDITIONAL WHENEVER *ABS.BETA.G. BETMAX, TRANSFER TO FIFTH BETA=BETA+DBETA BBSKIP=OB TRANSFER TO FOURTH'FIFTH TRANSFER TO FIRST END OF PROGRAM

-87$COMPILE MAD, PRINT OBJECT, PUNCH OBJECT- COEFF001 DIMENSION RZERO(100),A(2000,DIM),BETA(100),.BETA1(100),BETA2( 1100) SMALLA(100), DUMMY(100l),DIM(3) INTEGER I,JiM,.N VECTOR VALUES TITLE=$S3,1HI,S5,6HSMALLA*$ VECTO.R VALUES ANSWER=$52.,I2,E13.6*$ START READ AND PRINT DATA N=DIM(2) M=DIM(2) THROUGH NEWR,FOR I=1,-1,I.G.N NEWR RZ.ERO(I)=RZERO(I)+DELR*(BETA('I)-BETA1(I))/(BETA2(I)-BETA1(I)) THROUGH ONE,FOR I=1,1,I.G.N THROUGH ONE FOR.J= 1,1 J..GN ONE A(I,J)=RZERO('I').P.J R=SLE.(N,M,A(1,1 ),SMALLA (),BETA(l),DUMMY(1), 0) PRINT RESULTS R PRINT RESULTS MN,RZERO(1)...RZERO(N),BETA(1)...BETA(N) PRINT FORMAT TITLE THROUGH ENDs FOR J=l,l,J.'G.N END PRINT FORMAT ANSWER,J,SMALLA(J) TRANSFER TO START END OF PROGRAM ~$COMPILE MAD, PRI.NT OBJECT,'PUNCH OBJECT EXPONOO1 INTEGER I,IMAX DIMENSION RZERO(100OOBETA(100 ),'BETA1(100),BETA2(100),ALPHA(10 10),C (100) VECTOR VALUES VALUE=$(S2,I3,S4,4E16.6)*$ VECTOR VALUES LABEL=$1H6,S.3,1HIS13,4HBETA,S115 HRZERO,S10,5H 1ALPHA,S13,1HC//*$ VECTOR VALUES TITLE=$S46-,40HCOMPUTED VALUES FOR THE EXPONENT 1OF BETA//*$ READ AND PRINT DATA THROUGH NEWR,F-OR I=,1I,-I.G.IMAX NEWR RZERO(I)=RZERO(I)+DELR*(BETA( I )-BETA( I )/(BETA2(I)-BETA1( I)) THROUGH ALPHAC,FOR I=1,1,I.E.IMAX ALPHA(I)'=(ELOG.(-RZERO(I+1))-ELOG.(RZERO(I)))/(ELOG.(BETA(I+1') 1)-ELOG..(BETA(I))) ALPHAC C(I)=RZERO(I)*(BETA(I).P.-ALPHA(I)) PRINT FORMAT TITLE PRINT FORMAT LABEL THROUGH ANSWER,FOR I=1,.1,I.E.IMAX ANSWER PRINT FORMAT VALUE, I,BETA(I),RZERO(I),ALPHA(I ),C(I ) END OF PROGRAM $ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT

-88RCOMPUTE INITIAL VALUES OF PSI, V, AND GAMMA EXTERNAL FUNCTION PROGRAM COMMON U,V,VSTAR,VPRIME,.W,GAMMAVSTILL RE, PSI,AAA 1,CEEDEE,,CAPB,CA.PCA,DXI2I,XISQIDXDZ,IDXICU,OVER2I,MINUSI,PL2US2I,OIDXI,OIDXIS, H,G,A, IDXI BB,DTAUD2,DXISQDXISQ2,DZE 3TSQDZETA2,DZETS2,DZSQI,DZETS4,PL2ARE,AOARE,DZETA6,DZET12,DXI 412,DXISQ1,DX DXISQHARE,EPSI,BETA,DXI,DZETA, I,J,ITMAXITER 5,BBJ1,BB2,BB3,BBF INE I START,JSTART, I ENDJEND 6,BB4,BB5,BB6,BB7,8B8 BB 9 7,DIM,ANSW1,ANSW2,ANSW3,BETMAX,DBETA,DTAU,OLDATA,TAU,TAUMAXD I 8TMAX, K,BEE, IMAX, JMAX, LIMAX,LJMAX, ENDBET 9,BBPSI, ENDTAU,SPACE DIMENSION U(1900,DIM) V(1900,DIM),VSTAR(1900,DIM),VPRIME(1900 1,DIM), W( 1900,DIM),GAMMA( 1900,DIM),GAMMAS( 1900,DIM),GAMMAP( 190 20,DIM) PSI (1900,DIM),DELTA(1900,DIM)AAA(250),CEE(250),DEE(25 30),CAPB(250),CAPC (250),DXI2I (250),IDXISQ(250),IDXDZ(250), IDXI 4CU(250), OVER2I(250),MINUSI.(250),PLUS2I(250).,OIDXI(250-),OIDXIS 5(250),H(50),G(50)50),A(11)IDXI(250) 6,DIM(3) EQIVALENCE (VSTAR (O),GAMMAS(O) )', (VPRIME(O),GAMMAP (O)), l(GAMMAP(O),DELTA(C0)) INTEGER IMAX,JMAX,DITMAX, LIMAX, LJMAX INTEGER ISTART, JSTART, IEND, JEND INTEGER.-I,J,K,ITMAX,ITER VECT"OR VALUES SPACE=$5HDXI= -E13.6,7HDZETA= E13.6*$ VECTOR VALUES ENDTAU=$5HTAU= E13.6,6HDTAU -E13.6*$ VECTOR VALUES OLDATA=$ (213,5E13.6 ) *$ VECTOR VALUES ENDBET=$6HBETA= Ei3.6*$ VECTOR VALUES ANSW1=$(S2, I3,S1,I3,5El6.6)*$ VECTOR VALUES ANSW2=$(2I3,5E13.6)*$ VECTOR -VALUES ANSW3=$1H6,S3,1HI,S3,1HJS10,3HPSI,S13,1HU,S14, 11HV,S15,1HW,S14,5HGAMMA//*$ BOOLEAN BB.,BB1,BB2,BB3,BB4,BB5,BB6',BB7,BB8,BB9 BOOLEAN BBPSI BOOLEAN VSTILL BOOLEAN BBFINE ENTRY TO BEGIN. THROUGH ONE,FOR J=O,1,J.G.JMAX V( 0 J)=O.O U(O,J)=O.O W(IMAX,J)=-H(J) W(O,J)=O0O U( IMAX,J)=0.O GAMMA(O,J)=O0O ONE PSI(OJ)=0.O THROUGH TWO,FOR I=O,1,I-.G.IMAX V( I,jMAX')=OO GAMMA(I,JMAX)=0. PSI(I,JMAX)=IDXISQ(I)*H(JMAX)/2.0 V(I,O)=IDXI(I) WHENEVER VSTILL,V(I,0)=0.0 PSI (: IO)=O.-O U( I 0)=0.0. U(I,JMAX)=0.0O W(I,JMAX)=-0.886 TWO W(I,0)=0.0 THROUGH THREE, FOR J=i,1,J.E.JMAX THROUGH THREE,FOR I=1,1,I.G.IMAX PSI(I,J)=IDXISQ(I)*H(J)/2.0

-89THREE V(I,J)=IDXI (I )*G(J) THROUGH FOUR, FOR J=1,1,J.E.JMAX THROUGH FOUR, FOR I=1,1,I.E.IMAX FOUR GAMMA(I,J)=(PSI'( I+1,J)-2.0*PSI(I,J)+PSI(I-1,J))/IDXICU(I) 1-(PSI(I 2+1,J)-PSI(I-1,J))*DXI2I (I)/IDXI (I )+(PSI (IJ+1)-2.0*PSI (IJ)+P 3SI(I,J-1))*DZSQI/IDXI (I) THROUGH UIMAX,FOR VALUES OF I=IMAX THROUGH UIMAX FOR J= 1 J. E.JMAX WHENEVER J.E.1 U( I 1 )=OIDXI ( I)*(-2.0*PSI ( I,O)-3.0*PSI ( I1 )+6.0*PSI ( I,2)-PSI ( 1 I,3 )) /DZETA6 OR WHENEVER J.G.1.AND.J.L.LJMAX U( I,J)=OIDXI (I )*(PSI (I,J-2)-8.0*PSI (I,J-1)+8.0*PSI (I J+)-PSI 1 ( I,J+2) )/DZET12 OTHERW I SE U( I,J)=OIDXI (I )*(2.0*PSI (I,J+1 )+3.0*PSI (I,J)-6.0*PSI (I,J-1 )+P 1SI (I,J-2))/DZETA6 U IMAX END OF CONDIT-IONAL WHENEVER BB2,.PRINT RESULTS PSI(0,O)...PSI (IMAX,JMAX.) i,V(O,O)...V(IMAX,JMAX) FUNCTION RETURN END OF FUNCTION $COMPILE MAD, PRINT OBJECT, PUNCH OBJECT AOPT.0O1 R COMPUTE A OPTIMt UM EXTERNAL FUNCTION PROGRAM COMMON U,V,VSTAR,VPRIME,W,GAMMA, VSTILL,RE, PSI,AAA 1,CEEDEE,DEECAPB,CAPC,DX 12 I,I DX ISQ,I DXDZ, IDXI CU,OVER2I,MINUSI,PL 2US2I,OIDXI,OIDXIS, H,G,A,IDXI,BB,DTAUD2,DXISQ,DXISQ2,DZE 3TSQ', DZETA2,DZETS2,DZSQI',DZETS4-,PL2ARE,AOARE, DZET.A6,.DZET12,DX I 412, DX.ISQ1, DXI2 DXISQH, ARE, EPSI, BETA,DX I.DZETA,I,J, ITMAX, ITER 5,BB1,BB2,B3,BBFI NE,ISTARTJ-ST-ART, I END,JEND 6,BB4,BB5,BB6,BB7,BB8,BB.9 7,DIM,ANSW1, ANNSW2,ANSW3,BETAX,DBETA,DTAU,OLDATA,TAU-, TAUMAX,D I 8TMAX,K,BEE,IMAX,JMAX,LIMAX,LJMAX,ENDBET 9,BBPSI,ENDTAU,SPACE DIMENSION U(1900,DIM ) V("1900,DIM),VSTAR(1900,DIM),VPRIME ( 1900 1,DIM),W( 1900,DIM),GAAMA( 1900 IM),GAMMAS( 1900,D IM),GAMMAP( 190 20,DIM),PSI ( 1900DIM),DE'LTA( 1900,DIM),AAA ( 250),CEE(250),DEE(25 30 —),CAPB(250),CAPC(250),DXI2I(-250 )-,IDXISQ(250') I DXDZ(250),IDXI 4CU(250),OVER2I(250),MINUSI (250) PLUS2I (250),OIDXI (250),OIDXIS 5(250),H(50), G(50),A(11'),IDXI(250) 6,DIM( 3) EQUIVALENCE (VSTAR(O),GAMMAS( O ) ), ( VPRIME ( O. ),GAMMAP(O) ), 1(GAMMAP(O),DELTA(O))INTEGER IMAX,JMAX DIT MAX,'LIMAX, LJMAX INTEGER I,J';, (t TMAX, I TER INTEGER I.START, JSTART, IEND, JEND BOOLEAN VSTILL 3OOLEAN BBFINE VECTOR VALUES ENDTAU=$5HTAU= E13.-.6,6HDTAU= E 13.6*$ VECTOR VALUES SPACE=$SHDXI: E13.6,7HDZETA= E13.6*$ VECTOR VALUES OLDATA=$( 21I3 5E13.6 ) *$ VECTOR VALUES ENDBET=$6HBETA= E13..6*$

-90VECTOR VALUES ANSW1=$(S2,I3,Sl1I3,5E16.6)*$ VECTOR VALUES ANSW2=$(2I3,5E13.6)*$ VECTOR VALUES ANSW3=$1H6,S3,1HI,S3,1HJ,S1O,3HPSI,S13,1HU,S14, 1 1HV,S15,1HW,S14,5HGAMMA//*$ BOOLEAN BB,BB1,BB2,BB3,BB4,BB5,BB6,BB7,BB8,BB9 BOOLEAN BBPSI ENTRY TO AOPT. K=O AOARE=A(K)/PL2ARE THROUGH KZERO, FOR ITER=1,1,ITERoG.ITMAX THROUGH KZERO,FORJ=1,,J.E.JMAX THROUGH KZERO, FOR I=l,1,I.G.LIMAX DELTA(I,J)=AOARE*(-IDXICU(I)*GAMMA(IJ)*RE.P.,-1.5+PSI(I+1,J) 1*MINUSI(I)+ 2PSI:(I-l,J)*PLUS2I(I)+ARE*(PSI(I,J+1)+PSI(I,J-1))-PL2ARE*PSI(I 3,J)) KZERO PSI (I,J)=PSI(I,J)+DELTA(I,J) FUNCTION RETURN END OF FUNCTION $'COMPILE MAD, PRINT OBJECT, PUNCH OBJECT AOPT.101 R COMPUTE A OPTIMUM EXTERNAL FUNCTION PROGRAM COMMON U,V,VSTAR,VPRIME,W,GAMMA,VSTI'LLs RE, PSI,AAA 1,CEE,DEE,CAPB,CAPC,DXI2I',IDXISQ,IDXDZ,IDXICU,OVER2'I,MINUSI,PL 2US2I,OIDXI,'OIDXIS, H,G,A, IDXI,BB,DTAUD2,DXISQ,DXISQ2,DZE 3TSQ,DZETA2,DZETS2,DZSQI,DZETS4, PL2ARE,AOARE, DZETA6,DZET12,DXI 412,DXISQ1,DXI2,DXISQH,ARE,EPSI,.BETA,DXI,DZETA,I,J,ITMAX,ITER 5,BB1,BB2,BB3,BBFINEISTA JS,START,IEND,JEND 6,BB4,BB5,BB6,BB7,BB8,BB9 7,DIM.,ANSW1,ANSW2,ANSW3,BETMAX,DBETA,DTAU,OLDATA,TAU,TAUMAX,DI 8TMAX,K, BEE, IMAXMAX, MAXI MAX,.LJ MAX, ENDBET 9,BBPSI,ENDTAUSPACE DIMENSION U(1900,DIM), V(1900,DIM),VSTAR(1900,DIM),VPRIME(1900 1,DIM),W(1900,DIM),GAMMA(1900,DIM).,GAMMAS(1900,DIM),GAMMAP(190 20,DIM),PSI(1900,DIM),DELTA(1900,DIM).,AAA(250),'CEE(250),DEE(25 30),CAPB(250),CAPC(250),DXI2I(250), IDXISQ(250'),I.DXDZ(250'),IDXI 4CU( 250),OVER2I (.250),MINUSI (250),PLUS2I (250),OIDXI( 250),OIDXIS 5(250),H(50)G(50 ),A(11), IDXI(250) 6,DIM(3) EQUIVALENCE (VSTAR(0),GAMMAS(O)),(VPRIME(0),GAMMAP(0)), l(GAMMAP(0),DELTA.('0)) INTEGER IMAX,JMAX,DITMAX, L IMAX,LJMAX INTEGER I,J,K,ITMAX,ITER INTEGER ISTART, JSTART, IEND,'JEND BOOLEAN BBFINE VECTOR VALUES ENDTAU=$5HTAU= E13.6,6HDTAU= E13.6*$ VECTOR VALUES SPACE=$5HDXI= E13.6,7HDZETA= E13.6*$ VECTOR VALUES OL.DATA=$(2I3,5E13 3.6)*$ VECTOR VALUES ENDBET=$6HBETA= E13.6*$ VECTOR VALUES ANSW1=$(S2,I3,S1,I3,5El6.6)*$ VECTOR VALUES ANSW2=$(2I3,5E13.6)*$ VECTOR VALUES ANSW3=$1H6,S3,1HI,S3,1HJ,S10,3HPSI,S13,1HU,S14, 1HV,S15,1HW,S 14, 5HGAMMA//*$

-91-,BOOLEAN BB, BB1,'BB2,BB3, BB4, BB5, BB6,BB7, BB8,BB9 BOOLEAN VST ILL BOOLEAN BBPS I ENTRY TO AOPT. THROUGH STORE, FOR I=O,1,I.G.IMAX THROUGH' STORE, FOR J=O,1,J.G.JMAX STORE VSTAR f I,J)=PSI (^I J) THROUGH AFTER, FOR K=O,1,K.G.10 AOARE=A(K) /PL2ARE THROUGH CHANGE, FOR I=O,1,IG. IMAX THROUGH CHANGE, FOR J=O,1,J.G.JMAX CHANGE PSI (I,J)=VSTAR(I,J) THROUGH NEXTIT. FOR ITER=l,lITER.G.ITMAX THROUGH NEXTX, FOR I=1,1,I.G.LIMAX THROUGH NEXTX, FOR J1,1l,J.E.JMAX DELTA(I,J)=AOARE*(-IDXICU(I)*GAMMA(I,J)*RE.P.-1.5+PSI (I+l,J) 1*MINUSI (I) + 2PSI ( I-1,J)*PLUS2I ( I )+ARE*(PSI ( I,J+1 )+PSI ( I,J-1) )-PL2ARE*PSI ( I 3,J)) NEXTX PS I ( I, J J)=PSI(I, J)+DELTA( I,J) NEXTIT WHENEVER BB4, PRINT RESULTS ITER,A(K),PSI(O,O)...PSI (IMAX, 1JMAX ),DELTA (0,0)...DELTA(IMAXJMAX) AFTER PRINT RESULTS A(K),PSI (O,0)..PSI (IMAX,JMAX ),DELTA(OO) *... 1DELTA( IMAXJMAX) WHENEVER BBPSI, PRI'NT RESULTS A(K),PSI (0,O)...PSI(IiMAX,JMAX), 1DELTA(O0,0)..DELTA(IMAX,JMAX) WHENEVER BB6 THROUGH MANYIT,. FOR I=O,,I.G. IMAX THROUGH MANYIT, FOR J=O,1,J.G.JMAX MANYIT PSI (I,J )=VSTAR( I,J) THROUGH. AFTERA, FOR VALUES OF K=5 AOARE=A (K) /PL2ARE ITMAX= I TMAX+DI TNIAX THROUGH. NEXTXA, FOR ITER=11, ITER.G.ITMAX THROUGH NEXTXA, FOR I=1,1,I.E.-IMAX THROUGH NEXTXA, FOR J=1,1,J.E.JMAX DELTA( I,J) =AOARE.*(-IDXICU( I )'GAMMA( I,J)*RE.P.-1.5+PSI (I+l,J) I*MINUSI(I)+ 2PSI (I-1,J)*PLUS2I ( I )+ARE*(PSI ( I,J+ )+PSI I,J-l ) )-PL.2ARE*PSI ( I 3,J)) NEXTXA PSI(I,J)=PSI(I,J)+DELTA(I-,J) AFTERA PRINT RESULTS A(K), PSI(-0,O )...PSI(IMAX JMAX),DELTA( O,O)... 1DELTA ( I MAX,MAX) END OF CONDITIONAL FUNCTION RETURN END OF FUNCTION $COMPILE MAD, PRINT OBJECT, PUNCH OBJECT SUBUW00] RCOMPUTE NEW U AND W VALUES, COMPUTE NEW GAMMA(I,O) VALUES EXTERNAL FUNCTION PROGRAM COMMON U,V,VSTAR,VPRIME,W,GAMMA,VSTILL, RE, PSI,AAA 1,CEE,DEE,CAPB,CAPC,DXI2I, I DXISQ,IDXDZ,IDXICU,OVER2 I,MINUS I,PL 2US2I,OIDXI,OIDXIS, H,G,A,I.DXI,BB,DTAUD2,DXISQ,DXISQ2,DZE 3TSQ,DZETA2,DZETS2,DZSQI,DZETS4,PL2ARE,AOARE,DZETA6,DZET12,DXI 412,DXISQ1,DXI2,DXISQH,ARE,EPSI,BETA,DXI,DZETA,I,J,ITMAX,ITER 5,BB1,BB2,BB3,BBFINE,ISTARTJSTART,IENDJEND 6,BB4,BB5,BB6,BB7,BB8,BB9

-927,D IM,ANSW1, ANSW2,ANSW3,BETMAX.DBETA,DTAU,OLDATA,TAU,TAUMAX,.DI 8TMAX,KBEE, I MAX, JMAX,LIMAX,LJMAX,ENDBET 9,BBPSI ENDTAU,SPACE DIMENSION U( 1900, DIM),V(190,DIM),VSTAR(1900,DIM),VPRIME(1900 1 DIM),W(1900, DIM), GAMMA (1900,DIM).,GAMM AS(1900,DIM),GAMMAP(190 20,DIM),PSI(1900,DIM),DELTA(1.900,DIM),AAA(250),CEE(250) DEE(25 30),CAPB(250),CAPC(250),DXI2I(250) IDXISQ(250) IDXDZ(250),IDXI 4CU( 250),OVER2I( 250)MINUSI.(250),PLUS2I (250),OIDXI (250),OIDXIS 5(250),H(50),G(50),A(11),IDXI(250) 6,DIM( 3) EQUIVALENCE (VSTAR(0 ),GAMMAS(O) ), (VPRIME(O) GAMMAP (O ) ), 1(GAMMAP(O),DELTA(O) ) INTEGER IMAX,JMAX,DI TMAX, LIMAX, LJMAX INTEGER ISTART, JSTART, IEND, JEND INTEGER I,J,K,ITMAX,ITER VECTOR VALUES SPACE=$5HDXI= E13.6,7HDZETA= E13.6*$ VECTOR VALUES ENDTAU=$5HTAU= E13.6,6HDTAU= E13.6*$ VECTOR VALUES OLDATA=$(2I13,5E13.6)*$ VECTOR VALUES ENDBET=$6HBETA= E13,6*$ VECTOR VALUES ANSW1=$(S2, I3,SlI3,5E16.6)*$ VECTOR VALUES ANSW2=$ (2I3, 5E13.6) *$' VECTOR VALUES ANSW3=$1H6,S3,1HI,53,1HJ,S10,3HPSI,S13,1HU,.S14, 11HV, S15,1HW, S14,5HGAMMA//*$ BOOLEAN VSTILL BOOLEAN BB,BB1, BB2 BB3BB3,BB4,B-B5, BB6,BB7,BB8,BB9 BOOLEAN BBFINE BOOLEAN BBPSI ENTRY TO SUBUWG. WHENEVER BBFINE,TRANSFER TO SKIP THROUGH ORIGIN, FOR I =1,1,I.G. IMAX ORIGIN GAMMA( I,O)=OIDXI (I )*( 8.0*PSI ( I.'1 )-7.0*PSI (I,O)-PSI (I 2) )/DZET 1S2 2*RE.P. 1.5 SKIP THROUGH NEXT1, FOR I=1,1,I.E.IMAX THROUGH NEXT1, FOR J=1:,1,J.E.-JMAX WHENEVER J.E.1 U(I,1)=OIDXI (I)*(-2.0*PSI (IO 0)-3.0*PSI(I,1)+6.O*PSI (I,2)-PSI( 1 I,3) ) /DZETA6 2*RE.P. 1.5 OR WHENEVER J.G. 1.AND.J.L.LJMAX U( I, J)=OIDXI(I )*(PS.I-(I,J-2)-8.0*PSI (.I,J-1.)+8.0*PSI (I. J+1)-PSI 1(I,J+2) ) /DZET12 2'*RE.P..1 5' OTHERWISE U( I,J)=OIDXI (I)*(2.0O*PSI (I,J+1)+3.0*PSI (I,J)-6.0*PSI (I J-1)+P 1SI (I,J-2) ) /DZETA6 2*RE. P.1.'5 END OF CONDITIONAL WHENEVER I.E.1 W(1,J)=-OIDXI(I)*(-2.0*PSI (0,J)-3.0*PSI(1,J)+6.0*PSI (2,J')-PSI 1(3,J))/(6.O*DXI) 2*RE.P. 1.5 OR WHENEVER I.G.1.AND. I.LoLIMAX W(IJ)=-OIDXI(I)*.(PSI (I-2,J)-8O0*PSI(I-1,J)+8.0*PSI(I+1IJ)-PS 1I(I+2,J) )/DXI12 2*RE.P.1.5 OTHERW I SE W( I J)=-OIDXI (I )*(2.0*PSI( I+I1,J)+3.O*PSI (I,J)-6.0*PSI ( I-1,J)+ n n T rS q J | i E t Pl rs r

-93 - 2*RE.P,1.5 NEXT1 END OF CONDITIONAL THROUGH WZERO,FOR J= 1 1, J GJMAX WZERO W(O,J)=(PSI (3J)-2.0*PSI (2,J)-11.0*PSI (1,J) )*DX ISQ1/5.0 1*RE.P.1.5 WHENEVER BB2, PRINT RESULTS U(0,O)..oU(IMAX,JMAX),W(OO),0) 1W(IMAX,JMAX),GAMMA(O,0)...GAMMA(IMAX,JMAX) WHENEVER TAU.G.0.0 WHENEVER BB7 PRINT RESULTS U(O,0) o.U(IMAX,JMAX),W(0,0)..W( IMAX,JMAX) END OF CONDITIONAL END OF.CONDITIO'NAL FUNCTION RETURN END OF FUNCTION $ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT SUBV.001' RCOMPUTE VSTAR AND VPRIME EXTERNAL FUNCTION PROGRAM COMMON.U,V,VSTAR,VPRIME,W,GAMMA,VSTILL, RE, PSI,AAA 1,CEE,DEE,CAPB,CAPC,DXI2I I-DXISQ, IDXDZ, IDXICU,OVER2I,MINUSI,PL 2US2I,OIDXI,OIDXIS, H,G,A,IDXI,BB,DTAUD2,DXISQ,DXISQ2,DZE 3TSQ,DZETA2,DZETS2,DZSQI,DZETS4,PL2ARE,AOARE,DZETA6,DZET1 2,DXI 412,DXISQ1,DXI2,DXISQH,AREEPSI,BETA,DXIDZETA,I,J, ITMAX ITER 5,BB1,BB2BB3,BBF INE, ISTART,JSTART, I END,JEND 6,BB4,BB5,BB6,BB7,BB8,BB9 7,DI M,ANSW1, ANSW2,ANSW3,BETMAXDBET TA,'DTAU,OLDA'TA,TAU,TAUMAX,D I 8TMAX,K, BEE IMAX,JMAX,LIMAX, LJMAX,ENDBET 9,BBPS I, ENDTAU,SPACE DIMENSION U(1900,DIMI)V(1900,DIM),VSTAR(1900,DIM),VPRIME(1900 1,DIM), W( 1900, DIM ), GAMMA ( 1900 D I M),GAMMAS( 1900,DIM), GAMMAP( 190 20,DIM),PSI (1900DIM), DELTA(1900, DIM),AAA(250),CEE(250),DEE(25 30),CAPB(250),CAPC(250), DXI2I (250), IDXISQ(250),IDXDZ(250),IDXI 4CU(250),OVER2I (250), MINUSI (250),PLUS2I (250),OIDXI (250),OIDXIS 5(250),H(50),G(50),A( 11), IDXI.'(250) EQUDIM(3)VALENCE ( EQUIVALENCE (VSTAR(O),GAMMAS(O)),(VPRIME(O),GAMMAP(O)), 1 (GAMMAP(0),DELTA(O0) ) INTEGER IMAX,JMAX,D I TMA X D I MAX, LJMAX INTEGER I,J,K, ITMAX, ITER VECTOR VALUES SPACE=$5HDXI= E13.6,7HDZETA= E13.6*$ VECTOR VALUES ENDTAU=$5HTAU= E13.6,6HDTAU= E13.6*$ VECTOR VALUES OLDATA=$(2I3,5E13.6.)*$ VECTOR VALUES ENDBET=$6HBETA= E13.6*$ VECTOR VALUES ANSW1=$(S2, I 3,S1, I 3,5E16.6) *$ VECTOR VALUES ANSW2=$'(2I3, 5E13.6) *$ VECTOR VALUES ANSW3=$1H6,S3,1HI,S3,1HJ,S10,3HPSI,513,1HU,514, 1 1HV,S 15, 1HW,S14,5HGAMMA//*$ BOOLEAN VSTILL BOOLEAN BB,B2 BB,B2 4,BB5,BB6,BB7,BB8 BB9 INTEGER ISTART, JSTART, IEND, JEND BOOLEAN BBFINE BOOLEAN BBPSI ENTRY TO SUBV. THROUGH REPEAT, FOR I =0 1, I.G. IMAX THROUGH REPEAT, FOR J=O,1,J.G.JMAX VSTAR(IMAX,J)=V(IMAXJ)

-94VSTAR( I90)=V(ItO) VSTAR(, JMAX)=V(' I JMAX) REPEAT VSTAR(OJ) =V(OJ ) THROUGH VEE1, FOR J=1,,J.G.LJM.AX THROUGH VEE3, FOR I=11I.G.G.LIMAX WHENEVER I E 1 CAPB (1)=DTAUD2+DX I SQ2*RE CAPC ( 1)=(V(1,J-1)*(W(,JJ )/DZETA2+DZSQI*RE ) +V(,J ) * ( DTAUD21U(1,J )/DXI- ( DZETS4+DXISQ1)RE)+V ( 1J+1 ) * (-W (,J)/DZETA2+DZSQI 2*RE) )/CAPB ((1) CEE (1 ) =U( 1,J )/DX I2-( DXISQI+DX I SQH) *RE OR WHENEVER I.L.LIMAX.AND. I.G.1 AAA(I)=-U(I J)/DXI2-('DXISQ1-DXI2I (I))*RE CEE( I )=U( I J)/DXI2-(DXISQ+DXI2I (I))*RE DEE(I)=V( I,J-1)*(W(I, I J)/DZETA2+DZSQI*RE)+V( I J)*(DTAUD21U( I J)/IDXI (I )-( D-ZETS4+1./IDXISQ( I) )*RE)+V( I,J+1)*(-W( I,J)/ 2DZETA2+DZSQI*RE) OTHERWISE AAA(I)=-U( I,J)/DXI2-(DXISQ1-DXI2I (I ) )*RE DEE ( I)=-V ( I + J)*UI, )(UIJ ) /D'X I 2- ( DXISQ1+DX II 2 I ) ) *RE ) +V( I,J-l ) * 1(W(I,J)/DZETA2+DZSQI*-RE)+V(I,',J)*(DTAUD2-U(I,J)/IDXI(I)-(DZETS 24+o.0/IDXISQ(-I.))*RE)+V(I,J+1)*(-W(I,J)/DZETA2+DZSQI*RE) END OF CONDITIONAL BEE=DTA-UD2+DX I SQ2*RE WHENEVER I.NE.1 CAPB(I)=BEE-AAA(I)CEE (I-1)/CAPB (I- 1) CAPC(I) =(DEE (I)-AAA (I)*CAPC(I-1) )/CAPB(I) VEE3 END OF CONDITIONAL THROUGH VEE1, FOR I=LIMAX,-1,I.L.1 WHENEVER I E. LIMAX,VSTAR(I J)=CAPC ( I ) VEE1 WHENEV-ER I. L.LIMAX,VSTAR ( I, J) =CAPC ( I )-(CEE( I )*VSTAR( I+1, J ) ) / 1CAPB(I) WHENEVER BB3, PRINT RESULTS VSTAR(O,O)..VSTAR(IMAX,JMAX) THROUGH REPEET, FOR I=O,1,I.G*IMAX THROUGH REPEET, FOR J=O,19J.G.JMAX VPRIME( I,0)=V( I,0) VPRIME( I,JMAX)=V( I,JMAX) VPR I ME (O,J ) =V( 0 J ) REPEET VPRIME(IMAX,J)=VSTAR(IMAX,.J) THROUGH VEE2, FOR I=1,1,.I.G.LIMAX THROUGH VEE4, FOR J=1,1,J.G.LJMAX WHENEVER J.E.1 CAPB(J)=DTAUD2+DZETS4*RE CAPC(1)=(-V(ItO)*(-W(I,I)/DZETA2-DZSQI*RE)+VSTAR(I-1,1)*(U(I, 11)/DXI2+(DXISQ1-DXI2I(I))*RE)+VSTAR(I,1)*(DTAUD2-DXISQ2*RE) 2+VSTAR(I+I,1)*(-U(I, 1)/DXI2+(DXISQ1+DXI2I(I)) *RE))/CAPB(1) CEE(1)=W (,11)/DZETA2-DZSQI*RE OR WHENEVER JL.L.LJMAX.AND. J.G.1 AAA(J)=-W(I,J)/DZETA2-DZSQI*RE CEE(J)=W(IJ)/DZETA2-DZSQI *RE DEE(J)=VSTAR(I-1,J)*(U(I,J)/DXI2+(DXISQ1-DXI2I(I))*RE) 1+VSTAR( I,J-)*(DTAUD2-RE*D 2XISQ2)+VSTAR(I+1,J)*(-U( I,J)/DXI2+(DXISQ1+DXI2I(I))*RE)+V(IJ 3)*(-U(IJ)/IIDXI.(I)-RE/I-DXISQ(I)) OTHERWISE AAA(J)=-W(I,J)/DZETA2-DZSQI*RE DEE (J) =-V( I,J+ 1 )*(W( I J)/DZETA2-DZSQ I*RE) +VSTAR ( I-1,J)*(U (.I,J 1)/DXI2+(DXISQl-DXI2I(I)) *RE)+VSTAR(I,J)*(DTAUD2-DXISQ2*RE+ 2VSTAR(I+1,J)*(-U(I,J)/DXI2+(DXISQ1+DXI2I(I))*RE)+V(-IJ)*(-U(I

-953,J)/IDXI (I)-OIDXIS(I) *RE) END OF CONDITIONAL BEE=DTAUD2+DZETS4*RE WHENEVER J.NE.1 CAPB (J) =BEE-AAA ( J )*CEE (J-1 )/CAPB.( J-1) CAPC(J)=(DEE(J)-AAA (J)CAPC (,J-1))/CAPB(J) VEE4 END OF CONDITIONAL THROUGH VEE2, FOR J=LJMAX,-1,J.L.1 WHENEVER J.E.LJMAX, VPRIME(IJ)=CAPC(J) VEE2 WHENEVER J.L.LJMAX, VPRIME(I,J)=CAPC(J)-CEE(J)VPRIIME(I,J+1)/ LCAPBE(J) WHENE.VER BB3, PRINT RESULTS VPRIME(O,O)...VPRIME(IMAX,JMAX) BB3=BB THROUGH VEEF, FOR I=0L9lI.Go.IMAX THROUGH VEEF, FOR J=OlJ.GoJiMAX VEEF V('I,J)=VPRIME( I,J) FUNCTION RETURN END OF FUNCTION $COMPILE MAD,PRINT OBJECT, PUNCH OBJECT SUBG.O01 RCOMPUTE GAMMAS AND' GAMMAP EXTERNAL FUNCTION PROGRAM COMMON U,V,VSTAR,VPRIME,W,GAMMA,VSTILL, RE, PSI,AAA 1,CEEDEE,CAPB,CAPCDXI2I9IDXISQ,IDXDZ,IDXICU,OVER2I,MINUSIPL 2US2I,OIDXI,OIDXIS, H,G,AIDXI,BB,DTAUD2,DXISQ,DXISQ2,DZE 3TSQDZETA2,DZETS2,DZSQI,DZ-TS4,PL2AREAOAREDZETA6,DZET12,DXI 412,DXISQ1,DXI2,DX2,DXISQH,ARE,EPSIBETADXI,DZETA,I,J,ITMAX,ITER 5,BB1,BB2,BB3,BBFINE, STARTJSTART,IEND,JEND 6,BB4,BB5,BB6,BB 7,BB89,BB9 7,DIM,ANSW1,ANSW2,ANSV39,BETMAXDBETADTAUOLDATA,TAU,TAUMAXDI 8TMAXqK,BEEIAXJMKAXLIMAXEELJAX IMAX,,ENDBET 9,BBP.S I,'NDTAU,SPACE DIMENSION U( 1900,DIM),V( 190'0,DIM),VSTAR( 1900-,DIM),VPRIME( 1900 1,DI),W( 1900 DIM),GAMMA( 1900, DIM),GAMMAS(1900 DIM),GAMAP( 190 20,DIM),PSI(1900,DIM),DELTA(1900',DIM),AAA(250),CEE(.250),DEE(25 30),CAPB(250),CAPC(250),DXI2I (250) IDXISQ(250), IDXDZ(250) IDXI 4CU(250),OVER2I(250),MINUSI (250),PLUS.2I (250) OIDXI (25.0),OIDXIS 5(250),H(50),G(50),A(11), IDXI(250) 6,DIM( 3) EQUIVALENCE (VSTAR (),GAMMAS(O) ), (VPRIME(0),GAMMAP (0 ) ), 1(GAMMAP(0),DELTA(0) ) INTEGER IMAX,JMAXDITMAXLIMAXLJMAX INTEGER IJ,K, ITMAX, ITER INTEGER ISTART, JSTART, IEND, JEN"D BOOLEAN BBFINE VECTOR VALUES SPACE=$SHDXI= E13.6,7HDZETA= E13.6*$ VECTOR VALUES ENDTAU=$5HTAU= E13.6,6HDTAU= E13.6*$ VECTOR VALUES OLDATA=$(2I3,5E13.6 )*$ VECTOR VALUES ENDBET=$6HBETA= E13.6*$ VECTOR VALUES ANSW1=$(S2,I3,Sl I3,5E1E6.6)*$ VECTOR VALUES ANSW2=$(2I3, 5E13 6) *$ VECTOR VALUES ANSW3=$1H6,S3,1HI,S3,lHJ,S10,3HPSI,513,1HU,S14, 11HV,S15,1HW,S 14,5HGA MMA// *$ BOOLEAN VSTILL BOOLEAN BBe, BB1,BB2, BB3, BB4,BB5 BB6, BB7,BBB,BB9

-96BOOLEAN BBPSI ENTRY TO SUBG. THROUGH- REPEAT, FOR I=OI,1,I.G. IMAX THROUGH REPEAT,,FOR J=O,1,J.G.JMAX GAMMAS( 0J) =GAMMA(OJ) GAMMAS ( I MAX, J ) =GAMMA ( IMAX, J ) GAMMAS(I sJMAX)=GAMMA( IJMAX) REPEAT GAMMAS(I90)=GAMMA(I O) THROUGH GEE1. FOR J=,119J.G.LJMAX THROUGH GEE39 FOR. I=ll,I.G.LIMAX WHENEVER IoE.1 CAPB ( 1) =DTA.UD2+DX I SQ2*RE CAPC ( 1) = (GAMMA ( 1,J-1 ) * ( W( 1,J ) /DZETA2+DZSQI *RE )+GAMMA ( 1,J ) *( 1DTAUD2+U ( 1 J /DX I-( DZETS4+DX I SQ1 )*RE) +GAMMA (1 J+1 ) * ( -W 1,J ) / 2DZETA2+DZSQI*RE)+(VPRIME(1,J)/IDXDZ(I))*(VPRIME(1lJ+1)-VPRIME 3(19J-1) ) )/CAPB( 1) CEE (1)=U (1 J) /DX I-( DXISQ1+DXI SQH) *E OR WHENEVER I.L.LIMAX.AND. I.G.1 AAA(I)=-U( I,J)/DXI2DXI2-DXISQ1-DXI2I (I))*RE CEE( I )=U( I J)/DX I2-(DXISQ1+DXI 2I (I )*RE DEE(I)=GAMMA(I,J-1)*(W(II,J)/DZETA2+DZSQI*RE)+GAMMA(I,J)*(DTAU 1D2+U( I,J)/IDXI (I)-(DZETS4+OIDXIS(I))*RE)+GAMMA(I,J+1)*(-W( I J 2)/DZETA2+DZSQI*RE)+(VPRIME( I,J)/IDXDZ( I) )*(VPRIME(I,J+1)-VPRI 3ME( I,J-1 ) ) OTHERWISE AAA(I)=-U( I,J)/DXI2-(DXISQ1-DXI2I (I))*RE DEE(I)=-GAMMA(I+1,J)*(U(I,J)/DXI2-(DXISQ1+DXI2I(I))*RE)+ 1GAMMA ( I,J-1')*(W ( I,J)/DZETA2+DZSQI*RE)+GAMMA ( IJ)*(DTAUD2+U( 2I,J)*OIDXI(I)-(DZETS4+OIDXIS(I))*RE)+GAMMA(I,J+1)*(-W(I,J)/ 3DZETA2+DZSQI*RE)+(VPRIME(I,J)/IDXDZ(I) )*(VPRIME(IJ+1) 4-VPRIME('I,J-1)) END OF CONDITIONAL BEE=DTAUD2+DXISQ2*RE WHENEVER I.NE.1 CAPB( I )=BEE-AAA( I )*CEE( I-1)/CAPB( I-1 ) CAPC (I)=(DEE(I)-AAA(I)*CAPC(I-1))/CAPB(I) GEE3 END OF CONDITIONAL THROUGH GEEl, FOR I=LIMAX,-lI.L.1 WHENEVER I.E.LIMAX GAMMAS ( I,J)=CAPC (I) OR WHENEVER I.L.LIMAX GAMMAS( I,J)=CAPC(I)-CEE( I )*GAMMAS( I+1,J)/CAPB(I) GEE1 END OF CONDITIONAL WHENEVER BB1, PRINT RESULTS GAMMAS(O,O)...GAMMAS(IMAXJMAX) THROUGH REPEET, FOR I=O,1,I.G.IMAX THROUGH REPEET, FOR J=O,1,J.G.JMAX GAMMAP (OJ) =GAMMAS (0,3J) GAMMAP(IO)=GAMMAS(I O) GAMMAP(IMAX,J)=GAMMAS(IMAX,J) REPEET GAMMAP(I,JMAX)=GAMMAS(I,JMAX) THROUGH GEE2, FOR I=1,1,I.G.LIMAX THROUGH GEE4, FOR J=ltl, J.G.LJMAX WHENEVER J.E.1 CAPB ( 1 ) =DTAUD2+DZETS4*RE CEE(1)=W(I,1)/DZETA2-DZSQI*RE CAPC(1)=(GAMMAS(I-1,1)*(U(I,1)/DXI2+(DXISQ1-DXI2I (I))*RE)+ 1GAMMAS(I, 1.)*(DTAUD2-DXI SQ2*RE)+GAMMAS(I+1,1 )*(-U(I, 1)/DXI2+( 2DXISQL+DXI2I (I))*RE)+GAMMA(I,1)*(U'(.I,1)/IDXI (I)-OIDXIS(I)*RE) 3+(V( I,1)/IDXDZ( I) )*(V( I,2)-V( I,-O) )+GAMMAP( I,O)*(DZSQI*RE+W( I,

-97-. 41)/DZETA2))/CAPB(1)' OR WHENEVER JoL.LJMAX.AND.J.G.1 AAA(J)=-W(I,J)/DZETA2-DZSQI*RE CEE(J)=W(I,J)/D-ZETA2-DZSQI*RE DEE(J)=GAMMAS(I-1,J)*(U(I,J)/DXI2+(DXISQ1-DXI2I(I))*RE)+ 1GAMMAS(I,J)*(DTAUD2-DXISQ2*RE)+GAMMAS(I+1,J)*(-U(I,J)/DXI2+( 2DXISQ1+DXI2I(I))*RE)+GAMMA( I,J)*(U( I,J)/IDXI(I )-OIDXIS( I )*RE) 3+(V(I,J)/IDXDZ(I))*(V(I,J+1)-V(I,J-1)) OTHERWISE CEE(J)=W(I,J)/DZETA2-DZSUI*RL AAA(J)=-W(I,J)/DZETA2-DZSQI*RE DEE(J)=GAMMAS(I-1,J)*(U'(I,J)/DXI2+(DXISQ1-DXI2I(I))*RE)+ 1GAMMAS( I,J)*(DTAUD2-DXISQ2*RE)+GAMMAS(I+1, J)*(-U ( I,J)/DXI2+( 2DXISQ1+DXI2I (I))*RE)+GAMMA( I,J)*(U( I,J)/IDXI(I)-OIDXIS( I )*RE) 3+(V (I,J)/IDXDZ(I))*(V(I,J+ )-V(I,J-1)) 4-GAMMA(I,J+1)*(WlI,J')/DZETA2-RE*DZSQI) END OF CONDITIONAL BEE=DTAUD2+DZETS4*RE W.HENEVER J.NE.1, CAPB(J)=BEE-AAA (J)*CEE(J-l)/CAPB (J-1) CAPC(J)=(DEE(J)-AAA(J)*CAPC(J-1))/CAPB(J) GEE4 END OF CONDITIONAL THROUGH GEE2, FOR J=LJMAX,-i,J.L.1 WHENEVER J.E.LJMAX GAMMAP( I,J)=CAPC(J) OR WHENEVER J.L.LJMAX GAMMAP J=CAPCJCEEJ*GAMMAP(I,J+1)/C.APB(J) GEE2 END OF CONDITIONAL WHENEVER B'B1, PRINT RESULTS GAMMAP(O,O)...GAMMAP(IMAX,JMAX) BB1=BB THROUGH GEEF, FOR I=O,1,I.G.IMAX THROUGH GEEF, FOR J=O,liJ..G.JMAX GEEF GAMMA(I,J)=GAMMAP(I,J) FUNCTION RETURN END OF FUNCTION

BIBLIOGRAPHY 1. Barua, S. N,, " A Source in a Rotating Fluid," Q, Journ. Mech, Appl. Math,, 8 (1955) 22. 2, Birkhoff, G,, Personal communication. 3. Birkhoff, G,, Varga, R, S., Young, D,, "Alternating Direction Implicit Methods," Advances in Computers, 3_, Academic Press Inc., New York, 1962. 4, Bddewadt, U, T., "Die Drehstre'mung {tber festen Grunde, ZM, 20, (1940) 241. 5. Cochran' WO Go, "The Flow Due to a Rotating Disc," Proc. Cambridge Phil, Soco, 30 (1934) 365o 6. Fromm, J, E,, Harlow, F, H,, "Numerical Solution of the Problem of Vortex Street Development," Physics of Fluids,, (July, 1963) 7o 7, Long, R, R,, "A Vortex in an Infinite Viscous Fluid," J. Fluid Mech,, 11, 1961, 8, Long, R, R., "Sources and Sink at the Axis of a Rotating Liquid," Q, Journ,o Mech, Appl. Math,, ~, 1956, 9. Rosenblat, S,, "Flow Between Torsionally Oscillating Discs," J, Fluid Mech,, 8 July, 1960, 10, Schlicting, H,, Boundary Layer Theory, McGraw-Hill, New York, 1960o 11, Stewartson, K,, Proc. Camb. Phil, Soce, L49 3335 12, Stewartson, K,, "On Almost Rigid Rotations," J. Fluid Mech,,, 1957o 135 Taylor, G, I,, "The Boundary Layer in the Converging Nozzle of a Swirl Atomizer," Q, Journ, Mech, Appl, Math., L, 1950, 14, Von Karman, T,, "Laminare und Turbulente Reibung," Z 1, 1, (1921) 233-252; NACA T, M,, (1946) 1092, 15e Wilkes, J. 0,, The Finite Difference Computation of Natural Convection in an Enclosed Rectangular Cavity, Ph.D. Thesis, The University of Michigan, 1963, -98